How the huge energy of quantum vacuum gravitates to drive the slow accelerating expansion of the Universe

TL;DR

This paper proposes a new perspective on quantum vacuum gravity, suggesting it leads to a naturally slow accelerating universe without the need for a cosmological constant or dark energy, potentially resolving the old cosmological constant problem.

Contribution

It introduces a modified gravitational behavior of quantum vacuum energy that results in a slow acceleration, challenging previous assumptions and eliminating the need for fine-tuned dark energy.

Findings

Quantum vacuum gravitates differently than previously thought.

The universe's expansion rate is naturally slow, proportional to an exponential decay involving the cutoff scale.

The model potentially resolves the old cosmological constant problem without fine-tuning.

Abstract

We investigate the gravitational property of the quantum vacuum by treating its large energy density predicted by quantum field theory seriously and assuming that it does gravitate to obey the equivalence principle of general relativity. We find that the quantum vacuum would gravitate differently from what people previously thought. The consequence of this difference is an accelerating universe with a small Hubble expansion rate $H\propto \Lambda e^{-\beta\sqrt{G}\Lambda}\to 0$ instead of the previous prediction $H=\sqrt{8\pi G\rho^{vac}/3}\propto\sqrt{G}\Lambda^2\to\infty$ which was unbounded, as the high energy cutoff $\Lambda$ is taken to infinity. In this sense, at least the "old" cosmological constant problem would be resolved. Moreover, it gives the observed slow rate of the accelerating expansion as $\Lambda$ is taken to be some large value of the order of Planck energy or higher. This result suggests that there is no necessity to introduce the cosmological constant, which is required to be fine tuned to an accuracy of $10^{-120}$, or other forms of dark energy, which are required to have peculiar negative pressure, to explain the observed accelerating expansion of the Universe.

Equations

$$G_{\mu \nu }+{\lambda }_b{g}_{\mu \nu }=8\pi G{T}_{\mu \nu }^{\mathrm{vac}},$$

$$T_{\mu \nu }^{\mathrm{vac}}(t,\mathbf{x})=-{\rho }^{\mathrm{vac}}{\eta }_{\mu \nu }.$$

$$T_{\mu \nu }^{\mathrm{vac}}(t,\mathbf{x})=-{\rho }^{\mathrm{vac}}{g}_{\mu \nu }(t,\mathbf{x}).$$

$${\nabla }^{\mu }{T}_{\mu \nu }^{\mathrm{vac}}=0.$$

$$G_{\mu \nu }+{\lambda }_{\mathrm{eff}}{g}_{\mu \nu }=0,$$

$${\lambda }_{\mathrm{eff}}={\lambda }_b+8\pi G{\rho }^{\mathrm{vac}};$$

$$G_{\mu \nu }=-8\pi G{\rho }_{\mathrm{eff}}^{\mathrm{vac}}{g}_{\mu \nu },$$

$${\rho }_{\mathrm{eff}}^{\mathrm{vac}}={\rho }^{\mathrm{vac}}+\frac{{\lambda }_b}{8\pi G}.$$

$$d{s}^2=-d{t}^2+a^2(t)\left(d{x}^2+d{y}^2+d{z}^2\right).$$

$$3$$

$$\ddot{a}$$

$$a(t)=a(0)e^{Ht},$$

$${\lambda }_{\mathrm{eff}}=3{\Omega }_{\lambda }{H}_0^2\approx 4.32\times 10^{-84}(\mathrm{GeV}{)}^2,$$

$${\rho }_{\mathrm{eff}}^{\mathrm{vac}}={\Omega }_{\lambda }{\rho }_{\mathrm{crit}}\approx 2.57\times 10^{-47}(\mathrm{GeV}{)}^4,$$

$$\varphi (t,\mathbf{x})=\int \frac{d^{3}k}{(2\pi {)}^{3/2}}\frac{1}{\sqrt{2\omega }}\left(a_{\mathbf{k}}{e}^{-i(\omega t-\mathbf{k}\cdot \mathbf{x})}+a_{\mathbf{k}}^{†}{e}^{+i(\omega t-\mathbf{k}\cdot \mathbf{x})}\right),$$

$$a_{\mathbf{k}}|0⟩=0,\text{mbox}forall\mathbf{k},$$

$$H=\int d^{3}x{T}_{00}=\frac{1}{2}\int d^{3}k\omega \left(a_{\mathbf{k}}{a}_{\mathbf{k}}^{†}+a_{\mathbf{k}}^{†}{a}_{\mathbf{k}}\right),$$

$$T_{00}=\frac{1}{2}\left({\dot{\varphi }}^2+(\nabla \varphi {)}^2\right).$$

$$\displaystyle T_{00}(t,\mathbf{x})=\frac{1}{2}\int\frac{d^{3}kd^{3}k^{\prime}}{(2\pi)^{3}}\frac{1}{2}\left(\sqrt{|\mathbf{k}||\mathbf{k}^{\prime}|}+\frac{\mathbf{k}\cdot\mathbf{k}^{\prime}}{\sqrt{|\mathbf{k}||\mathbf{k}^{\prime}|}}\right)\Bigg{(}a_{\mathbf{k}}a_{\mathbf{k}^{\prime}}^{{\dagger}}e^{-i\left[\left(|\mathbf{k}|-|\mathbf{k}^{\prime}|\right)t-\left(\mathbf{k}-\mathbf{k}^{\prime}\right)\cdot\mathbf{x}\right]}$$

$$\displaystyle+a_{\mathbf{k}}^{{\dagger}}a_{\mathbf{k}^{\prime}}e^{+i\left[\left(|\mathbf{k}|-|\mathbf{k}^{\prime}|\right)t-\left(\mathbf{k}-\mathbf{k}^{\prime}\right)\cdot\mathbf{x}\right]}-a_{\mathbf{k}}a_{\mathbf{k}^{\prime}}e^{-i\left[\left(|\mathbf{k}|+|\mathbf{k}^{\prime}|\right)t-\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{x}\right]}-a_{\mathbf{k}}^{{\dagger}}a_{\mathbf{k}^{\prime}}^{{\dagger}}e^{+i\left[\left(|\mathbf{k}|+|\mathbf{k}^{\prime}|\right)t-\left(\mathbf{k}+\mathbf{k}^{\prime}\right)\cdot\mathbf{x}\right]}\Bigg{)},$$

$$\left\langle\big{(}T_{00}-\left\langle T_{00}\right\rangle\big{)}^{2}\right\rangle=\frac{2}{3}\langle T_{00}\rangle^{2},$$

$$⟨T_{00}⟩=\frac{{\Lambda }^4}{16{\pi }^2},$$

$$\Delta\rho^{2}\left(\Delta x\right)=\frac{\left\langle\left\{\big{(}T_{00}\left(t,\mathbf{x}\right)-T_{00}\left(t,\mathbf{x}^{\prime}\right)\big{)}^{2}\right\}\right\rangle}{\frac{4}{3}\left\langle T_{00}(t,\mathbf{x})\right\rangle^{2}},$$

$$d{s}^2=-d{t}^2+a^2(t,\mathbf{x})(d{x}^2+d{y}^2+d{z}^2).$$

$$G_{00}$$

$$G_{ii}$$

$$G_{0i}$$

$$G_{ij}$$

$$⟨T_{ij}⟩=0,i,j=1,2,3,i\ne j.$$

$$⟨T_{ij}^2⟩\sim {⟨T_{00}⟩}^2.$$

$$L(t)={\int }_{{\mathbf{x}}_1}^{{\mathbf{x}}_2}\sqrt{a^2(t,\mathbf{x})}dl$$

$$H(t)=\frac{\dot{L}}{L}=\frac{{\int }_{{\mathbf{x}}_1}^{{\mathbf{x}}_2}\frac{\dot{a}}{a}(t,\mathbf{x})\sqrt{a^2(t,\mathbf{x})}dl}{{\int }_{{\mathbf{x}}_1}^{{\mathbf{x}}_2}\sqrt{a^2(t,\mathbf{x})}dl},$$

$$H(t)=\frac{\dot{a}}{a}(t).$$

Full text

I Introduction

The two pillars that much of modern physics is based on are Quantum Mechanics (QM) and General Relativity (GR). QM is the most successful scientific theory in history, which has never been found to fail in repetitive experiments. GR is also a successful theory which has so far managed to survive every test lrr-2006-3 . In particular, the last major prediction of GR–the gravitational waves, has finally been directly detected on Sept 2015 Abbott:2016blz . However, these two theories seem to be incompatible at a fundamental level (see e.g. Wald:1984rg ). The unification of both theories is a big challenge to modern theoretical physicists.

While the test of the combination of QM and GR is still difficult in lab, our Universe already provides one of the biggest confrontations between both theories: the Cosmological Constant Problem RevModPhys.61.1 . Quantum field theory (QFT) predicts a huge vacuum energy density from various sources. Meanwhile, the equivalence principle of GR requires that every form of energy gravitates in the same way. When combining these concepts together, it is widely supposed that the vacuum energy gravitates as a cosmological constant. However, the observed effective cosmological constant is so small compared with the QFT’s prediction that an unknown bare cosmological constant (6) has to cancel this huge contribution from the vacuum to better than at least 50 to 120 decimal places! It is an extremely difficult fine-tune problem that gets even worse when the higher loop corrections are included Padilla:2015aaa .

In 1998, the discovery of the accelerating expansion of the Universe 1538-3881-116-3-1009 ; 0004-637X-517-2-565 has further strengthened the importance of this problem. Before this, one only needs to worry about the “old” cosmological constant problem of explaining why the effective cosmological constant is not large. Now, one also has to face the challenge of the “new” cosmological constant problem of explaining why it has such a specific small value from the observation, which is the same order of magnitude as the present mass density of the Universe (coincidence problem).

This problem is widely regarded as one of the major obstacles to further progress in fundamental physics (for example, see Witten 2001 witten ). Its importance has been emphasized by various authors from different aspects. For example, it has been described as a “veritable crisis” (Weinberg 1989, RevModPhys.61.1 p.1), an “unexplained puzzle” (Kolb and Turner 1993 , Kolb:EU90 p.198), “the most striking problem in contemporary physics” (Dolgov 1997 Dolgov:1997za p.1) and even “the mother of all physics problems” , “the worst prediction ever”(Susskind 2015 citeulike:430764 chapter two). While it might be possible that people working on a particular problem tend to emphasize or even exaggerate its importance, those authors all agree that this is a problem that needs to be solved, although there is little agreement on what is the right direction to find the solution pittphilsci398 .

In this paper, we make a proposal for addressing the cosmological constant problem. We treat the divergent vacuum energy density predicted by QFT seriously without trying renormalization and assume that it does gravitate to obey the equivalence principle of GR. We notice that the magnitude of the vacuum fluctuation itself also fluctuates, which leads to a constantly fluctuating and extremely inhomogeneous vacuum energy density. As a result, the quantum vacuum gravitates differently from a cosmological constant. Instead, at each spatial point, the spacetime sourced by the vacuum oscillates alternatively between expansion and contraction, and the phases of the oscillations at neighboring points are different. In this manner of vacuum gravitation, although the gravitational effect produced by the vacuum energy is still huge at sufficiently small scales (Planck scale), its effect at macroscopic scales is largely canceled. Moreover, due to the weak parametric resonance of those oscillations, the expansion outweighs contraction a little bit during each oscillation. This effect accumulates at sufficiently large scales (cosmological scale), resulting in an observable effect—the slow accelerating expansion of the Universe.

Our proposal harkens back to Wheeler’s spacetime foam Wheeler:1957mu ; misner1973gravitation and suggests that it is this foamy structure which leads to the cosmological constant we see today.

The paper is organized as follows: in section II, we first review several key steps in formulating the cosmological constant problem; in section III, we point out that the vacuum energy density is not a constant but is constantly fluctuating and extremely inhomogeneous; in section IV, we investigate the differences made by the extreme inhomogeneity of the quantum vacuum by introducing a simple model; in section V, we give the solutions to this model by solving the Einstein field equations and show how metric fluctuations leads to the slow accelerating expansion of the Universe; in section VI, we explain the meaning of our results; in section VII, we investigate the back reaction effect of the resulting spacetime on the matter fields propagating on it; in section VIII, we generalize our results to more general metrics; in section IX, we discuss some questions raised and a couple of new concepts suggested by the different way of vacuum gravitation.

The units and metric signature are set to be $c=\hbar=1$ and $(-,+,+,+)$ throughout except otherwise specified.

II The formulation of the cosmological constant problem

The cosmological constant problem arises when trying to combine GR and QFT to investigate the gravitational property of the vacuum:

[TABLE]

where $G_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}$ is the Einstein tensor and the parameter $\lambda_{b}$ is the bare cosmological constant.

One crucial step in formulating the cosmological constant problem is assuming that the vacuum energy density is equivalent to a cosmological constant. First, it is argued that the vacuum is Lorentz invariant and thus every observer would see the same vacuum. In Minkowski spacetime, $\eta_{\mu\nu}$ is the only Lorentz invariant $(0,2)$ tensor up to a constant. Thus the vacuum stress-energy tensor must be proportional to $\eta_{\mu\nu}$ (see, e.g. carroll2004spacetime , pittphilsci398 )

[TABLE]

This relation is then straightforwardly generalized to curved spacetime:

[TABLE]

If $T_{\mu\nu}^{\mathrm{vac}}$ does take the above form (3), the vacuum energy density $\rho^{\mathrm{vac}}$ has to be a constant, which is the requirement of the conservation of the stress-energy tensor

[TABLE]

The effect of a stress-energy tensor of the form (3) is equivalent to that of a cosmological constant, as can be seen by moving the term $8\pi GT_{\mu\nu}^{\mathrm{vac}}$ in (1) to the left-hand side

[TABLE]

where,

[TABLE]

Or equivalently by moving the term $\lambda_{b}g_{\mu\nu}$ in (1) to the right-hand side

[TABLE]

where,

[TABLE]

So anything that contributes to the energy density of the vacuum acts like a cosmological constant and thus contributes to the effective cosmological constant. Or equivalently we can say that the bare cosmological constant acts like a source of vacuum energy and thus contributes to the total effective vacuum energy density. This equivalence is the origin of the identification of the cosmological constant with the vacuum energy density.

Following the above formulations, the effective cosmological constant $\lambda_{\mathrm{eff}}$ or the total effective vacuum energy density $\rho_{\mathrm{eff}}^{\mathrm{vac}}$ are the quantities that can be constrained and measured through experiments. While solar system and galactic observations have placed a small upper bound on $\lambda_{\mathrm{eff}}$, large scale cosmological observations provide the most accurate measurement. It is interpreted as a form of “dark energy”, which drives the observed accelerating expansion of the Universe 1538-3881-116-3-1009 ; 0004-637X-517-2-565 .

Based on the assumption of homogeneity and isotropy of the Universe, the metric has the cosmology’s standard FLRW form, which is, for the spatially flat case,

[TABLE]

Then by applying the equations (5) or (7) for the above special metric (9), one obtains the contributions to the Hubble expansion rate $H=\dot{a}/a$ and the acceleration of the scale factor $\ddot{a}$ from $\lambda_{\mathrm{eff}}$ and/or $\rho_{\mathrm{eff}}^{\mathrm{vac}}$ are

[TABLE]

The solution to the dynamic equation (11) is

[TABLE]

where $H$ is determined by the initial value constraint equation (10).

According to the Lambda-CDM model of the big bang cosmology, the effective cosmological constant is responsible for the accelerating expansion of the Universe as shown in (11) and contributes about $69\%$ to the current Hubble expansion rate Ade:2015xua :

[TABLE]

or

[TABLE]

where $\Omega_{\lambda}=0.69$ is the dark energy density parameter, $H_{0}$ is the current observed Hubble constant and $\rho_{\mathrm{crit}}=\frac{3H_{0}^{2}}{8\pi G}$ is the critical density.

Unfortunately the predicted energy density of the vacuum from QFT is much larger than this. It receives contributions from various sources, including the zero point energies ($\sim 10^{72}(\mathrm{GeV})^{4}$) of all fundamental quantum fields due to vacuum fluctuations, the phase transitions due to the spontaneous symmetry breaking of electroweak theory ($\sim 10^{9}(\mathrm{GeV})^{4}$) and any other known and unknown phase transitions in the early Universe (e.g. from chiral symmetry breaking in QCD ($\sim 10^{-2}(\mathrm{GeV})^{4}$), grand unification ($\sim 10^{64}(\mathrm{GeV})^{4}$) etc)pittphilsci398 ; Carroll:2000fy . Each contribution is larger than the observed value (14) by $50$ to $120$ orders of magnitude. There is no mechanism in the standard model which suggests any relations between the individual contributions, so it is customary to assume that the total vacuum energy density is at least as large as any of the individual contributions pittphilsci398 . One thus has to fine tune the unknown bare cosmological constant $\lambda_{b}$ to a precision of at least 50 decimal places to cancel the excess vacuum energy density.

III The fluctuating quantum vacuum energy density

The vacuum energy density is treated as a constant in the usual formulation of the cosmological constant problem. While this is true for the expectation value, it is not true for the actual energy density.

That’s because the vacuum is not an eigenstate of the local energy density operator $T_{00}$, although it is an eigenstate of the global Hamiltonian operator $\operatorname{H}=\int d^{3}x\,T_{00}$. This implies that the total vacuum energy all over the space is constant but its density fluctuates at individual points.

To see this more clearly, consider a quantized real massless scalar field $\phi$ in Minkowski spacetime as an example:

[TABLE]

where the temporal frequency $\omega$ and the spatial frequency $\mathbf{k}$ in (15) are related to each other by $\omega=|\mathbf{k}|$.

The vacuum state $|0\rangle$, which is defined as

[TABLE]

is an eigenstate of the Hamiltonian operator

[TABLE]

where $T_{00}$ is defined as

[TABLE]

But, $|0\rangle$ is not an eigenstate of the energy density operator

[TABLE]

because of the terms of the form $a_{\mathbf{k}}a_{\mathbf{k}^{\prime}}$ and $a_{\mathbf{k}}^{{\dagger}}a_{\mathbf{k}^{\prime}}^{{\dagger}}$.

Direct calculation shows the magnitude of the fluctuation of the vacuum energy density diverges as the same order as the energy density itself,

[TABLE]

where

[TABLE]

where $\Lambda$ is the effective QFT’s high energy cutoff. (For more details on this calculation, see equation (220) in Appendix A.) Thus the energy density fluctuates as violently as its own magnitude. With such huge fluctuations, the vacuum energy density $\rho^{\mathrm{vac}}$ is not a constant in space or time.

Furthermore, the energy density of the vacuum is not only not a constant in time at a fixed spatial point, it also varies from place to place. In other words, the energy density of vacuum is varying wildly at every spatial point and the variation is not in phase for different spatial points. This results in an extremely inhomogeneous vacuum. The extreme inhomogeneity can be illustrated by directly calculating the expectation value of the square of difference between energy density at different spatial points,

[TABLE]

where $\Delta x=|\mathbf{x}-\mathbf{x}^{\prime}|$ and we have normalized $\Delta\rho^{2}$ by dividing its asymptotic value $\frac{4}{3}\langle T_{00}\rangle^{2}$ (the curly bracket $\{\}$ is the symmetrization operator which is defined by (216)). The behavior of $\Delta\rho^{2}$ for the scalar field (15) in Minkowski vacuum is plotted in FIG. 1, which shows that the magnitude of the energy density difference between two spacial points quickly goes up to the order of $\langle T_{00}\rangle$ itself as their distance increases by only the order of $1/\Lambda$. (For more details on the calculations and how the energy density fluctuates all over the spacetime, see Appendix A.)

As the vacuum is clearly not homogeneous, equation (10) is not valid as it depends on a homogeneous and isotropic matter field and metric. Therefore a new method of relating vacuum energy density to the observed Hubble expansion rate is required.

IV Differences made by the inhomogeneous vacuum—a simple model

The extreme inhomogeneity of the vacuum means its gravitational effect cannot be treated perturbatively, so another method is required. As solutions to the fully general Einstein equations are difficult to obtain, we will first look at a highly simplified model.

IV.1 Beyond the FLRW metric

To describe the gravitational property of the inhomogeneous quantum vacuum, we must allow inhomogeneity in the metric. This is accomplished by allowing the scale factor $a(t)$ in the FLRW metric (9) to have spatial dependence,

[TABLE]

The full Einstein field equations for the coordinate (23) are

[TABLE]

where $\nabla=(\partial_{1},\partial_{2},\partial_{3})$ is the ordinary gradient operator with respect to the spatial coordinates $x,y,z$.

By choosing the above simplest inhomogeneous metric (23), we are assuming a mini-superspace type model, and will choose which of these equations do apply later. This treatment might result in inconsistencies as general vacuum fluctuations of the matter fields posses rich structures that they may not produce spacetime described by the metric (23). To fully describe the resulting inhomogeneous spacetime, one needs a more general metric. However, as a first approximation, using (23) is relatively easy to calculate and leads to interesting results. We are also going to do the calculations for a more general metric in section VIII.

IV.2 The fluctuating spacetime

The role played by the value of vacuum energy density in the above equations (24), (25), (26) and (27) is different from (10). The value of vacuum energy density is no longer directly related to the Hubble rate $H$ through the equation (10). This is evident from the $00$ component of the Einstein equation (24). The equation (10) is only the special case of (24) when the spatial derivatives $\nabla a$ and $\nabla^{2}a$ are zero, which requires that the matter distribution is strictly homogeneous and isotropic. However, as shown in the last section, the quantum vacuum is extremely inhomogeneous and necessarily anisotropic, which requires $\nabla a$ and $\nabla^{2}a$ be huge. This can be seen through the $ij$ component of the Einstein equation (27). In fact, due to symmetry properties of the quantum vacuum, we have the expectation value of shear stress $T_{ij}$ on the right side of (27)

[TABLE]

Meanwhile, $T_{ij}$ must fluctuate since the quantum vacuum is not its eigenstate either, and the magnitude of the fluctuation is on the same order of the vacuum energy density

[TABLE]

This means that the $T_{ij}$ is constantly fluctuating around zero with a huge magnitude of the order of vacuum energy density. As a result, in (27), the spatial derivatives of $a(t,\mathbf{x})$ must also constantly fluctuate with huge magnitudes.

More importantly, since the scale factor $a(t,\mathbf{x})$ is spatially dependent, the physical distance $L$ between two spatial points with comoving coordinates $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$ is no longer related to their comoving distance $\Delta x=|\mathbf{x}_{1}-\mathbf{x}_{2}|$ by the simple equation $L(t)=a(t)\Delta x$ and the observed global Hubble rate $H$ is no longer equal to the local Hubble rate $\dot{a}/a$. Instead, the physical distance and the global Hubble rate are defined as

[TABLE]

and

[TABLE]

where the line element $dl=\sqrt{dx^{2}+dy^{2}+dz^{2}}$.

Equation (31) shows the key difference between the gravitational behavior of quantum vacuum predicted by the homogeneous FLRW metric (9) and by the inhomogeneous metric (23).

For the homogeneous metric (9), the scale factor $a$ is spatially independent and (31) just reduces to

[TABLE]

In this case, there are only two distinct choices for Hubble rates on a spatial slice $t=Const$ under the initial value constraint equation (10)

[TABLE]

which implies that all points in space have to be simultaneously expanding or contracting at the same constant rate (Here we do not include the cosmological constant $\lambda$).

But for the inhomogeneous metric (23), the scale factor $a$ is spatially dependent and there is much more freedom in choosing different local Hubble rates at different spatial points of the slice $t=Const$ under the corresponding initial value constraint equation (24).

In fact, the local Hubble rates must be constantly changing over spatial directions within very small length scales. This can be seen from the initial value constraint equations (26), which can be rewritten as

[TABLE]

where $\mathbf{J}=\left(T_{01},T_{02},T_{03}\right)$ is vacuum energy flux111One might notice that (34) requires $\nabla\times\mathbf{J}=0$, which means that to produce the metric of the form (23), the energy flux of the matter field needs to be curl free. As mentioned in the last paragraph of section IV.1, this is not true for general matter fields, but here as a first approximation we will use (34) to estimate the magnitude of change in $\dot{a}/a$..

The solution to (34) or (26) is

[TABLE]

where $\mathbf{dl}^{\prime}=(dx^{\prime},dy^{\prime},dz^{\prime})$ and $\mathbf{x}_{0}$ is an arbitrary spatial point. The above solution (35) shows that the difference in the local Hubble rates $\dot{a}/a$ between $\mathbf{x}_{0}$ and $\mathbf{x}_{1}$ is determined by the spatial accumulations (integral) of the vacuum energy flux $\mathbf{J}$. Similar to the shear stress, $\mathbf{J}$ has zero expectation value

[TABLE]

but huge fluctuations

[TABLE]

which implies that the local Hubble rates differ from point to point due to the fluctuations. The average of the absolute value of $\dot{a}/a$ can be estimated with the constraint equation (24)

[TABLE]

By using (35), we find that the difference in local Hubble rates becomes comparable with itself for points separated by only a distance of the order $\Delta x\sim\frac{1}{\sqrt{G}\Lambda^{2}}$ as $\Lambda\to+\infty$:

[TABLE]

Up to this point, we have used the equations (24), (26) and (27). These equations are all initial value constraint equations which do not contain the scale factor’s second order time derivative $\ddot{a}$. To get the information about the time evolution of $a(t,\mathbf{x})$, we also need to use (25). A linear combination of equations (24) and (25) gives,

[TABLE]

where all the spatial derivatives of $a$ cancel and only the second order time derivative left. Therefore we reach the following dynamic evolution equation for $a(t,\mathbf{x})$:

[TABLE]

where

[TABLE]

If $\Omega^{2}>0$, which is true if the matter fields satisfy normal energy conditions, (41) describes a harmonic oscillator with time dependent frequency. The most basic behavior of a harmonic oscillator is that it oscillates back and forth around its equilibrium point, which implies that the local Hubble rates $\dot{a}/a$ are periodically changing signs over time. By using equation (39) you can find that $\dot{a}/a$ must also have this periodic sign change in a given spatial direction.

Physically, these fluctuating features of $\dot{a}/a$ imply that, at any instant of time, if the space is expanding in a small region, it has to be contracting in neighboring regions; and at any spatial point, if the space is expanding now, it has to be contracting later.

These features result in huge cancellations when calculating the averaged $H$ through (31). The observable overall net Hubble rate can be small although the absolute value of the local Hubble rate $|\dot{a}/a|$ at each individual point has to be huge to satisfy the constraint equation (24). In other words, while the instantaneous rates of expansion or contraction at a fixed spatial point can be large, their effects can be canceled in a way that the physical distance (30) would not grow $10^{120}$ times larger than what is observed.

This picture of fluctuating spacetime is not completely new. It is similar to the concept of spacetime foam devised by John Wheeler Wheeler:1957mu ; misner1973gravitation that in a quantum theory of gravity spacetime would have a foamy, jittery nature and would consist of many small, ever-changing, regions in which spacetime are not definite, but fluctuates. His reason for this “foamy” picture is the same as ours—at sufficiently small scales the energy of vacuum fluctuations would be large enough to cause significant departures from the smooth spacetime we see at macroscopic scales.

The solution for $a(t,\mathbf{x})$ will be given by equations (51), (58) and (59) in the next section V to describe this foamy structure more precisely.

IV.3 Methods and assumptions in solving the system

In principle, we need a full quantum theory of gravity to solve the evolution details of this quantum gravitational system. Unfortunately, no satisfactory theory of quantum gravity exists yet.

In this paper, we are not trying to quantize gravity. Instead, we are still keeping the spacetime metric $a(t,\mathbf{x})$ as classical, but quantizing the fields propagating on it. The key difference from the usual semiclassical gravity is that we go one more step—instead of assuming the semiclassical Einstein equation, where the curvature of the spacetime is sourced by the expectation value of the quantum field stress energy tensor, we also take the huge fluctuations of the stress energy tensor into account. In our method, the sources of gravity are stochastic classical fields whose stochastic properties are determined by their quantum fluctuations.

The evolution details of the scale factor $a(t,\mathbf{x})$ described by equation(41) depends on the property of the time dependent frequency $\Omega(t,\mathbf{x})$ given by (42). For both simplicity and clarity, in the following sections we investigate the properties of $\Omega$ by considering the contribution from a real massless scalar field $\phi$. In this case, the stress energy tensor for a general spacetime metric $g_{\mu\nu}$ is

[TABLE]

Direct calculation using the inhomogeneous metric (23) gives that

[TABLE]

where all the spatial derivatives and explicit dependence on the metric $a$ are canceled. Thus we obtain

[TABLE]

which is not explicitly dependent on the metric $a(t,\mathbf{x})$.

However, the resulting spacetime sourced by this massless scalar field $\phi$ does have back reaction effect on $\phi$ itself. This is because $\phi$ obeys the equation of motion in curved spacetime

[TABLE]

which reduces to

[TABLE]

for the special metric (23).

Incorporating the back reaction effect by solving both the Einstein equations (24), (25), (26), (27) for the metric $a$ and the equation of motion (47) for the field $\phi$ at the same time is difficult. Fortunately, solving the system in this way is unnecessary. Physically, the quantum vacuum locally behaves as a huge energy reservoir, so that the back reaction effect on it should be small and can be neglected. In our method, we will first assume that the quantized field $\phi$ is still taking the flat spacetime form of (15) for field modes below the effective QFT’s high frequency cutoff $\Lambda$. We use (15) to calculate the stochastic property of the time dependent frequency $\Omega$ and then solve (41) to get the resulting curved spacetime described by the metric $a(t,\mathbf{x})$. This will be done in the next section V.

We then investigate the back reaction effect in section VIIby quantizing the field $\phi$ in the resulting curved spacetime. It turns out that, while the resulting spacetime is fluctuating, the fluctuation happens at scales which are much smaller than the length scale $1/\Lambda$. Therefore the corrections to the field modes with frequencies below the cutoff $\Lambda$ is quite small and thus the flat spacetime quantization (15) is valid to high precision. (See equations (126) (or (158)), (127) and (129) for quantitatively how high this precision is.) In this way we justify neglecting the aforementioned back reaction.

Empirically, this must be true since ordinary QFT has achieved great successes by assuming flat Minkowski background and using the expansion (15). So if our method is correct, (15) has to be still valid even the background spacetime is no longer flat but wildly fluctuating at small scales. In other words, the resulting spacetime should still looks like Minkowskian for low frequency field modes. Long wavelength fields ride over the Wheeler’s foam as if it is not there. This is similar to the behavior of very long wavelength water waves which do not notice the rapidly fluctuating atomic soup over which they ride.

V The solution for $a(t,\mathbf{x})$

In this section we give the solution for the local scale factor $a(t,\mathbf{x})$.

V.1 Parametric resonance

One important feature of a harmonic oscillator with time dependent frequency is that it may exhibit parametric resonance behavior.

If the $\Omega(t,\mathbf{x})$ is strictly periodic in time with a period $T$, the property of the solutions of (41) has been thoroughly studied by Floquet theory teschlordinary . Under certain conditions (for example, the condition (78)), the parametric resonance phenomenan occurs and the general solution of (41) is (see e.g. Eq(27.6) in Chapter V of LANDAU197658 )

[TABLE]

where $H_{\mathbf{x}}>0$, $c_{1}$ and $c_{2}$ are constants. The $P_{1}$ and $P_{2}$ are purely periodic functions of time with period $T$. They are in general functions oscillating around zero. The amplitude of the first term in (48) increases exponentially with time while the second term decreases exponentially. Therefore the first term will become dominant and the solution will approach a pure exponential evolution

[TABLE]

where we have absorbed the constant $c_{1}$ into $P(t,\mathbf{x})$ by letting $P(t,\mathbf{x})=c_{1}P_{1}(t,\mathbf{x})$.

Physically, the exponential evolution of the amplitude of $a(t,\mathbf{x})$ is easy to understand. If $\Omega$ is strictly periodic, the system will finally reach a steady pattern of evolution (when the second term in (48) has been highly suppressed). In this pattern, after each period of evolution of the system, $a$ increases by a fixed ratio, i.e. $a(t+T,\mathbf{x})=\mu_{\mathbf{x}}a(t,\mathbf{x})$, which results in the exponential increase since after $n$ cycles, $a(t+nT,\mathbf{x})=\mu_{\mathbf{x}}^{n}a(t,\mathbf{x})$. Here the $\mu_{\mathbf{x}}$ is related to the $H_{\mathbf{x}}$ by $H_{\mathbf{x}}=\frac{\ln\mu_{\mathbf{x}}}{T}$.

Due to the stochastic nature of quantum fluctuations, the $\Omega(t,\mathbf{x})$ in (41) is not strictly periodic. However, its behavior is still similar to a periodic function. In fact, $\Omega$ exhibits quasiperiodic behavior in the sense that it is always varying around its mean value back and forth on an approximately fixed time scale. To see this, we calculate the following normalized covariance:

[TABLE]

where $\Delta t=t_{1}-t_{2}$ and we have dropped the label $\mathbf{x}$ in the second line of the above definition (V.1) since the final result is independent with $\mathbf{x}$.

Explicit expression for $\chi$ as a function of $\Delta t$ is given by (226), which is plotted in FIG. 2. It describes how $\Omega^{2}$ at different times change around their mean values together. We say that two $\Omega^{2}$ separated by time difference $\Delta t$ are positively (negatively) correlated if $\chi(\Delta t)>0(<0)$, since it means that they are most likely to be at the same (opposite) side of their mean value $\langle\Omega^{2}\rangle$.

FIG. 2 and (226) show that $\Omega^{2}$ at different times are strongly correlated at close range. Especially, the negative correlation is strongest when $\Delta t\sim 2/\Lambda$, which implies that if at $t=0$ the $\Omega^{2}$ is above its mean value $\left\langle\Omega^{2}\right\rangle$, then at $t\sim 2/\Lambda$, it is most likely below $\left\langle\Omega^{2}\right\rangle$. So basically, $\Omega^{2}$ varies around its mean value quasiperiodically on the time scale $T\sim 1/\Lambda$.

This quasiperiodic behavior of $\Omega$ should also lead to parametric resonance behavior seen in (49), instead with a difference in that $H_{\mathbf{x}}$ would become time dependent, i.e. the solution would take the following form

[TABLE]

where $P(t,\mathbf{x})$ here is no longer a strictly periodic function as in (49) but a quasiperiodic function with the same quasiperiod of the order $1/\Lambda$ as the time dependent frequency $\Omega(t,\mathbf{x})$. (The solution (58) for $P(t,\mathbf{x})$ in the next subsection V.2 reveals this property.)

The physical mechanism is similar. The system will also reach a final steady evolution pattern. In this pattern, after each quasiperiod of evolution of the system, $a$ will increase by an approximately fixed ratio. Suppose that during the $i$th cycle of quasiperiod $T_{i}$, $a$ increases by a factor $\mu_{i\mathbf{x}}$, i.e. $a(t+T_{i},\mathbf{x})=\mu_{i\mathbf{x}}a(t,\mathbf{x})$. Then after the $n$ cycles, we have $a(t+\displaystyle\sum_{i=1}^{n}T_{i},\mathbf{x})=\left(\prod_{i=1}^{n}\mu_{i\mathbf{x}}\right)a(t,\mathbf{x})$. Because the quasiperiods $T_{i}$ and the factors $\mu_{i\mathbf{x}}$ are generally different from each other, the exponent in (51) would need to take the form of integration.

The detailed oscillating behavior of $P(t,\mathbf{x})$ is not observable at macroscopic scales. However, the factor of the exponential increase $e^{\int_{0}^{t}H_{\mathbf{x}}(t^{\prime})dt^{\prime}}$ can be observed. In fact, inserting (51) into (30), the observable physical distance would become

[TABLE]

where

[TABLE]

and the global Hubble expansion rate $H$ is

[TABLE]

In the next two subsections, we are going to give the solution for $P(t,\mathbf{x})$ and the global Hubble expansion rate $H$.

V.2 The solution for $P(t,\mathbf{x})$

The magnitude of the time dependent frequency $\Omega$ is of the order $\sim\sqrt{G\left\langle T_{00}\right\rangle}\sim\sqrt{G}\Lambda^{2}$, while $\Omega$ itself varies roughly with a characteristic frequency $\Lambda$ (this has been shown by FIG. 2). Then according to (41), the scale factor $a$ would oscillate with a period that roughly goes as $T=2\pi/\Omega\sim 1/\sqrt{G}\Lambda^{2}\ll 1/\Lambda$, as $\Lambda\to\infty$, where $1/\Lambda$ is the time scale on which the $\Omega$ itself would change significantly.

So comparing to the oscillating period $T$ of the scale factor $a$, the variation of $\Omega$ itself is very slow, although the time $1/\Lambda$ is already very short for large $\Lambda$. Therefore, during one period of the oscillation of $a$, $\Omega$ is almost constant since it has not have a chance to change significantly during such a short time scale. In this sense the time dependent frequency $\Omega$ is slowly varying and the evolution of the scale factor $a$ is an adiabatic process.

The slow variation of $\Omega$ can be verified in a more formal way by calculating the expectation values of $\Omega^{2}=\frac{8\pi G}{3}\dot{\phi}^{2}$ and $\left(\frac{d\Omega}{dt}\right)^{2}=\frac{8\pi G}{3}\ddot{\phi}^{2}$. Using (15), we have

[TABLE]

[TABLE]

(55) just gives $\Omega\sim\sqrt{G}\Lambda^{2}$ as expected, (56) gives $d\Omega/dt\sim\sqrt{G}\Lambda^{3}$. Therefore, during one period of oscillation $T\sim 2\pi/\Omega\sim\frac{1}{\sqrt{G}\Lambda^{2}}$, we have, as $\Lambda\to+\infty$, the slow varying condition (see equation (49.1) in Chapter VII of LANDAU1976131 )

[TABLE]

is satisfied. Thus the system varies adiabatically since $\Omega$ varies only slightly during the one period of oscillation time $T$.

The leading order solution of the equation (41) for a harmonic oscillator with the slowly varying frequency $\Omega$ can be obtained by a first order WKB approximation. This adiabatic approximation neglects the small exponential factor in (51). It gives the solution $P(t,\mathbf{x})$ which is describing the oscillating behavior of $a(t,\mathbf{x})$. The result is,

[TABLE]

The $P(t,\mathbf{x})$ above is a quasiperiodic function with the same quasiperiod of the order $1/\Lambda$ as the time dependent frequency $\Omega(t,\mathbf{x})$ just as expected. The two constants of integration $A_{0}$ and $\theta_{\mathbf{x}}$ in (58) can be determined by the initial values $a(0,\mathbf{x})$ and $\dot{a}(0,\mathbf{x})$.

The quantum vacuum is fluctuating everywhere, but its statistical property must be still the same everywhere. Correspondingly, the statistical property of $P(t,\mathbf{x})$ must also be the same everywhere, which requires that the constant $A_{0}$ to be independent with respect to the spatial coordinate $\mathbf{x}$. In addition, the constant $A_{0}$ can be chosen as any nonzero value since the scale factor $a$ multiplying by any nonzero constant describes physically equivalent spacetimes.

The initial phase $\theta_{\mathbf{x}}$ at different places must be dependent on $\mathbf{x}$. In applying the initial value constraint equation (35), neglecting the small exponential factor in (51) and neglecting the relatively small time derivative terms of the slowly varying frequency $\Omega$, we obtain the result,

[TABLE]

where $\theta_{\mathbf{x}_{0}}$ is the initial phase of the scale factor $a$ at an arbitrary spatial point $\mathbf{x}_{0}$.

In solutions (58) and (59) we see the fluctuating nature of spacetime at very small scales as described in the previous section IV.2. In particular, (59) shows that the phases of $a(t,\mathbf{x})$ vary on a given initial Cauchy slice; some locations contract while others expand. In this new physical picture the catastrophic vacuum energy density is confined to very small scales.

V.3 The global Hubble expansion rate $H$

As the system is adiabatic, the parametric resonance effect is weak. The adiabatic solution (58) in the last subsection does not include the parametric resonance and thus misses the small exponential factor expected in (51). In this subsection we go beyond the adiabatic approximation and investigate the exact strength of the weak parametric resonance.

When considering the weak parametric resonance effect, the constant $A_{0}$ in (58) would become time and space dependent and take the following form

[TABLE]

in order to satisfy (51).

To determine how the $H_{\mathbf{x}}(t)$ depends on the spacetime dependent frequency $\Omega(t,\mathbf{x})$, we consider the adiabatic invariant of a harmonic oscillator with time dependent frequency, which is defined as

[TABLE]

where

[TABLE]

Replace the constant $A_{0}$ in (58) by $A(t,\mathbf{x})$ and then plug it into the above expression (61) we get that

[TABLE]

where we have neglected the time derivatives of $A$ and $\Omega$ in the above equation (63), which are higher order infinitesimals. $I$ is invariant in the first order adiabatic approximation. When going to higher orders, $I$ will slowly change with time. Through the relation (63) between $I$ and $A$ we can obtain how the $A(t,\mathbf{x})$ changes by investigating how accurately the adiabatic invariant is preserved and how it changes with time.

It has been proved by Robnik and Romanovski adibaticinvariant ; 0305-4470-39-1-L05 that, in full generality (no restrictions on the function $\Omega(t,\mathbf{x})$), the final value of the adiabatic invariant for the average energy $\bar{I}=\bar{E}/\Omega$ is always greater or equal to the initial value $I_{0}=E_{0}/\Omega_{0}$ (see the references adibaticinvariant ; 0305-4470-39-1-L05 for precise definition about the average energy). In other words, the average value of the adiabatic invariant $\bar{I}=\bar{E}/\Omega$ for the mean value of the energy never decreases, which is a kind of irreversibility statement. It is conserved only for infinitely slow process, i.e. an ideal adiabatic process.

Therefore, in the case of our quasiperiodic frequency $\Omega(t,\mathbf{x})$ in (41), $\bar{I}$ will also always increase. Moreover, it will increase by a fixed factor after each quasiperiod of evolution, which results in an exponentially increasing $\bar{I}$. This is in fact evident because of the weak parametric resonance effect. In the following we investigate this exponential behavior in detail.

First we construct the evolution equation for the adiabatic invariant $I$. Do the canonical transformation

[TABLE]

Then the evolution equations for $a$ and its conjugate momentum $\dot{a}$ transfer to the evolution equation for the new action variable $I$ and the angle variable $\varphi$,

[TABLE]

Integrating (66) yields

[TABLE]

where

[TABLE]

The $H_{\mathbf{x}}(t^{\prime})$ in the above equation (69) is just the same with the $H_{\mathbf{x}}(t^{\prime})$ defined in (51) and (60), which can be seen by applying equation (63). Thus equation (69) constructed the dependence of $H_{\mathbf{x}}(t^{\prime})$ on the time dependent frequency $\Omega(t^{\prime},\mathbf{x})$.

The observable global Hubble expansion rate $H$ is the average of $H_{\mathbf{x}}(t^{\prime})$ over time, which was defined by equation (54). Plugging (69) into (54) gives,

[TABLE]

When the slow varying condition (57) holds, from equation (67) we know that $d\varphi/dt$ is positive, i.e. $\varphi$ is a monotonic function in time. Thus we can change the integral in (70) from the integration over $t^{\prime}$ to integration over $\varphi^{\prime}$:

[TABLE]

where $\varphi_{0}=\varphi(0)$ and $\varphi=\varphi(t)$.

To evaluate $H$, we formally treat $\varphi$ as a complex variable and close the contour integral in the upper half plane. The integrand in (71) has no singularities for real $\varphi$ if the slow varying condition (57) holds. Equation (67) implies that $\varphi\sim\Omega t\sim\sqrt{G}\Lambda^{2}t$, so the length of the interval $\varphi-\varphi_{0}\sim\sqrt{G}\Lambda^{2}t$ goes to infinity as $\Lambda\to+\infty$. Hence the principle contribution to the integral in (71) comes from the residue values at singularities $\varphi_{(k)}$ inside the contour:

[TABLE]

Each term in (72) gives a contribution containing a factor $\exp\left(-2\operatorname{Im}\varphi_{(k)}\right)$. So the dominant contribution in (72) comes from the singularities near the real axis, i.e. those with the smallest positive imaginary part. To keep the calculation simple, we retain only those terms. Since $\Omega(t)$ varies quasiperiodically with a characteristic time $\tau\sim 1/\Lambda$, the number of singularities near the real axis would roughly be on the order $t/\tau\sim\Lambda t$. Therefore the $H$ in (72) is roughly

[TABLE]

Let $t_{(k)}$ be the (complex) “instant” corresponding to the singularity $\varphi_{(k)}$: $\varphi_{(k)}=\varphi(t_{(k)})\sim\Omega\,t_{(k)}$. In general, $|t_{(k)}|$ has the same order of magnitude as the characteristic time $\tau\sim 1/\Lambda$ of variation of the $\Omega$. Remember that $\Omega\sim\sqrt{G}\Lambda^{2}$, thus the order of magnitude of the exponent in (73) is

[TABLE]

Therefore, inserting (74) into (73) gives

[TABLE]

where $\alpha$ and $\beta$ are two dimensionless constants which depend on the variation details of the time dependent frequency $\Omega(t,\mathbf{x})$. Therefore $H$ becomes exponentially small in the limit of taking $\Lambda$ to infinity. This is a manifestation of the well-established result that the error in adiabatic invariant is exponentially small for analytic $\Omega$ LANDAU1976131 ; adibaticinvariant . In fact, the technique we used in deriving (75) is very similar to the one used in deriving the error in adiabatic invariant in the pages “$160-161$” of LANDAU1976131 .

V.4 A more intuitive explanation

So far we have obtained our key result (75) for the global Hubble expansion rate $H$. To understand the mechanism of weak parametric resonance better, we give a more intuitive explanation in this subsection.

Consider the following simplest parametric oscillator:

[TABLE]

where

[TABLE]

The behavior of the above harmonic oscillator with time dependent frequency has been thoroughly studied (see e.g. eq(27.7) in Chapter V of LANDAU197658 ). The parametric resonance occurs when the frequency $\gamma$ with which $\omega(t)$ varies is close to any value $2\omega_{0}/n$, i.e.

[TABLE]

where $n$ is an integer. The strength of the parametric resonance is strongest if $\gamma$ is nearly twice $\omega_{0}$, i.e. if $n=1$. As $n$ increases to infinity, the strength of the parametric resonance decreases to zero. This is easy to understand since as $n$ increases, the varying frequency $\gamma$ of $\omega(t)$ becomes slower compared to the oscillator’s natural frequency $\omega_{0}$ and as $n\to\infty$, (76) reduces to an ordinary harmonic oscillator with constant frequency which has no parametric resonance behavior.

Now let us go back to Eq.(41) for $a(t,\mathbf{x})$. The time dependent frequency $\Omega(t,\mathbf{x})$ in (41) is more complicated than the $\omega(t)$ given in our example (77). However, it can be written in a similar form:

[TABLE]

where

[TABLE]

and $f(\gamma)$, $g(\gamma)$ are operator coefficients, whose exact form are given by (229) and (230) in Appendix A. The behavior of $\Omega^{2}(t,\mathbf{x})$ for an arbitrary $\mathbf{x}$ is the same with $\Omega^{2}(t,\mathbf{0})$ except phase differences. The power spectrum density of the varying part of $\Omega^{2}(t,\mathbf{0})$ (except for the constant $\Omega_{0}^{2}$ part) given by (232) is plotted in FIG. 3.

Unlike the case (77) where the $\omega(t)$ varies with a single frequency $\gamma$, the $\Omega(t,\mathbf{0})$ in (79) varies with frequencies continuously distributed in the range $(0,2\Lambda)$ with a peak around $1.7\Lambda$ (see FIG. 3). From (80) we have that, as taking the cutoff frequency $\Lambda$ to infinity, $\Omega_{0}\sim\sqrt{G}\Lambda^{2}\gg 2\Lambda$. Because of the continuity of the spectrum of $\Omega$, we can always find integers $n$ such that if

[TABLE]

then

[TABLE]

So $\Omega(t,\mathbf{x})$ always contains frequencies $2\Omega_{0}/n$ that may excite resonances. From (81) we see that $n\to\infty$ as taking the cutoff $\Lambda$ to infinity. While as $n$ increases, the relative magnitude of the resonance frequency $2\Omega_{0}/n$ decreases comparing to the $a(t,\mathbf{x})$’s natural frequency $\Omega_{0}$. Then for reasons similar to the simplest parametric oscillator (76), the strength of the parametric resonance of (41) would also decrease to zero. This weak parametric resonance effect leads to the global Hubble expansion rate

[TABLE]

V.5 Numerical verification

In this subsection, we do a numerical calculation for the evolution equation (41) to verify our result (75), which describes the dependence of $H$ on cutoff $\Lambda$.

In this subsection, Planck units will be used, so all instances of Newton’s constant are set to unity, $G=1$.

The main idea is to rewrite the time dependent frequency $\Omega(t)$ in phase space. (To see more details about this numeric method, please check Appendix B. Here we only list the most crucial results.) For a real massless scalar field, we have

[TABLE]

This is the Weyl transformation of the operator $\hat{\Omega}^{2}(t)$. Here $\{x_{\mathbf{k}},p_{\mathbf{k}}\}$ are phase space points of a particular field mode with momentum $\mathbf{k}$. Approximately, for a particular choice of $\{x_{\mathbf{k}}\},\{p_{\mathbf{k}}\}$, we can get an classic equation for $a$:

[TABLE]

The observed value $a_{o}(t)$ is the average of $a(\{x_{\mathbf{k}}\},\{p_{\mathbf{k}}\},t)$ over the Wigner pseudo distribution function $W(\{x_{\mathbf{k}}\},\{p_{\mathbf{k}}\},t)$, which is based on the wave function of the quantum field:

[TABLE]

If the quantum field is in its ground state, we have

[TABLE]

which means $\{x_{\mathbf{k}}\},\{p_{\mathbf{k}}\}$ are all Gaussian variables. Based on this observation, our method to simulate this equation is as following: i) at first we generate a set of random Gaussian numbers for $\{x_{\mathbf{k}}\},\{p_{\mathbf{k}}\}$ ; ii) we solve the equation (85) for this particular set of numbers; iii) then we repeat the process for another set of random numbers until a certain amount of repetitions; iv) The result $a_{o}(t)$ is the average over all samples we have generated. We choose the repetition amount to be big enough for the results to converge. The result of a single scalar field case is illustrated in Fig. 4.

We can find that the slope of $\log|a_{o}(t)|\sim t$ is decreasing as we increase the cutoff $\Lambda$ as we expect. For the single field, $\Omega^{2}(t)=\frac{8\pi}{3}\dot{\phi}^{2}$ repeatedly reaches zero since classically $\dot{\phi}$ is continuous and oscillates from positive to negative. Around these zero points the slow varying condition for $\Omega^{2}$ is not satisfied. But because the time duration of reaching zero is very short, this would not cause the adiabatic expanding scheme to breaks down. This point reveals in the numerical calculation.

The real Universe contains many different quantum fields. In Fig. 5, we show the result when we include two independent massless scalar fields in which $\Omega^{2}(t)=\frac{8\pi}{3}(\dot{\phi}_{1}^{2}+\dot{\phi}_{2}^{2})$. In this case $\Omega^{2}$ would not reach zero at almost all times since this can happen only when both $\dot{\phi}_{1}$ and $\dot{\phi}_{2}$ pass zero, which is unlikely to happen frequently.

In the two field case, we plot the $\log(H/\Lambda)\sim\Lambda$ graph to verify the quantitative relation (75). The result is illustrated in Fig.6. We can see that for $\Lambda\geq 10$, the result shows decent linearity, which is what we expected since the derivation of (75) is only valid for large $\Lambda$. In this case, the two constants $\alpha=e^{4.6}\approx 100$ and $\beta=0.12$ according to the numeric calculation.

VI Meaning of our results

It is interesting to notice that both (12) and (52) give the exponential evolution and predict an accelerated expanding Universe. However, the underlying mechanisms are completely different, which leads to opposite results on the predicted magnitude of the observable Hubble expansion rate $H$.

The solution (12) is based on the assumption that quantum vacuum energy density is constant all over the spacetime, which is a necessary requirement if one suppose that vacuum acts as a cosmological constant. This assumption leads to a huge Hubble expansion rate

[TABLE]

as taking the high energy cutoff $\Lambda$ to infinity.

Our proposal (52) is based on the fact that quantum vacuum energy density is constantly fluctuating and extremely inhomogeneous all over the whole spacetime. This fact leads to a small Hubble expansion rate given by (75) which goes to zero as taking the high energy cutoff $\Lambda$ to infinity.

If we can literally take the cutoff $\Lambda$ in (75) to infinity, then $H=0$. In this sense, at least the “old” cosmological constant problem would be resolved.

In principle, this effective theory is valid only up to a large but finite cutoff $\Lambda$, which leads to a tiny but nonzero $H$. Since $H\to 0$ as $\Lambda\to+\infty$, there always exists a very large cutoff value of $\Lambda$ such that $H=\sqrt{\Omega_{\lambda}}H_{0}\approx 1.2\times 10^{-42}\,\mbox{GeV}$ to match the observation, where $H_{0}$ is current observed Hubble constant.

So our result suggests that there is no necessity to introduce the cosmological constant, which is required to be fine tuned to an accuracy of $10^{-120}$, or other forms of dark energy, which are required to have peculiar negative pressure, to explain the observed accelerating expansion of the Universe.

The exact value of $\Lambda$ cannot be determined since we do not know the values of the two dimensionless parameters $\alpha$ and $\beta$ in (75). In principle, we need the knowledge of all fundamental fields in the Universe to determine $\alpha$ and $\beta$, this deserves further investigations in the future and might provide some hint on elementary particle physics.

The value of $\Lambda$ should be on the order of Planck energy or higher. According to the numerical calculation in the last subsection, $\Lambda\sim 1000E_{P}$ if we consider contributions to $\Omega^{2}$ from only two scalar fields. If more fundamental fields are included, we expect the value of $\beta$ would increase and thus decrease the value of $\Lambda$ needed. This is because that it increases the mean value of $\Omega$ and, as a consequence, reduces the ratio between the variation of $\Omega$ over its mean value $\langle\Omega\rangle$ that $\Omega$ varies slower. A slower $\Omega$ leads to smaller $H$ since the parametric resonance is weaker.

VII The back reaction

In this section, we investigate the back reaction effect by quantizing the field $\phi$ in the resulting curved spacetime to justify our method of using the quantized field expansion (15) in Minkowski spacetime as an approximation.

The standard way to quantize the scalar field $\phi$ in a generic curved spacetime $g_{\mu\nu}$ is by first defining the following inner product on a spacelike hypersurface $\Sigma$ with induced metric $h_{ij}$ and unit normal vector $n^{\mu}$ (see e.g. carroll2004spacetime ; 9780511622632 ):

[TABLE]

where $h=\det h_{ij}$ and $\phi_{1}$, $\phi_{2}$ are solutions to the equation (46). The above inner product is independent of the choice of $\Sigma$.

One then choose a complete set of mode solutions $u_{\mathbf{k}}$ of (46) which are orthonormal in the product (89):

[TABLE]

Then the field $\phi$ may be expanded as

[TABLE]

For the flat Minkowski spacetime, i.e. $g_{\mu\nu}=\eta_{\mu\nu}$, (46) reduces to the usual wave equation

[TABLE]

In this case, the mode solutions are usually chosen as

[TABLE]

where $\omega=|\mathbf{k}|$. Plugging (95) into (93) just gives the usual quantum field expansion (15).

For our specific metric (23), (46) reduces to (47). In this case, since the rate of accelerating expansion is extremely small, the back reaction effect due to the macroscopic expansion of the Universe is only important on large cosmological time scales. For this reason, we only worry about the back reaction due to the wildly fluctuating spacetime at small scales. i.e. we neglect the small exponential factor in (51) and use the form of the $a$ based on the solution (58):

[TABLE]

where

[TABLE]

Then (47) becomes

[TABLE]

In order to understand the effect from back reaction, we need to find out how the mode solutions of the above equation (98) in the resulting curved spacetime change from the mode solutions (95) of the equation (94) in the flat Minkowski spacetime.

Physically, the correction to (95) should be small for wave modes with frequencies lower than the cutoff frequency $\Lambda$. That is because the wave length of those field modes is larger than $2\pi/\Lambda$, while our spacetime fluctuates on the length scale $2\pi/\Omega\sim 1/(\sqrt{G}\Lambda^{2})\ll 2\pi/\Lambda$. The relatively long wave length modes should not be sensitive to what is happening on small scales. This is analogous to the situation of sound waves traveling in the medium such as air or water or solids. The medium is constantly fluctuating at atomic scales, but this fluctuation does not affect the propagation of the sound wave whose wavelength is much larger than the atomic scale. Similarly, the propagation of the field modes in the “medium”–the spacetime, which is constantly fluctuating on scales much smaller than the wavelength of the field modes, should also not be affected.

Mathematical demonstration will be given in the following subsections.

VII.1 A simplified toy model

It is complicated to obtain the mode solutions of (98) for a generic stochastic function $\Theta(t,\mathbf{x})$ whose stochastic property is determined by the quantum nature of the field $\phi$. To illustrate the underlying physical mechanism more clearly, we start with a simplified toy model by restricting the phase angle $\Theta(t,\mathbf{x})$ defined by (97) to take the following form:

[TABLE]

where both $\Omega$ and $\mathbf{K}$ are constants and they have the same order of magnitude $\Omega\sim|\mathbf{K}|\sim\sqrt{G}\Lambda^{2}$.

Of course this toy model does not describe the real spacetime sourced by the quantum vacuum since the $\Omega$ is by no means a constant but always varying, although the varying is slow compared to it own magnitude. However, this toy model possesses the key property needed — the spacetime is constantly fluctuating. It will be convenient for visualizing the back reaction effect from a fluctuating spacetime.

After setting the $\Omega\equiv Constant$ and the phase angle $\Theta(t,\mathbf{x})$ to be the form of (99), the equation of motion (98) for $\phi$ becomes

[TABLE]

where we have set $A_{0}=\sqrt{2\Omega}$ such that the average of the coefficient $\frac{A_{0}^{2}}{\Omega}\cos^{2}\Theta$ before $\ddot{\phi}$ is $1$ for convenience.

In the flat spacetime case (94), each mode solution $u_{\mathbf{k}}$ in (95) contains only one single frequency. However, for the above fluctuating spacetime case (100), high frequencies mixes with low frequencies and each mode solution must contain multiple frequencies. In fact, since (100) describes a strictly periodic system with time period $\pi/\Omega$ and spatial period $\pi/|\mathbf{K}|$, each mode solution $u_{\mathbf{k}}$ must change from (95) to the following form:

[TABLE]

where $c_{m}$ are constants.

Inserting (101) into (100) and using the orthogonality of $e^{2im(\Omega t+\mathbf{K}\cdot\mathbf{x})}$, we obtain the following infinite system of linear equations:

[TABLE]

In the above calculations, we have used the Fourier series expansion

[TABLE]

to expand the term $\tan(\Omega t+\mathbf{K}\cdot\mathbf{x})$ in (100).

For the equations of $m\leq-1$, we successively add the ($m+1$)th equation to the $m$th equation by the order from $m=-\infty$ to $m=-1$; and for the equations of $m\geq 1$, we successively add the ($m-1$)th equation to the $m$th equation by the order from $m=+\infty$ to $m=1$. Most terms can be eliminated by these elementary row operations and the above infinite system of linear equations (VII.1) becomes

[TABLE]

To characterize the property of the solutions of this system more clearly, we define the following parameters for convenience:

[TABLE]

As mentioned before that our effective theory has a cutoff $\Lambda$ such that only modes with $\omega,|\mathbf{k}|\leq\Lambda$ are relevant, which are much smaller than $\Omega\sim|\mathbf{K}|\sim\sqrt{G}\Lambda^{2}$ as $\Lambda$ grows large. Therefore, we are only interested in the solutions of (VII.1) or (104) when $\omega,|\mathbf{k}|\ll\Omega$, i.e. when $\epsilon,\upsilon\to 0$.

Dividing both sides of (104) by $\Omega^{2}$ and doing some necessary algebraic manipulations, (104) can be rewritten as

[TABLE]

As $\epsilon,\upsilon\to 0$, the leading order asymptotic solution for $\{c_{n}\}$ of the above system of linear equations (106) depends only on the leading order of the coefficients before $\{c_{n}\}$. By keeping only the leading term for each coefficient, we obtain that the leading order solution of (106) for $\{c_{n}\}$ satisfies the following infinite system of linear equations:

[TABLE]

We will denote the infinite matrix in the above equation (107) by $B$ and its elements by $b_{mn}$ with $-\infty<m,n<+\infty$. In order to have a nonzero solution, the determinant of $B$ must be zero. This gives us the dispersion relation that $\epsilon$ and $\upsilon$ must satisfy in the asymptotic regime $\epsilon,\upsilon\to 0$.

The determinant can be calculated by Laplace expansion:

[TABLE]

where $M_{0n}$ is the $0,n$ minor of $B$, i.e. the infinite determinant that results from deleting the [math]th row and the $n$th column of B. Due to the symmetry property of $B$, we have that, for each $n\neq 0$,

[TABLE]

which implies that all the terms inside the summation symbol $\sum$ of (108) exactly cancel. Therefore, only the first term in (108) survive and thus we have that

[TABLE]

which leads to

[TABLE]

or equivalently

[TABLE]

This proves that the usual dispersion relation still holds for low frequency field modes.

After setting $\epsilon^{2}=\upsilon^{2}$, we start solving the infinite system (107).

First, we rewrite (107) as the following form:

[TABLE]

Notice that the matrix elements of $B$ has the following symmetry properties:

[TABLE]

The above symmetry properties leads to the following relation

[TABLE]

which implies that we only need to solve $c_{n}$ for $n>0$ to solve the whole system.

For convenience, we define the following new variables $x_{n}$ by

[TABLE]

Then using the relation (116), the infinite system of linear equations (113) simplifies to the following infinite recurrence equations:

[TABLE]

where the dependence on $\epsilon$ in the equation (107) or (113) has been eliminated by introducing the new variables $x_{n},n\neq 0$ through (117) and the solution for $x_{n}$ depends only on $\delta$.

In order to find the general formula for the sequence $\{x_{m}\}$, we define the following new variables:

[TABLE]

Then the recurrence equations (120) become

[TABLE]

Sequences satisfying (122) must take the following form:

[TABLE]

where $D$ and $\psi$ are two constants and $\vartheta$ is determined by

[TABLE]

Combining (123) and (121), the general formula for $x_{m}$ can be obtained by iteration

[TABLE]

Replacing the $c_{m}$ in (101) by $x_{m}$ through (117) we obtain that, as $\epsilon\to 0$, the mode solution $u_{\mathbf{k}}(t,\mathbf{x})$ is asymptotic to

[TABLE]

where $x_{m}$ is determined by (116), (117), (118), (119), (120) and (VII.1).

Using the orthogonality of $e^{i2m\left(\Omega t+\mathbf{K}\cdot\mathbf{x}\right)}$, the relative magnitude of the correction to $u_{\mathbf{k}}$ from the usual plane wave mode $e^{-i(\omega t-\mathbf{k}\cdot\mathbf{x})}$ in Minkowski spacetime can be characterized by applying Parseval’s identity:

[TABLE]

From the solution (VII.1) we know that as $m\to\infty$,

[TABLE]

Thus the summation inside the bracket of (127) converges and the correction

[TABLE]

Thus we have demonstrated that the low frequency wave modes ($\omega\leq\Lambda$) are almost not affected by the fluctuating spacetime with much higher frequency ($\Omega\sim\sqrt{G}\Lambda^{2}$).

VII.2 General case

The methods used and results obtained in the last subsection for the particular simplified toy model (100) can be generalized to the generic case (98). To start, we rewrite (98) to the following form:

[TABLE]

where

[TABLE]

For convenience, we choose the constant $A_{0}$ such that the average of $f_{1}$

[TABLE]

Unlike the toy model (100) we used in the last subsection, (130) is not strictly periodic. However, (130) is quasiperiodic and its quasiperiod is the same as the period of (100). This property is reflected in the Fourier transforms $f_{1}(\omega,\mathbf{k})$, $f_{2}(\omega,\mathbf{k})$ and $\mathbf{f}_{3}(\omega,\mathbf{k})$ of the functions $f_{1}(t,\mathbf{x})$, $f_{2}(t,\mathbf{x})$ and $\mathbf{f}_{3}(t,\mathbf{x})$ respectively which are defined by

[TABLE]

For the function $f_{1}(t,\mathbf{x})$ defined by (131), after setting the constant $A_{0}$ by (135) and considering the slow varying property of $\Omega(t,\mathbf{x})$ and $\Theta(t,\mathbf{x})$ in both temporal and spatial directions, its leading order goes as

[TABLE]

which implies that the Fourier transform $f_{1}(\omega,\mathbf{k})$ would have two peaks centered at

[TABLE]

For the function $f_{2}(t,\mathbf{x})$ defined by (132), the second term which includes the factor $\sin 2\Theta$ is dominant since the first term which includes the factor $\dot{\Omega}/\Omega^{2}$ goes as $\sim 1/\Lambda\to 0$ due to the slow varying condition described by (55) and (56). Thus, its leading order goes as

[TABLE]

which implies that the Fourier transform $f_{2}(\omega,\mathbf{k})$ would also have two peaks centered at

[TABLE]

Similarly, for the function $\mathbf{f}_{3}(t,\mathbf{x})$ defined by (133), the second term which includes the factor $\tan\Theta$ is dominant since the absolute value of the first term which includes the factor $\nabla\Omega/(\Omega K_{0})$ also goes as $\sim 1/\Lambda\to 0$ due to the slow varying property of $\Omega$ in spatial directions. Thus, its leading order goes as

[TABLE]

Then using the Fourier series expansion (103) for $\tan\Theta$, we know that the Fourier transform $\mathbf{f}_{3}(\omega,\mathbf{k})$ would have infinitely many peaks centered at

[TABLE]

(For a rough calculation of the above Fourier transforms, see Appendix C)

In addition, we have the zero frequency component (see (269) in Appendix C)

[TABLE]

In summary, the system described by (130) is very similar to the system described by the simplified toy model (100). The only difference is that the Fourier transforms of the coefficients $f_{1}$, $f_{2}$, and $\mathbf{f}_{3}$ in (130) spread around center points given by (140), (142) and (144) while the Fourier transforms of the corresponding coefficients in (100) are ideal delta functions exactly located at same points given by (140), (142) and (144).

Therefore, the mode solution of (130) would take the form similar to (101):

[TABLE]

where $u_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ is non-negligible only when $\omega^{\prime},\mathbf{k}^{\prime}$ are taking values around the centers given by (140), (142) and (144).

Inserting (146) into (130) and replacing the coefficients $f_{1}(t,\mathbf{x})$, $f_{2}(t,\mathbf{x})$ and $\mathbf{f}_{3}(t,\mathbf{x})$ in (130) by the equations (136), (137) and (138) and then using the orthogonality of $e^{i(\omega^{\prime}t+\mathbf{k}^{\prime}\cdot\mathbf{x})}$, we obtain the following uncountably infinite system of linear equations which are similar to (VII.1):

[TABLE]

where we have defined the notation $u_{\mathbf{k}}(0,\mathbf{0})=c_{0}\delta(0,\mathbf{0})$ for convenience.

To characterize the property of the solutions of this system more clearly, we define the following parameters similar to (105) for convenience:

[TABLE]

Dividing both sides of (VII.2) by $\Omega_{0}^{2}$ gives

[TABLE]

Similar to the toy model case, as $\epsilon,\upsilon\to 0$, the leading order solution of (VII.2) for $u_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ satisfies the following uncountably infinite system of linear equations:

[TABLE]

where we have used the property (145) in obtaining (150) from (VII.2).

The above uncountably infinite system of linear equations (150) can also be written formally in matrix form similar to (107). We use similar notations that denoting the matrix here by $B$ and its elements by $b_{(\omega^{\prime},\mathbf{k}^{\prime}),(\omega^{\prime\prime},\mathbf{k}^{\prime\prime})}$ for convenience.

In order to have nonzero solutions, the determinant of the uncountably infinite matrix $B$ has to be zero, which gives the dispersion relations that $\epsilon$ and $\upsilon$ must be satisfied in the asymptotic region $\epsilon,\upsilon\to 0$.

The “determinant” of $B$ can be formally calculated through Laplace expansion similar to (108):

[TABLE]

where $M_{(0,\mathbf{0}),(\omega^{\prime\prime},\mathbf{k}^{\prime\prime})}$ is the $(0,\mathbf{0}),(\omega^{\prime\prime},\mathbf{k}^{\prime\prime})$ minor of $B$, i.e. the ‘determinant’ resulting from deleting the $(0,\mathbf{0})$th row and $(\omega^{\prime\prime},\mathbf{k}^{\prime\prime})$th column of $B$.

Notice that since $f_{1}(t,\mathbf{x})$, $f_{2}(t,\mathbf{x})$ and $\mathbf{f}_{3}(t,\mathbf{x})$ are all real, their Fourier transforms $f_{1}(\omega,\mathbf{k})$, $f_{2}(\omega,\mathbf{k})$ and $\mathbf{f}_{3}(\omega,\mathbf{k})$ defined by (136), (137) and (138) must satisfy the following relations:

[TABLE]

where the $*$ means complex conjugate.

The above symmetry property (152) leads to

[TABLE]

which implies that all the terms inside the integral symbol $\int$ of (151) exactly cancel. Therefore, only the first term in (151) survives and thus we have

[TABLE]

which gives again the usual dispersion relation

[TABLE]

After setting the dispersion relation (155), we only need to solve the $(\omega^{\prime},\mathbf{k}^{\prime})\neq(0,\mathbf{0})$th equations in (150) since $\operatorname{det}B=0$ implies that the $(\omega^{\prime},\mathbf{k}^{\prime})=(0,\mathbf{0})$th equation is redundant.

For convenience, we define new variables $x_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ similar to the $x_{n}$ defined in (117):

[TABLE]

Then (150) can be rewritten as

[TABLE]

Replacing the $u_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ in (146) by $x_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ through (156) we obtain that, as $\epsilon\to 0$, the mode solution $u_{\mathbf{k}}(t,\mathbf{x})$ is asymptotic to

[TABLE]

where $x_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ is determined by (157).

Analogous to (VII.1) in the simplified toy model, $x_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ would also go as

[TABLE]

when $\omega^{\prime},\mathbf{k}^{\prime}$ taking values around the centers

[TABLE]

($x_{\mathbf{k}}(\omega^{\prime},\mathbf{k}^{\prime})$ is negligible if $\omega^{\prime},\mathbf{k}^{\prime}$ is far away from these centers).

Due to Parseval’s theorem, (159) implies that the integral inside the bracket of (158) converges which is similar to (127) and thus the correction to $u_{\mathbf{k}}(t,\mathbf{x})$ also goes as $\epsilon$.

Therefore, when we quantize the scalar field $\phi$ in our wildly fluctuating spacetime by expanding it in terms of the annihilation and creation operators according to (93), the leading order would still be the form of the Minkowski quantum field expansion (15). The correction to the dispersion relation $\omega^{2}=\mathbf{k}^{2}$ and the plane wave mode $e^{-i(\omega t-\mathbf{k}\cdot\mathbf{x})}$ are on the order $\sim\epsilon$. In addition, the extra wave modes which mixing in (126) or (158) are all modes with frequencies higher than $\Omega_{0}\sim\sqrt{G}\Lambda^{2}$, which is much larger than our effective QFT’s cutoff $\Lambda$. These extremely high frequency modes beyond the cutoff are irrelevant to our low energy physics. This also explains why the ordinary QFT works by assuming fixed Minkowski spacetime. The small scale structure averages out in its effect on the long wavelength low energy fields.

In summary, we have argued that although our spacetime sourced by the quantum vacuum is highly curved and wildly fluctuating, the back reaction of the resulting spacetime on the quantum field sitting on it is small. This justifies our method of neglecting back reaction and using the quantum field expansion (15) in Minkowski spacetime at the beginning.

VIII the more general inhomogeneous metrics

In previous sections we assume the simplest inhomogeneous metric (23) to describe the spacetime resulting from the inhomogeneous vacuum. In this section, we try to generalize the result to more general inhomogeneous metrics.

The quantum fluctuations of the vacuum is not completely arbitrary, the magnitude of the fluctuations must be the same everywhere and in every direction, i.e. the spacetime is still stochastically homogeneous and isotropic. Thus we can always choose a special gauge and construct the following general synchronous coordinate:

[TABLE]

For the above metric (161), we employ the initial value formulation of general relativity. In this formulation, the Einstein equation is equivalent to six equations for the evolution of the second fundamental form

[TABLE]

plus the usual four constraint equations,

[TABLE]

[TABLE]

where $k_{ab}=\frac{1}{2}\dot{h}_{ab}$, $k^{ab}=h^{ac}h^{bd}k_{cd}$, $trk=h^{ab}k_{ab}$, $\rho=T_{00}$, $j_{b}=h_{b}^{a}T_{0a}$, $trT=h^{ab}T_{ab}$, $R^{(3)}$ is the 3-dimensional spatial curvature and $D_{a}$ is the derivative operator associated with $h_{ab}$.

Taking trace on both sides of (162) and then combining with (163) gives:

[TABLE]

It is interesting to notice that there are no spatial derivatives included on the left hand of the above equation (165). The key evolution equation (41) for $a(t,\mathbf{x})$ we used in previous sections is just the special case of the above equation (165).

Direct calculation using the expression (43) shows that, the contribution from a real massless scalar field to the right-hand side of (165) is

[TABLE]

where all the spatial derivatives of $\phi$ and all the explicit dependence on the metric $g_{\mu\nu}$ in the definition of stress energy tensor (43) are canceled. It is also interesting to notice that the above exact expression (166) is exactly the same with the corresponding expression (44) for the simplest inhomogeneous metric (23) case.

We first consider the following special case:

[TABLE]

The spacetime described by the above coordinate (167) possesses more freedoms than (23) and thus would exhibit richer structures. In this case, the expansion rate at the same point becomes directionally dependent. Along the three principle axes $\hat{x}$, $\hat{y}$ and $\hat{z}$, which are eigenvectors of the symmetric matrix $h_{ab}$ in (167), the expansion rates $\dot{a}/a$, $\dot{b}/b$ and $\dot{c}/c$ can be different. This means that, at one same point, the space can be expanding in one or two directions and contracting on the other two or one directions.

Under the coordinate system (167), equation (165) becomes

[TABLE]

which is a generalization of the key evolution equation (41) we used in previous sections.

Let

[TABLE]

then (168) immediately leads to

[TABLE]

As the functions $a$, $b$ and $c$ are alternately symmetric, their expectation values must be equal

[TABLE]

Unlike equation (41), here $\Omega_{i}^{2}(t,\mathbf{x})$ does not necessarily go exactly the same as $\frac{4\pi G}{3}\left(\rho+trT\right)$. Since $4\pi G(\rho+trT)$ is slowly varying, $\Omega_{i}^{2}$ must also be slowly varying functions. Otherwise we would have two or three fast varying functions sum together and precisely cancel each other to give a slowly varying function, which is almost impossible in the system with such huge quantum fluctuations. Thus the evolution of $a$, $b$ and $c$ are still adiabatic processes and the conclusion we obtained in previous sections still holds. We can obtain solutions similar to (51) that

[TABLE]

and on average

[TABLE]

where

[TABLE]

Therefore, the observable physical volume is,

[TABLE]

where $h=\det h_{ab}=a^{2}b^{2}c^{2}$.

Next we investigate the most general case

[TABLE]

In this case, the three orthogonal eigenvectors of the symmetric matrix $h_{ab}$ can rotate in space. This gives more freedom and structure to the spacetime evolution than in the case described by the coordinate system (167). For example, an initial sphere will distort toward an ellipsoid with principle axes given by eigenvectors of $h_{ab}$, with rates given by time derivatives $\dot{\lambda}_{i}/\lambda_{i}$ of the corresponding eigenvalues $\lambda_{i}^{2}(t,\mathbf{x}),i=1,2,3$.

Expanding the dynamic equation (165) using the metric (178) gives

[TABLE]

where $h=\mathrm{det}(h_{ab})$ is the determinant of the matrix (178), $h^{\ast}_{ab}$ is the matrix’s $(a,b)$ cofactor and $F$ is a nonlinear function of the metric components $h_{ab}$ and their first time derivatives $\dot{h}_{ab}$. (179) is difficult to handle. Further investigations are needed in the future.

However, the results we obtained for the coordinate (167) suggest that, for the most general case (178), the eigenvalues $\lambda_{i}^{2}(t,\mathbf{x})$ should also evolve adiabatically similar to $a^{2}$, $b^{2}$ and $c^{2}$. In other words, we expect that the results (172), (173), (174) and (175) can be generalized to $\lambda_{i}$ in the most general case and the physical volume of space would expand as

[TABLE]

where $H$ is determined by (75).

IX Discussion

So far, we have presented a new mechanism of vacuum gravitation and showed that it leads to a slow accelerating expansion instead of a catastrophic explosion of the Universe. In this section, we discuss some questions raised and a couple of new concepts suggested by this different way of vacuum gravitating.

IX.1 Lorentz invariant cutoffs

A potential concern is that the cutoff $\Lambda$ we are using is not Lorentz invariant. However, the results would not change if using Lorentz invariant cutoffs instead. This is because the $\Lambda$ is just used for comparing the magnitude of different infinities, whose leading order dependencies on $\Lambda$ can also be obtained just by dimensional analysis. In previous sections, we have taken this simple non-Lorentz invariant cutoff $\Lambda$ for convenience. In this subsection, we use more complicated but Lorentz invariant Pauli-Villars type cutoffs to show directly that the results do not change.

First, we calculate the two point function

[TABLE]

where $\Delta t=t_{1}-t_{2}$. We then replace the photon propagator

[TABLE]

in (181) by

[TABLE]

where $\Lambda_{1}$ can be thought of as a fictitious heavy photon, which can serve as a Lorentz invariant cutoff.

Then (181) becomes

[TABLE]

where we have used the identity

[TABLE]

(184) is still logarithmically divergent when setting $\Delta t=0$. To make it converge, we employ another Pauli-Villars type cutoff $\Lambda_{2}$ by replacing the

[TABLE]

in (184) with

[TABLE]

Then (184) becomes

[TABLE]

where we have used the identity (185) again in obtaining (188).

Performing the integration $\int d^{3}k$ first, (188) becomes

[TABLE]

where

[TABLE]

Then let $\omega=\tilde{\Lambda}(\alpha_{1},\alpha_{2})u$, (189) becomes

[TABLE]

where the integrand is,

[TABLE]

If $\alpha_{2}\neq\alpha_{2}^{\prime}$, the integrand is multivalued in complex plane and has four branch points

[TABLE]

.

Assuming $\Delta t>0$, the integral contour goes around the lower half plane. Without loss of generality assuming $\alpha_{2}>\alpha_{2}^{\prime}$, we can choose the branch cut being the line connecting $u_{1}=\sqrt{1-\frac{i\epsilon}{\tilde{\Lambda}^{2}}}$ and $u_{3}=\frac{\tilde{\Lambda}^{\prime}}{\tilde{\Lambda}}\sqrt{1-\frac{i\epsilon}{\tilde{\Lambda}^{\prime 2}}}$. The integral $I_{0}$ is determined by the integral along this branch cut. Since $I_{0}\rightarrow\frac{1}{\sqrt{u-u_{1,3}}}$ when $u\rightarrow u_{1,3}$, we can easily see that the integral along the branch cut converges to a finite number.

If $\alpha_{2}=\alpha_{2}^{\prime}$, $I_{0}$ becomes $\left(u^{2}-1+\frac{i\epsilon}{\tilde{\Lambda}^{2}(\alpha_{1},\alpha_{2})}\right)^{-5/2}/8$. The integral can be solved easily in this case, i.e. $\int_{0}^{\infty}duI_{0}=\frac{1}{12}$ when $\epsilon\rightarrow 0$ , which is also finite.

Therefore, after setting $\Delta t=0$, we have

[TABLE]

Similarly, we have

[TABLE]

which gives the same result as (55).

When calculating the quantity $\left\langle\ddot{\phi}^{2}\right\rangle$ from (191) using the same technique as in (IX.1), we would get the integral $\int du\,u^{4}I_{0}$, which is logarithmically divergent. To make it converge, we can employ another two Pauli-Villars type cutoffs $\Lambda_{3}$ and $\Lambda_{4}$ just as what we did from (186) to (187). The exact dependence of $\left\langle\ddot{\phi}^{2}\right\rangle$ on $\Lambda_{1}$, $\Lambda_{2}$, $\Lambda_{3}$ and $\Lambda_{4}$ is complicated, but the result for its leading order goes as

[TABLE]

which gives the same result as (56).

So far, we have obtained that the leading order dependencies of $\left\langle\dot{\phi}^{2}\right\rangle$ and $\left\langle\ddot{\phi}^{2}\right\rangle$ on the non-Lorentz invariant sharp cutoff $\Lambda$ and on the Lorentz invariant Pauli-Villars type cutoffs $\Lambda_{1}$, $\Lambda_{2}$, ($\Lambda_{3}$, $\Lambda_{4}$) are the same. As we mentioned in the beginning of this subsection, these results are natural since they could have been guessed by dimensional analysis. Therefore, the slow varying condition (57) is still satisfied by using Lorentz invariant cutoffs that our results would not change.

IX.2 The singularities at $a(t,\mathbf{x})=0$

In our way of vacuum gravitating, the space is alternatively expanding and contracting at each spatial point, and, during each such cycle, the expansion outweighs the contraction a little bit due to the weak parametric resonance effect. This process gives a slowly increasing amplitude $A(t,\mathbf{x})$ of the scale factor $a(t,\mathbf{x})$, whose observable effect is just the accelerating expansion of our Universe.

Probably one of the biggest concerns about this physical picture is the appearance of singularities at points $a(t,\mathbf{x})=0$—according to the solution (58), the scale factor $a(t,\mathbf{x})$ must go through zero whenever the space at $\mathbf{x}$ switches from contraction phase to expansion phase.

Singularities are a generic feature of the solution of Einstein field equations under rather general energy conditions (e.g. strong, weak, dominant etc.), which is guaranteed by Penrose-Hawking singularity theorems PhysRevLett.14.57 ; Hawking511 ; Hawking490 ; Hawking187 ; Hawking:1973uf ; Hawking:1969sw . In this paper, since we investigate the gravitational property of quantum vacuum without modifying either QFT or GR, the appearance of singularities is inevitable—QFT predicts a huge vacuum energy, and according to GR, huge energy must collapse to form singularity.

It is usually thought that the Einstein field equations break down at singularities and thus the spacetime evolution will stop once the singularity is formed. However, it is not the case for our solution to the key dynamic evolution equation (41). (41) describes the oscillating motion a harmonic oscillator. It is natural for a harmonic oscillator to pass its equilibrium point $a(t,\mathbf{x})=0$ at maximum speed without stopping. So in our solution, the singularity immediately disappears after it forms and the spacetime continues to evolve without stopping. Singularities just serve as the turning points at which the space switches from contraction phase to expansion phase.

IX.2.1 Resolving singularity by multiplying $a$

In order to understand why in our solution the singularity is not the end of spacetime evolution, it is helpful to review one crucial step in deriving (41) from (40). Rigorously speaking, we can only obtain the following equation from (40):

[TABLE]

which is not equivalent to (41). To get (41), we need one more step—multiply both sides of (197) by $a$.

Mathematically, $a$ is not allowed to be zero in (197) since it is in the denominator. In fact, when writing down the Einstein field equations (24), (25), (26) and (27), it has been presumed that $a\neq 0$ since if $a=0$, the metric would become degenerate ($g=\operatorname{det}(g_{\mu\nu})=-a^{6}=0$), the curvature would become infinite and the Einstein tensor are simply not defined there.

But, after the inequivalent algebraic manipulation of multiplying both sides of (197) by $a$, $a$ is allowed to evolve to zero in the resulting equation (41) since there is nothing wrong for a harmonic oscillator to go through its equilibrium point. In this sense, we have smoothly extended the solution beyond the singularity by the mathematical operation of multiplying both sides of (197) by $a$ (or more generally by some power of the metric determinant).

The idea of resolving a singularity by mulptiplying Einstein equations with some power of the determinant of the metric is not new. Einstein himself had proposed this idea with his collaborator Rosen in 1935 (for which they credited this idea to Mayer) PhysRev.48.73 . Ashtekar used a similar trick in his method of “new variables” to develop an equivalent Hamiltonian formulation of GR PhysRevD.36.1587 . It is also proposed by Stoica that the equations obtained after multiplying the usual Einstein equations by some power of the metric determinant are actually more fundamental than the usual Einstein equations Stoica:2013wx ; doi:10.1142/S0219887814500418 ; Stoica:2012my ; Stoica:2012gb ; Stoica:2011xf ; Stoica:2011nm ; Stoica:2014tpa ; Stoica:2015yfa ; Stoica:2015hba . In this sense, we argue that our spacetime with singularities due to the metric becoming degenerate ($a=0$) is a legitimate solution of GR.

IX.2.2 Singularities do not cause problems

While singularities are natural and inevitable in solutions to Einstein’s equations, we must discuss the consequences they bring to this calculation.

Will the singularities cause serious problems? At least in our case we do not feel they cause problems. To see this, we investigate how the singularities affect the propagation of the field modes in our toy model (100).

In this toy model, the spacetime have singularities at the hypersurfaces

[TABLE]

Using the relation $x_{m}=-x_{-m}$, which is evident from (116) and (117), the asymptotic mode solution (126) becomes

[TABLE]

At the singularities (198), the terms $\sin 2m(\Omega t+\mathbf{K}\cdot\mathbf{x})$ of (199) are all zero and thus we have

[TABLE]

So $u_{\mathbf{k}}$ is normal at singularities which shows that the field can naturally pass the singularities without problems.

One might still worry about the divergences of the time or spatial derivatives of $u_{\mathbf{k}}$ at the singularities (198). However, these divergences arise from those small high frequency corrections (terms inside the summation symbol $\Sigma$ of (199)) with frequencies $2m\Omega$ which are much higher than our effective QFT’s cutoff $\Lambda$. When looking at low energy scales ($\leq\Lambda$), $u_{\mathbf{k}}$ behaves the same as the mode solution (95) when the background spacetime is flat; only when looking at high energy scales ($\geq\sqrt{G}\Lambda^{2}$) which are much higher than the cutoff scale $\Lambda$, those small high frequency corrections are noticeable.

In this sense, the singularities do not cause problems at the observable low energy regime—after all, the singularities only appear (and immediately disappear) above Planck energy scales, which should not affect the low energy physics whose energy scale is far below Planck.

IX.3 Similarity of effects of vacuum energy in non-gravitational system and gravitational system

Vacuum fluctuations and their associated vacuum energies are direct consequences of the Heisenberg’s uncertainty principle of quantum mechanics. Although it is still controversial PhysRevD.72.021301 , various observable effects are often ascribed to the existence of vacuum energies and have been experimentally verified, which strongly suggests the reality of vacuum fluctuations. These vacuum fluctuation effects include the spontaneous emission scully1997quantum , the Lamb shift PhysRev.72.241 , the anomalous magnetic moment of the electron PhysRev.73.416 ; PhysRevA.83.052122 and the Casimir effect Casimir:1948dh ; PhysRevLett.78.5 ; PhysRevLett.81.4549 ; PhysRevLett.88.041804 . The reality of the vacuum energy associated to the spontaneous symmetry breaking of electroweak theory has also been confirmed by the discovery of the Higgs boson at the LHC Aad20121 ; Chatrchyan201230 .

If we assume that the vacuum fluctuations do exist as evidenced by the above listed observable effects, then according to the equivalence principle, the associated vacuum energies would gravitate as well as all other forms of energy. This has been experimentally demonstrated by, for example, the gravitational test of Lamb shift energy Polchinski:2006gy ; Masso:2009zd ; Braginskii:1971tn . The gravitational property of Casimir energy has not been tested experimentally, but has been demonstrated theoretically with the conclusion that it does gravitate according to equivalence principle Fulling:2007xa ; Milton:2007ar ; Milton:2007hd ; Milton:2008ks .

However, in the literature, the value of vacuum energy density is usually thought to play a different role in non-gravitational systems and in gravitational systems. The actual value of the vacuum energy density is generally regarded as irrelevant in non-gravitational contexts based on the argument that only energy differences from the vacuum are measurable; while when gravity is present, the actual value of the energy matters, not just the differences, since the source for the gravitational field is the entire energy momentum tensor that its large value may be potentially disastrous.

We argue differently in this section with the following points: (i) the value of vacuum energy density can also be relevant in non-gravitational contexts; (ii) the huge value of vacuum energy density is not a direct observable and that it is not disastrous in a theory of gravity. Moreover, there is essentially no difference between the roles played by vacuum energy in non-gravitational systems and in gravitational systems. In other words, although technically more complicated when gravity is included, the gravitational effect of the vacuum energy on spacetime metric is intrinsically the same as its effect on material bodies when gravity is excluded.

IX.3.1 Value of vacuum energy is relevant in Casimir effect

Let us first consider the Casimir effect. The Casimir force is usually derived by calculating the change in vacuum energy due to the presence of the conducting plates, which acts as mirrors to reflect electromagnetic waves (We will call them mirrors in the following). This derivation is straightforward, but loses some important physical details about what is going on in the system GONZALEZ1985228 ; PhysRevA.38.1621 . Due to quantum fluctuations, the zero point fields constantly impinge on both sides of the mirror and then reflect back, which transmit momentum to the mirror and thus result in forces on both sides of the mirror. The Casimir stress (force per unit area) is just the difference between the pressure exerted by the electromagnetic field vacuum from inside and outside

[TABLE]

where we have set that the two parallel mirrors are normal to the $z$ axis. Since the vacuum fluctuations between the two mirrors are different from the vacuum fluctuations outside, the expectation values of $T_{zz}^{\mathrm{inside}}$ and $T_{zz}^{\mathrm{outside}}$ would be different and thus gives a net average force. Although both $T_{zz}^{\mathrm{inside}}$ and $T_{zz}^{\mathrm{outside}}$ are divergent, this average force is finite since the quartic divergent Minkowski zero point fluctuations are canceled after the subtraction in (201) and one obtains the well known Casimir stress PhysRevA.38.1621

[TABLE]

Thus the effect of the value of zero point energy disappears in the calculations. It is for this reason that although the Casimir effect is usually regarded as evidence of the reality of zero point energy, the actual value of its energy density is thought to be irrelevant in this effect.

However, the value of zero point energy density does have an effect. Note that (202) only gives the expectation value of the Casimir stress $S$, but $S$ is never a constant, it fluctuates. That’s because the amount of momentum carried by the zero point fields which impinge on both sides of the mirror is constantly fluctuating due to the fact that the vacuum is not an eigenstate of the $zz$ component of the stress energy tensor $T_{zz}$. The magnitude of the fluctuation of each $T_{zz}$ is large and diverges as the same order of the vacuum energy density $\left\langle T_{00}\right\rangle$. For a perfect mirror, since the fields on the two sides fluctuate independently of each other, the mean-squared stresses on the two sides simply add, resulting in the magnitude of the fluctuation of the net stress also diverges as

[TABLE]

For more realistic imperfect mirrors which become transparent for frequencies higher than its plasma frequency $\Lambda$, the $\left\langle T_{00}\right\rangle$ in (203) contains contributions only from field modes of frequencies lower than $\Lambda$ and the mean squared value of the net stress $S$ goes as

[TABLE]

The plasma frequency $\Lambda$ in (204) acts as an effective cutoff which depends on the microstructure of the mirror. It is similar but distinct from the effective QFT’s cutoff $\Lambda$ in (75), which depends on the microstructure of spacetime.

Therefore, the value of zero point energy density is still physically significant even in non-gravitational system. Its value appears in (203) and (204) to characterize the strength of Casimir stress fluctuation, which implies that the net Casimir stress is constantly fluctuating with huge magnitudes around its small mean value (202). Due to this huge fluctuation, at almost any instant, the magnitude of the stress at each single point of the mirror is as large as the value of the zero point energy density.

However, this effect is strong only at small scales. Its measurable effect becomes small at larger scales. In practice, the measurements must be taken over some finite time interval $T$ and some finite surface area of order $l^{2}$. More precisely, what the force detector measures is the time and surface average

[TABLE]

where the averaging function $f$ satisfies

[TABLE]

The exact shape of the averaging function depends on the measuring apparatus. On physical grounds one can choose $f$ to be a single peak over a time interval $T$ comparable to the experimental resolving time and over a spatial region of area $l^{2}$ comparable to the resolution of the measuring device. Although the magnitude of the fluctuations of the net stress $S$ is formally infinite as shown in (203), the magnitude of the measurable fluctuations of its average $\bar{S}$ is finite. This is because the effect of the vacuum fluctuations at small scales is significantly weakened when averaging over larger scales. The calculations have been done by G Barton in 0305-4470-24-5-014 with the conclusion that, for the realistic case where $l\ll cT$, the mean squared deviation

[TABLE]

where the “$\mathrm{constant}$” here is a pure number as could have been foreseen on dimensional grounds. The above equation (207) shows that $\left\langle\Delta\bar{S}^{2}\right\rangle$ increases as $T$ decreases, which means that the better the measuring device, the stronger fluctuation due to the effect of the value of the zero point energy density can be measured. And in principle, using a perfect instantaneous measuring device ($T\to 0$), one can measure the infinite fluctuations of the Casimir stress on a perfect mirror due to the infinite value of zero point energy density. In practice, however, $\left\langle\Delta\bar{S}^{2}\right\rangle$ is too small to be measured for a real force detector whose resolving time $T$ is too large 0305-4470-24-5-014 .

IX.3.2 Effect of vacuum energy on the motion of mirrors

The value of zero point energy density also has effects on the dynamic motion of small material bodies. Imagine that we place a single mirror of very small size in the vacuum and then release it. The mirror would experience a fluctuating force exerted by the quantum field vacuum and starts to move. The equation of motion of the mirror, which is called quantum Langevin equation, can be generally described by

[TABLE]

where $X$ is the mirror’s position, $\phi$ represents the field interacting with the mirror which is usually taken to be a scalar field for simplicity and we have set the mirror’s mass $M=1$ for convenience. The average force in this case would be zero because of symmetry

[TABLE]

and similar to the Casimir stress fluctuation (203), the force here also undergoes wild fluctuations with a magnitude

[TABLE]

The mathematically infinite fluctuating force $F$ gives infinite instantaneous accelerations of the mirror through (208). Similar to the case of infinite Casimir stress fluctuation (203), this infinite fluctuating force and infinite instantaneous acceleration make sense since they are also only significant at very small scales and will not result in infinite fluctuation of the mirror’s position at observable larger scales. In fact, the mirror would oscillate back and forth with very high speeds, but its range of motion is still small PhysRevD.89.085009 ; PhysRevD.92.063520 ; Gour:1998my ; Jaekel:1992ef ; Stargen:2016euf .

More precisely, suppose that the mirror is initially located at $X(0)=0$ with velocity $\dot{X}(0)=0$ and is then released at $t=0$. The magnitude of its acceleration $\ddot{X}(t)$ and velocity $\dot{X}(t)$, which can be characterized by the quantity $\left\langle\ddot{X}^{2}(t)\right\rangle$ and $\left\langle\dot{X}^{2}(t)\right\rangle$, is large. But, the magnitude of the range of the mirror’s fluctuating motion, which can be characterized by the observable mean squared displacement $\left\langle X^{2}(t)\right\rangle$, is still small.

In this sense, the value of vacuum energy density is still relevant even in non-gravitational physics. This value appears in the equation (210) to characterize the strength of the force fluctuations acting on the mirror at small scales and it may have small observable effects at larger scales such as diffusions predicted in Gour:1998my ; Jaekel:1992ef ; Stargen:2016euf .

IX.3.3 Analogies between the motion of mirror and the motion of $a(t,\mathbf{x})$

Although technically more complicated in gravity, the basic dynamic equation of motion (41) satisfied by the scale factor $a(t,\mathbf{x})$ is in fact very similar to the equation of motion (208) satisfied by the mirror’s position $X(t)$. Consider only the contribution from the massless scalar field $\phi$, (41) is just the following same form as the equation (208)

[TABLE]

where

[TABLE]

Also, the average of the fluctuating force $F\left(a,\dot{\phi}\right)$ is zero due to symmetry

[TABLE]

and its magnitude of fluctuation

[TABLE]

The above two statistical properties (213) and (214) satisfied by the “force” driving the “motion” of the scale factor $a$ are the same with the statistical properties (209) and (210) satisfied by the force driving the motion of the mirror. In this sense, the role played by the value of the vacuum energy density in gravitational system is similar to its role in non-gravitational system.

Concretely speaking, the vacuum energy density results in large instantaneous acceleration $\ddot{X}$ and velocity $\dot{X}$ of the mirror, but the observable position fluctuations of the mirror, which can be characterized by the quantity $\left\langle X^{2}\right\rangle$, is not large. Analogously, the vacuum energy density results in the large instantaneous “acceleration” $\ddot{a}$ and “velocity” $\dot{a}$ of the scale factor, but the observable physical distance defined by (30), whose value is determined by the quantity $\left\langle a^{2}\right\rangle$, is also not large. These properties about $a(t,\mathbf{x})$ are evident from the solutions (51), (58) and (75), from which we can see that the quantities $\left\langle\ddot{a}^{2}\right\rangle$ and $\left\langle\dot{a}^{2}\right\rangle$ are as large as $\left\langle T_{00}\right\rangle^{2}$ and $\left\langle T_{00}\right\rangle$ respectively, while the magnitude of the quantity $\left\langle a^{2}\right\rangle$ is on the order $1$.

In this sense, the role played by vacuum energy in gravitational system is similar to its role in the non-gravitational mirror systems—it appears both at (210) and (214) to show the strongness of vacuum fluctuations at microscopic scales (for mirrors, microscopic means atomic scale; for gravity, microscopic means Planck scale) and their observable effects are both small at macroscopic scales.

By this same kind of mechanism, the violent gravitational effect produced by the vacuum energy density is confined to Planck scales, and its effect at macroscopic scales—the accelerating expansion of the Universe, due to the weak parametric resonance is so small that, it is only observable after accumulations on the largest scale—the cosmological scale.

Acknowledgements

We would like to thank Andrei Barvinsky, Gordon W. Semenoff, Mark Van Raamsdonk, Philip C. E. Stamp, Daniel Carney, Yin-Zhe Ma, Daoyan Wang, Michael Desrochers, Giorgio Torrieri, Niayesh Afshordi, Roland de Putter, Jérôme Gleyzes and Olivier Doré for helpful discussions and criticisms. We especially thank Michael Desrochers for helping revising the manuscript. This research is partially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), the Templeton Foundation, and the Canadian Institute for Advanced Research. Zhen Zhu would also like to thank Philip C. E. Stamp for financial support.

Appendix A Real Massless Scalar Field

In this appendix we give the calculation details about how the quantum vacuum fluctuates all over the spacetime by using the massless scalar field (15) as an example.

We first define the covariance of the energy density operator at two spacetime points $x=(t,\mathbf{x})$ and $x^{\prime}=(t^{\prime},\mathbf{x}^{\prime})$

[TABLE]

where the curly bracket $\{\}$ in (A) is the symmetrization operator which is defined as, for any two operators $A$ and $B$,

[TABLE]

Inserting (15) and (18) into (A) gives the following result

[TABLE]

where $\Delta t=t-t^{\prime}$ and $\mathbf{\Delta x}=\mathbf{x-x^{\prime}}$ are time and space separation of the two spacetime points $x$ and $x^{\prime}$.

If $x$ and $x^{\prime}$ are timelikely separated, we can find a reference frame to set $\Delta x=|\mathbf{\Delta x}|=0$. In this case, evaluation of the integral in (LABEL:correlation_function_0) for a high frequency cutoff $|\mathbf{k}|=\Lambda$ gives

[TABLE]

If $x$ and $x^{\prime}$ are spacelikely separated, we can find a reference frame to set $\Delta t=0$. In this case, evaluation of the integral in (225) for a high frequency cutoff $|\mathbf{k}|=\Lambda$ gives

[TABLE]

As $\Delta t$ and $\Delta x$ goes to [math], both (218) and (219) reduces to the variance of the energy density,

[TABLE]

We then investigate the Pearson product-moment correlation coefficient

[TABLE]

where

[TABLE]

The correlation coefficient $\rho_{x,x^{\prime}}$ shows by its magnitude the strength of correlation between two random variables. $\rho_{x,x^{\prime}}$ is positive if the energy density $T_{00}$ at $x$ and $x^{\prime}$ are most possibly lying on the same side of the vacuum expectation value $\left\langle T_{00}\right\rangle=\Lambda^{4}/(16\pi^{2})$. Thus a positive correlation coefficient $\rho_{x,x^{\prime}}$ implies the energy density at $x$ and $x^{\prime}$ tend to be simultaneously greater than, or simultaneously less than the expectation value. Similarly, a negative $\rho_{x,x^{\prime}}$ implies the energy density tend to lie on opposite sides of the expectation value. We will call the energy density $T_{00}$ at $x$ and $x^{\prime}$ are positively correlated if $\rho_{x,x^{\prime}}>0$ or negatively correlated (anticorrelation) if $\rho_{x,x^{\prime}}<0$.

Because of transnational invariance, $\rho_{x,x^{\prime}}$ is only dependent on the temporal and spatial separation $\Delta t=t-t^{\prime},\mathbf{\Delta x}=\mathbf{x}-\mathbf{x}^{\prime}$. For the real massless scalar field (15), the behavior of the correlation coefficient $\rho_{x,x^{\prime}}$ as a function of temporal separation $\Lambda\Delta t$ for the case of $\Delta x=0$ and as a function of spatial separation $\Lambda\Delta x$ for the case of $\Delta t=0$ are plotted in FIG. 7 and 8 respectively.

In the temporal direction, i.e. the case of $\Delta x=0$ (Fig. 7), the correlation coefficient goes quickly from $1$ down to around $-0.9$ in a time scale around $\Delta t=1.9/\Lambda$ and then goes up to $0.7$ in a time scale around $\Delta t=3.8/\Lambda$ and then goes down and up alternatively from positive values to negative values with decreasing amplitudes. It roughly oscillates as $-\cos(2\Lambda\Delta t)/(\Lambda\Delta t)^{2}$ with a period $\pi/\Lambda$ as $\Delta t$ is large. Thus at the extremely small time scales $\Delta t\sim 1.9/\Lambda$, ($\Lambda\to+\infty$), the energy density are strongly anticorrelated. In other words, if at some time the value of the energy density is larger than its expectation value, for example, by an amount of $0.82\left\langle T_{00}\right\rangle$, after a short time $\Delta t=1.9/\Lambda$, its value is most likely to be smaller than the expectation value, for example, by an amount of $0.74\left\langle T_{00}\right\rangle$. The difference is $1.56\left\langle T_{00}\right\rangle$ only after such a short time.

In the spatial direction, i.e. the case of $\Delta t=0$ (Fig. 8), the correlation coefficient goes quickly from $1$ down to around $-0.14$ in a length scale around $\Delta x=3.24/\Lambda$ and then goes up to $0.03$ in a length scale around $\Delta x=5.4/\Lambda$ and then goes down and up alternatively from positive values to negative values with decreasing amplitudes. Compared to the temporal direction, the decay in the oscillation amplitude of the correlation coefficient is faster in spatial direction. It roughly oscillates as $2\cos(2\Lambda\Delta x)/(\Lambda\Delta x)^{4}$ with a period $\pi/\Lambda$ as $\Delta x$ is large. These properties show that the strength of the correlation between energy densities at close range in spatial direction is not as strong as in the temporal direction. For larger spatial separations, $\rho_{x,x^{\prime}}$ approaches zero and the vacuum energy density $T_{00}$ at different $x$ and $x^{\prime}$ fluctuate independently. These properties result in extreme spatial inhomogeneities of the quantum vacuum which can be characterized by the quantity $\Delta\rho^{2}$ defined by (22) in section III.

The quantity $\Delta\rho^{2}$ is related to $\rho_{x,x^{\prime}}$ by

[TABLE]

The behavior of $\Delta\rho^{2}$ has been plotted in FIG. 1.

Next we calculate the $\chi(\Delta t)$ defined by (V.1) in section V.1. Wick expansion of (V.1) gives

[TABLE]

where the correlation function can be calculated directly by inserting (15)

[TABLE]

Plugging (225) into (224) gives the following result

[TABLE]

The behavior of $\chi(\Delta t)$ has been plotted in FIG. 2. It is closely related to the correlation coefficient $\rho_{x,x^{\prime}}$ as a function of time difference $\Delta t$ in the case $\Delta x=0$ (FIG. 7).

Next we derive the equation (79) in section V.1. First, $\Omega^{2}(t,\mathbf{0})$ can be expanded as

[TABLE]

Specially, the vacuum state $|0\rangle$ is an eigenstate of the operator coefficients of the first two terms in the above expression (A). If $\mathbf{k}\neq\mathbf{k}^{\prime}$, the eigenvalues of the operator coefficients of the first two terms are zero. Thus in this case, the first two terms have to both take zero values. If $\mathbf{k}=\mathbf{k}^{\prime}$, the second term is zero since in this case $\omega=\omega^{\prime}$ and thus the factor $\sin(\omega-\omega^{\prime})t=0$. So only the first term survives and gives the expectation value of $\Omega^{2}(t,\mathbf{0})$:

[TABLE]

For the operator coefficients of the last two terms in the expression (A), the vacuum state $|0\rangle$ is not an eigenstate. So the last two terms are constantly fluctuating, and the time varying of $\Omega^{2}$ comes from these two terms.

After some algebraic manipulations, (A) can be rewritten as the form of (79) for the vacuum state $|0\rangle$, where

[TABLE]

[TABLE]

Evaluating the above integrals gives the expectation values

[TABLE]

and their fluctuations

[TABLE]

The above expression (232) gives the power spectrum density of the varying part of $\Omega^{2}(t,\mathbf{0})$ (except the constant $\Omega_{0}^{2}$ part), which has been plotted in FIG. 3.

Appendix B Wigner-Weyl Description of Quantum Mechanics and Numeric simulations

This chapter explain the principle of the numeric calculations in the main text. Same as the numeric part in the main text, we set $G=1$ in this section. Wigner functions and Weyl transforms of operators offer a formulation of quantum mechanics that is equivalent to the standard approach given by the Schrödinger equation. The Wigner distribution function is a quasi distribution function in the phase space. For a particular quantum wave function $\psi(x)$, its Wigner function is defined as

[TABLE]

The Weyl transform of an quantum operator $\hat{A}$ is defined as

[TABLE]

Then the expectation value of the operator $\hat{A}$ under the state $\psi(x)$ can be written as

[TABLE]

These two transformations give the Wigner-Weyl discription for quantum mechanics. The expectation values of physical quantities are obtained by averaging their Weyl transforms over phase space.

For a harmonic oscillator with frequency $\omega$ and $m=1$, the ground state Wigner function is a Gaussian distribution function for both $x$ and $p$

[TABLE]

We can easily check that the Weyl transform of an operator $H(\hat{x})$ ( or $H(\hat{p})$) is simply replaced the operator $\hat{x}$ by $x$ (or $\hat{p}$ by $p$). Other than that, another particular transform we are going to use in this write-up is

[TABLE]

We can see that the transform of the product does not necessarily equal to the product of transforms. In the following part we are going to get the general expression for the transform of the product.

Before that we notice that Weyl transform can be used to construct the original operator , i.e.

[TABLE]

Using this formula we can construct the transform of product of two states:

[TABLE]

Here we define

[TABLE]

We Taylor-expand $A(x+\frac{z_{1}}{2},p_{1})$ and $B(x+\frac{z_{2}}{2},p_{2})$ around $x$ and have

[TABLE]

and use the facts

[TABLE]

and

[TABLE]

Therefore, we can write the Weyl transform of operator $\hat{A}\hat{B}$ as

[TABLE]

The generalized FRW scale factor $a$ satisfies the equation

[TABLE]

in which

[TABLE]

Now we replace all the quantities by operators, assuming that operators still satisfy the previous equation

[TABLE]

with

[TABLE]

For a massless real scalar field, we can write it as

[TABLE]

in which

[TABLE]

are the generalized $\hat{x}$ $\hat{p}$ operators for each field modes.

We can write the Weyl transformation of the $\hat{\Omega}(t)$

[TABLE]

This expression is quadratic in $x_{k}$ and $p_{k}$, so if we apply it to (246), only $m+n\leq 2$ terms survive. Assuming $a(\{x_{k}\},\{p_{k}\},t)$ is the Weyl transform of operator $\hat{a}$, we have the equation for $a$ as

[TABLE]

The observed value $a$ is the average over Wigner function $W(\{x_{k}\},\{p_{k}\},t)$

[TABLE]

If the quantum field is in the ground state, then by (238)

[TABLE]

Local approximation

Generally the equation (255) depends on not only the value of $\Omega$ and $a$ on a particular phase space point $(x,p)$, but also on the neighboring values (i.e. derivatives). If our solution $a$ is ”smooth” enough in the phase space then we can neglect the last two derivative terms in the (255). It can be simplified to

[TABLE]

Assuming the length of the Universe is $L$. We can replace the integral by summations. For simplicity, we define

[TABLE]

The equation can be written as

[TABLE]

with

[TABLE]

Here $\vec{n}=(n_{1},n_{2},n_{3}),n_{1,2,3}\in\mathbb{Z}$ and $n=|\vec{n}|$. $\{\tilde{x}_{\vec{n}}\}$ $\{\tilde{p}_{\vec{n}}\}$ are random Gaussian variables with unit standard deviation. We can solve the equation for a randomly generated set of $\{\tilde{x}_{\vec{n}}\}$ and $\{\tilde{p}_{\vec{n}}\}$, and repeat. The result $a_{o}(t)$ is the average over all solutions as long as our sample size is big enough.

Appendix C Fourier transforms of the coefficients in (130)

In this appendix, we demonstrate the property of the spectrum of the coefficients in (130) given by (140), (142) and (144). Observing that the $\cos 2\Theta$, $\sin 2\Theta$ and $\tan\Theta$ in (139), (141) and (143) respectively can all be decomposed as Fourier series sum of the form $e^{i2n\Theta}$, where $n=\pm 1,\pm 2,\cdots$, we only need to analyze the spectrum of $e^{i2n\Theta}$.

For simplicity, we only analyze the time component Fourier transform of $e^{i2n\Theta}$. The spatial part has similar property. The phase angle $\Theta$ is determined by $\Omega$ through (97) while $\Omega$ is determined by (79). The power spectrum of $\Omega^{2}$ is given by (232) (illustrated in FIG. 3).

Calculation of the Fourier transform of $e^{i2n\Theta}$ exactly based on (79) is complicated. For simplicity, we assume that $\Omega$ taking the following simple form which is similar to (77)

[TABLE]

where $\gamma$ take the peak value of the power spectrum (232) which is around $\sim 1.7\Lambda$ (see FIG. 3) and $h<1$ to make sure that $\Omega>0$.

Then we have

[TABLE]

Using the Jacobi-Anger expansion we have

[TABLE]

where $J_{m}$ is the $m$th Bessel function of the first kind.

As $|m|\to\infty$, we have

[TABLE]

which drops faster than the exponential. Therefore, the Fourier transform of $e^{i2n\Theta}$ is centered around $2n\Omega_{0}$.

To estimate the magnitude of the Fourier coefficients of $e^{i2n\Theta}$ around zero frequency, we evaluate the Bessel function for

[TABLE]

In this case, the zero component Fourier coefficient is asymptotic to (see watson1995treatise )

[TABLE]

since $\nu$ is determined by $h=\operatorname{sech}\nu<1$ that we always have $\nu-\tanh\nu>0$.

When calculating the Fourier transform of $e^{i2n\Theta}$ exactly based on (79), the spectrum becomes continuous instead of discrete. But the distribution of the spectrum should be similar.