Study of Nonlinear Evolution of Spacetime Fluctuations in Quantum Gravity Inflation for Deriving Primordial Spectrum

Abstract

We study the evolution of quantum fluctuations of gravity around an inflationary solution in renormalizable quantum gravity, in which the initial scalar-fluctuation dominance is shown by the background-free nature expressed by a special conformal invariance. Inflation ignites at the Planck scale and continues until spacetime phase transition occurs at a dynamical scale of about ${10}^{17}$ GeV. We show that during inflation, the initially large scale-invariant fluctuations reduce in amplitude to the appropriate magnitude suggested by tiny CMB anisotropies. The goal of this research is to derive the spectra of scalar fluctuations at the phase transition point, that is, the primordial spectra. A system of nonlinear evolution equations for the fluctuations is derived from the quantum gravity effective action. The running coupling constant is then expressed by a time-dependent average following the spirit of the mean field approximation. In this paper, we determine and examine various nonlinear terms, not treated in previous studies such as the exponential factor of the conformal mode. These contributions occur during the early stage of inflation when the amplitude is still large. Moreover, in order to verify their effects concretely, we numerically solve the evolution equation by making a simplification to extract the most contributing parts of the terms in comoving momentum space. The result indicates that they serve to maintain the initial scale invariance over a wide range beyond the comoving Planck scale. This is a challenge toward the derivation of the precise primordial spectra, and we expect in the future that it will lead to the resolution of the tensions that have arisen in cosmology.

1. Introduction and Summary

Inflation is an unbelievable idea proposed by Guth [1], Sato [2], and Starobinsky [3] to solve the horizon and flatness problems, which claims that there was a period of exponential expansion in the universe before the big bang1. Remarkably, if we genuinely embrace this concept, a significant portion of the universe we observe today originated from a region even smaller than the Planck length.

Going over the Planck-scale wall requires quantization of gravity with not only renormalizability but also background freedom. Defying the common perception in this research area, we have proposed asymptotically background-free quantum gravity [5,6,7,8,9] as a theory with such properties and also a realistic inflationary dynamics. Its distinctive feature is that diffeomorphism invariance in the ultraviolet (UV) limit is described as a special conformal invariance [10,11,12,13,14,15,16]. Since it is a gauge symmetry, it represents background freedom in which all different conformally-flat spacetimes are gauge equivalent. This symmetry is called BRST conformal invariance,2 which was first found in two-dimensional quantum gravity [17,18,19,20] and then extended to four dimensions. The existence of this symmetry makes it possible to describe the world far beyond the Planck scale where scalar fluctuations dominate and no tensor ones.

Perturbation theory is formulated by introducing a dimensionless coupling constant t that represents deviation from the conformal invariance. The asymptotic background freedom represents that nonperturbative fluctuations of the conformal mode in the gravitational field become dominant in the UV limit of $t\to 0$. In this paper, we investigate how the universe evolves and deviates from such a scale-invariant state. The moment when the conformal invariance is completely broken at a novel energy scale of quantum gravity, ${\Lambda }_{\mathrm{QG}}$, is called the “spacetime phase transition”. If we set this scale below the Planck scale, then the theory has an inflationary solution with the expansion time constant of order of the Planck mass [21,22,23].

The goal of this study is to clarify how perturbations, or spacetime fluctuations, around the inflationary solution evolve, and to derive the primordial spectra that are the initial conditions of the current universe. It can be said that in the universe before inflation begins, quantum fluctuations of gravity are large and thus spacetime is not substantially fixed. As inflation begins, the fluctuations gradually decrease in amplitude, eventually reaching the magnitude suggested by observations of tiny anisotropies in the cosmic microwave background radiation (CMB) [24,25] until the time of the spacetime phase transition. Since fluctuations in gravity represent fluctuations in time and distance, the reduction of them implies that the universe becomes a real world where time and distance can be measured accurately3.

In order to obtain the primordial spectrum, various nonlinear effects in the evolution of the fluctuations must be taken into account properly. The running coupling constant is one of these nonlinear effects, which is incorporated into the equations of motion by approximating it as a time-dependent mean field, as done in previous studies [22,23]. In this way, we show that the amplitude of the fluctuations reduces during inflation.

To determine the spectral pattern more accurately, we need to evaluate contributions from nonlinear terms, such as the exponential factor of the conformal mode existing in the Einstein–Hilbert action. The fluctuations dealt with here are of sizes that momentum dependence disappears as soon as inflation begins, but the nonlinear effects cannot be ignored in the early stage where the fluctuations are still large4. One of the purposes of this research is to determine a system of nonlinear evolution equations of the fluctuations involving such effects.

Since it is a rather complicated equation system, we here simplify it by extracting most contributing parts of the nonlinear terms in comoving momentum space so that we can solve it numerically and quantify their effects, as will be described in Section 4 and Section 5. One of the results indicating that there is a large effect on the evolution of fluctuations is shown in Figure 1 and Figure 2, together with previous results from linear equations. Consequently, we confirm that the nonlinear terms have expected effects by which the initial scale invariance is maintained up to relatively-high momentum regions, as inferred in the previous works [23].

In determining the spectrum, the presence of the physical infrared cutoff ${\Lambda }_{\mathrm{QG}}$ plays an essential role. The sharp falloff in low multipole components of the CMB angular power spectrum is explained by this scale [21]. Its wavelength is given by a comoving scale of the correlation length ${\xi }_{\Lambda }=1/{\Lambda }_{\mathrm{QG}}$. This means that most of the universe we see today was within the size of ${\xi }_{\Lambda }$ before inflation.

Quantum gravity inflation suggests that the observed CMB anisotropy spectrum contains real quantum gravity effects. The guiding principles here are clear: diffeomorphism invariance and renormalizability. Thus the model gives life to the concept of inflation, aligning with Starobinsky’s perspective [3]. Moreover, we expect that the primordial spectrum derived from the first principle will settle the Hubble tension [27].

2. Renormalizable and Asymptotically Background-Free Quantum Gravity

The problems with Einstein’s theory of gravity are summarized in unrenormalizability and the existence of singularities. Therefore, a finite UV cutoff is usually introduced in the Planck scale to avoid such problems. This corresponds to thinking of spacetime discretized in the Planck length, and some researchers regard the length as a quantum of spacetime. However, introducing the UV cutoff breaks diffeomorphism invariance.

If we consider diffeomorphism invariance as a first principle, then we have to bring the cutoff to infinity, i.e., take the continuum limit. It implies that the theory is renormalizable. Being renormalizable is equivalent to having no singularity. The reason why the existence of singularities cannot be denied probabilistically is due to the fact that the Einstein–Hilbert action becomes finite for such singular spacetime configurations, because it does not contain the Riemann tensor $R_{\mu \nu \lambda \sigma }$ which controls the magnitude of curvature. A renormalizable quantum theory of gravity requires an action involving the Riemann tensor squared with correct sign. Since such an action is positively divergent for singularities, they are eliminated as unphysical5. Furthermore, since the action becomes dimensionless, the corresponding coupling constant is also dimensionless and renormalizable [28,29,30,31].

Despite these excellent properties, it is generally known that another problem arises in such a theory. This is the so-called ghost problem. However, note here that the existence of ghost modes itself is not a problem. In fact, Einstein’s theory of gravity has a ghost mode due to the indefiniteness of the Einstein–Hilbert action. Owing to that, there are non-trivial solutions such as the Friedmann solution even though we are considering an equation where the total Hamiltonian vanishes. If the action was positive definite, there should only be a trivial vacuum solution. Thus, the ghost mode is an indispensable element to make the Hamiltonian zero, i.e., to preserve diffeomorphism invariance. However, it causes problems if it appears locally as a real object. A problem with early renormalizable theories using curvature squared is that the ghost mode appears as a physical particle state in the UV limit.

Considering why such a problem occurs, it turns out that there is a problem in the method of perturbation expansion. The conventional perturbation expansion has been carried out in the graviton picture propagating in the flat spacetime of $R_{\mu \nu \lambda \sigma }=0$, but this expansion method cannot make the ghost mode unphysical. To begin with, assuming that flat spacetime appears asymptotically implies that the closer we approach the center of the black hole, the weaker the gravitational field. Thus, this is a physically unnatural setting. Even theoretically this perturbation method exhibits inconsistent UV behavior (see also Note 6).

How should we perform perturbation expansion to correctly describe features of the world beyond the Planck scale? A solution to this problem can be found by considering what is required in the early universe. The important fact here is that inflation is given by a spacetime configuration in which the Weyl tensor

$$\begin{array}{ccc}\hfill C_{\mu \nu \lambda \sigma }& =& R_{\mu \nu \lambda \sigma }-\frac{1}{2}{g}_{\mu [\lambda }{R}_{\sigma ]\nu }+\frac{1}{2}{g}_{\nu [\lambda }{R}_{\sigma ]\mu }+\frac{1}{6}{g}_{\mu [\lambda }{g}_{\sigma ]\nu }R\hfill \end{array}$$

vanishes, where $a_{[\mu }{b}_{\nu ]}=a_{\mu }{b}_{\nu }-a_{\nu }{b}_{\mu }$. Therefore, we perform perturbation expansion around the conformally-flat spacetime with $C_{\mu \nu \lambda \sigma }=0$.

This suggests that the conformal mode of the gravitational field, crucial for determining distances, receives special treatment. Here, we extract it in an exponential factor, ensuring its positivity

$$\begin{array}{c}\hfill g_{\mu \nu }=e^{2\varphi }{\overline{g}}_{\mu \nu },\end{array}$$

where the scalar-like field $\varphi$ is called the conformal-factor field. The remaining mode ${\overline{g}}_{\mu \nu }$ is expanded by the traceless tensor field $h_{\mu \nu }$ as

$$\begin{array}{c}\hfill {\overline{g}}_{\mu \nu }={\eta }_{\mu \lambda }{\left(e^h\right)}_{\nu }^{\lambda }={\eta }_{\mu \lambda }\left({\delta }_{\nu }^{\lambda }+h_{\nu }^{\lambda }+\frac{1}{2}{h}_{\sigma }^{\lambda }{h}_{\nu }^{\sigma }+\cdots \right),\end{array}$$

(1)

where $h_{\mu }^{\mu }=0$. Raising and lowering the legs of the traceless tensor field is done with the flat background metric ${\eta }_{\mu \nu }=(-1,1,1,1)$ as $h_{\mu \nu }={\eta }_{\mu \lambda }{h}_{\nu }^{\lambda }$. The coordinates of the flat metric are denoted as $x^{\mu }=(\eta ,\mathbf{x})$, and are henceforth called the comoving coordinate, following cosmology.

We introduce a coupling constant t as a dimensionless parameter that controls the above perturbation expansion. Noting that the field strength of the traceless tensor field is the Weyl tensor, the perturbation theory is defined by introducing the inverse of $t^2$ before the Weyl tensor squared with correct sign. The action of the whole system including it is given by [5,6,7,8,9]

$$\begin{array}{c}\hfill I=\int d^{4}x\sqrt{-g}\left\{-\frac{1}{t^2}{C}_{\mu \nu \lambda \sigma }^2-b{G}_4+\frac{1}{\hslash }\left(\frac{1}{16\pi G}R+{\mathcal{L}}_{\mathrm{M}}\right)\right\}.\end{array}$$

The first term is called the Weyl action and is conformally invariant, $\sqrt{-g}{C}_{\mu \nu \lambda \sigma }^2={\overline{C}}_{\mu \nu \lambda \sigma }^2$, where the quantity with the bar is defined by ${\overline{g}}_{\mu \nu }$. The second $G_4=R_{\mu \nu \lambda \sigma }{R}^{\mu \nu \lambda \sigma }-4R_{\mu \nu }{R}^{\mu \nu }+R^2$ is the Euler density and its volume integral is also conformally invariant, where b is not an independent coupling constant because this term does not contain a kinetic term, which is expanded in t. The third is the Einstein–Hilbert action and ${\mathcal{L}}_{\mathrm{M}}$ denotes the matter field actions which are conformally invariant in the UV limit. The cosmological constant term is ignored here because its contribution is negligible in the early universe. When performing perturbation calculations, we usually redefine the field as $h_{\mu \nu }\to t{h}_{\mu \nu }$ in the expansion formula (1), but here we proceed without redefining because we consider even regions where the coupling constant becomes large.

The coupling constant t represents the degree of deviation from the conformal flatness. The significant point of this expansion method is that in the UV limit of $t\to 0$, the traceless tensor field becomes small, but the conformal-factor field $\varphi$ remains fluctuating nonperturbatively. Thus, there is no picture of particles propagating in flat spacetime. The partition function under this expansion is given by [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20]:

$$\begin{array}{c}\hfill \int {\left[dg\right]}_g{e}^{iI\left(g\right)}=\int {\left[d\varphi dh\right]}_{\eta }{e}^{iS(\varphi ,\overline{g})+iI\left(g\right)},\end{array}$$

where the $\varphi$-dependent action S is the Wess–Zumino action [32] for conformal anomaly [33,34,35,36]. Unlike the conventional perturbation expansion around flat spacetime, S is necessary to preserve diffeomorphism invariance when replacing the path integral measure from the invariant ${\left[dg\right]}_g$ to the commonly used ${\left[d\varphi dh\right]}_{\eta }$ defined on the flat background6. This is a quantity that makes up the running coupling constant introduced later, which represents the breaking of conformal invariance. Although the word anomaly is used, it is not physically anomalous. If formulated using dimensional regularization that preserves diffeomorphism invariance manifestly, the information of S is automatically contained between 4 and $D(<4)$ dimensions [6,7,8,9,31,37,38,39,40,41,42].

The Wess–Zumino actions are responsible for fourth-derivative dynamics of the conformal-factor field $\varphi$. The most important term among them is the Riegert action [10], which remains even in $t\to 0$7:

$$\begin{array}{c}\hfill S_{\mathrm{R}}=\int d^{4}x\left\{-\frac{b_1}{{\left(4\pi \right)}^2}\left(2\varphi {\overline{∆}}_4\varphi +{\overline{E}}_4\varphi +\frac{1}{18}{\overline{R}}^2\right)\right\},\end{array}$$

(2)

where $E_4=G_4-2{\nabla }^{2}R/3$ and ${\nabla }^2={\nabla }^{\mu }{\nabla }_{\mu }$. The differential operator $\sqrt{-g}{∆}_4$, which is conformally invariant for scalars, is defined by ${∆}_4={\nabla }^4+2R^{\mu \nu }{\nabla }_{\mu }{\nabla }_{\nu }-2R{\nabla }^2/3+{\nabla }^{\mu }R{\nabla }_{\mu }/3$. The coefficient is expanded as

$$\begin{array}{c}\hfill b_1=b_{c}B,B=1-{\gamma }_1\frac{t^2}{4\pi }+o\left(t^4\right).\end{array}$$

(3)

The lowest is given by $b_c=(N_X+11N_W/2+62N_A)/360+769/180$ with a right positive value, where $N_X$, $N_W$, and $N_A$ are the numbers of scalar fields, Weyl fermions, and gauge fields coupled with gravity, respectively [5,12,31]. For example, $b_c=7.0$ for the Standard Model, $b_c=9.1$ for $SU\left(5\right)$ GUT, and $b_c=12.0$ for $SO\left(10\right)$ GUT. Here we adopt $b_c=7$.

Initial quantum gravity spectrum at $t=0$ is given by a two-point function of $\varphi$. It is a logarithmic function of distance [11,15,23]

$$\begin{array}{c}\hfill \langle \varphi (\eta ,\mathbf{x})\varphi (\eta ,\mathbf{y})\rangle =-\frac{1}{4b_c}log{(\mathbf{x}-\mathbf{y})}^2.\end{array}$$

(4)

Fourier transforming this in three-dimensional space provides a scale-invariant scalar spectrum whose amplitude is a positive-definite constant proportional to the reciprocal of $b_c$. This fact is a consequence of the gravitational field being a dimensionless field.

The dynamics at the UV limit, where the Riegert action plays a central role, is described by the BRST conformal field theory, i.e., a special conformal field theory with conformal invariance as gauge symmetry8. This is a hidden quantum field theory that emerges when rewriting quantum gravity into the theory on a particular background, here ${\eta }_{\mu \nu }$, and performing the path integral over $\varphi$. Unlike normal conformal invariance, this invariance shows background freedom that expresses gauge equivalence between different conformally-flat spacetimes. Under this symmetry, all ghost modes involved in the $\varphi$ and $h_{\nu }^{\mu }$ fields become unphysical—not gauge invariant. Physical states are given by scalar composite states (primary scalar states) only, no tensor-type states. Thus, the asymptotic background freedom indicates that scalar fluctuations predominate in the early universe, consistent with observations9. The argument based on the conformal algebra can be found everywhere, which corresponds to solving the Hamiltonian and momentum constraints that bind ghosts [11,12,13,14,15,16].

As mentioned before, the unphysical ghost mode is a necessary ingredient for constructing a state in which the total Hamiltonian vanishes. If all modes were positive-definite, only such a state would be the trivial vacuum. Thus, the ghost mode in the gravitational field should be a hidden entity that does not appear locally, but necessary to preserve diffeomorphism invariance and to give entropy of the universe10. This means that the quantum theory of gravity is not a theory for dealing with local entities such as gravitons, but a theory for describing state-changes in spacetime. Ghost modes do not appear directly in calculations that describe the evolution of quantum spacetime. Although they definitely contribute to the evolution, they are not directly observed. This is true not only for quantum gravity but also for classical Einstein gravity. Evolution in the Friedmann universe is exactly due to the ghost mode, and the inflationary universe is as well.

Infrared dynamics of the traceless tensor field are represented by the running coupling constant. Defining the beta function as $\mu dt/d\mu$ = $-{\beta }_0$$t^3$$({\beta }_0>0)$ with an arbitrary mass scale $\mu$ introduced upon quantization, it is given by11

$$\begin{array}{c}\hfill {\overline{t}}^2\left(Q\right)={\left[{\beta }_{0}log(Q^2/{\Lambda }_{\mathrm{QG}}^2)\right]}^{-1},\end{array}$$

(5)

where $Q^2=q^2/e^{2\varphi }$ and $q^2$ is the momentum squared defined on flat spacetime [6,22,37]. In cosmology, Q is called the physical momentum and q is called the comoving momentum. The energy scale ${\Lambda }_{\mathrm{QG}}(=\mu e^{-1/2{\beta }_0{t}^2})$ is a renormalization group (RG) invariant satisfying $d{\Lambda }_{\mathrm{QG}}/d\mu =0$, i.e., one of the physical constants that do not change in the process of cosmic evolution [44].

The effective action $\Gamma$ including quantum corrections is given by replacing $t^2$ with ${\overline{t}}^2\left(Q\right)$. Specifically, the Weyl part of the effective action is ${\Gamma }^{\mathrm{W}}=-[1/t^2-2{\beta }_0\varphi +{\beta }_{0}log(q^2/{\mu }^2)]{\overline{C}}_{\mu \nu \lambda \sigma }^2$, where the second is one of the Wess–Zumino actions and the third is a loop correction, and then we can find that the inside of the square brackets can be summarized in the form of $1/{\overline{t}}^2\left(Q\right)$ [6,22,37]. Thus, the $\varphi$-dependence of the momentum Q comes from the Wess–Zumino action S. At higher orders of t, multipoint interactions of the type ${\varphi }^n{\overline{C}}_{\mu \nu \lambda \sigma }^2$ arise, which are incorporated into Q in higher-order corrections of the running coupling constant [37]. The Riegert part of the effective action, ${\Gamma }^{\mathrm{R}}$, is also given by replacing the $t^2$-dependence in the coefficient (3) with ${\overline{t}}^2\left(Q\right)$. In this case, the Wess–Zumino action of the type ${\varphi }^{n+1}(2{\overline{∆}}_4\varphi +{\overline{E}}_4)$ contributes. Similarly, in the presence of gauge field $F_{\mu \nu }$, the Wess–Zumino action of the type ${\varphi }^n{\overline{F}}_{\mu \nu }^2$ arises and the gauge coupling constant is then replaced with its running coupling constant [37].

The effective action is a RG invariant, and the Planck mass and the cosmological constant appearing in it are also RG invariants, namely physical constants, such as the dynamical scale ${\Lambda }_{\mathrm{QG}}$ [45]. These three constants must be determined from observations. Here, the two mass scales $m_{\mathrm{pl}}$ and ${\Lambda }_{\mathrm{QG}}$ are involved in the inflationary dynamics, while the cosmological constant is negligibly small. The magnitude relation must be $m_{\mathrm{pl}}>{\Lambda }_{\mathrm{QG}}$ for inflation to occur. This relation implies that quantum gravity activates before reaching the Planck scale, i.e., the correlation length ${\xi }_{\Lambda }$ becomes larger than the Planck length. The length ${\xi }_{\Lambda }$ gives the size of localized quantum gravity excitations, leading to a dynamical discretization of spacetime [46].

The relation $m_{\mathrm{pl}}>{\Lambda }_{\mathrm{QG}}$ also serves as a unitarity condition to ensure that ghost particles do not appear in the world after the spacetime phase transition. In quantum spacetime, all ghost modes are constrained by BRST conformal invariance, but this cannot be applied after the transition. Therefore, the transition must occur below the Planck scale so that ghost gravitons with mass about $m_{\mathrm{pl}}$ do not appear.

The energy-momentum tensor is a formula obtained by applying a variation on the gravitational field to the effective action. Therefore, it is also a RG invariant and a finite operator, namely a normal product [38,39,40,41,42]. The equation of motion for quantum gravity is expressed as the vanishing energy-momentum tensor. This is a general conclusion drawn from diffeomorphism invariance and renormalizability, implying no correction due to zero-point energy, which is consistent with no UV cutoff [47]. In addition, in the theory of gravity without absolute time, conservation of the energy-momentum tensor is guaranteed by the fact that it is a RG invariant.

Until now, we have paid no attention to ℏ, but here we clearly specify where it appears when restored. Considering the action I by dividing it into the fourth-derivative gravitational part $I^{\left(4\right)}$ and the lower-derivative term $I_{\mathrm{EG}}$ consisting of the Einstein–Hilbert action and matter actions, the reciprocal of ℏ is entered only in front of $I_{\mathrm{EG}}$. Since the gravitational field is essentially a dimensionless field, $I^{\left(4\right)}$ is completely dimensionless, thus it does not contain ℏ. Similarly, ℏ is not included in the Wess–Zumino action S derived from the path integral measure. This means that all fourth-derivative gravitational actions are quantities that describe purely quantum dynamics, and the whole including them as weights, except $I_{\mathrm{EG}}$, could be regarded as a whole measure of the path integral. Hence, the path integral of quantum gravity is often symbolically expressed using only $I_{\mathrm{EG}}$ as $\int \mathcal{D}{g}_{\mu \nu }{e}^{i{I}_{\mathrm{EG}}/\hslash }$, but here it is shown that the measure is expressed exactly as $\mathcal{D}{g}_{\mu \nu }={\left[d\varphi dh\right]}_{\eta }{e}^{iS+i{I}^{\left(4\right)}}$. This is also the reason why the Weyl action and the Wess–Zumino action resulting from the measure are treated on the same footing in the following discussion. From the consideration of ℏ, it turns out that the ghost mode is essentially a quantum entity and does not appear as a classical one like particles. In the following, $\hslash =1$.

3. Inflationary Universe and Evolution Equation of Spacetime Fluctuations

The idea of inflation has been conceived to solve the horizon problem of why there were larger correlations than the horizon size in the early universe and the flatness problem of why curvature remains near zero even after more than 10 billion years [1,2,3]. Here we focus on another important role of inflation: generating initial conditions for evolution of the current universe, namely, the primordial spectra. They have to be almost scale-invariant and the magnitude of their amplitude has to be very small, at least for fluctuations with sizes involved in the structure formation of the universe. Inflation theory needs to describe these matters.

First, we show that the quantum gravity theory has an inflationary solution. The equation of motion that expresses evolution of the universe is derived by a variation of the effective action as follows [22]

$$\begin{array}{c}\hfill \delta \Gamma =\int d^{4}x({\mathbf{T}}_{\mu }^{\mu }\delta \varphi +\frac{1}{2}{\mathbf{T}}_{\nu }^{\mu }\delta h_{\mu }^{\nu })=0.\end{array}$$

(6)

The energy-momentum tensor is represented by the sum of the Riegert, Weyl, Einstein–Hilbert, and matter field sectors, as ${\mathbf{T}}_{\mu \nu }={\mathbf{T}}_{\mu \nu }^{\mathrm{R}}+{\mathbf{T}}_{\mu \nu }^{\mathrm{W}}+{\mathbf{T}}_{\mu \nu }^{\mathrm{EH}}+{\mathbf{T}}_{\mu \nu }^{\mathrm{M}}$, where ${\mathbf{T}}_{\mu \nu }={\eta }_{\mu \lambda }{\mathbf{T}}_{\nu }^{\lambda }$. Expressions for the energy-momentum tensors of each sector are summarized in Appendix A.

3.1. Inflationary Solution

The conformal-factor field $\varphi$ is divided into a spatially homogeneous component $\widehat{\varphi }$ and a fluctuation $\phi$ around it as

$$\begin{array}{c}\hfill \varphi (\eta ,\mathbf{x})=\widehat{\varphi }\left(\eta \right)+\phi (\eta ,\mathbf{x}).\end{array}$$

(7)

Homogeneous nonlinear equation for $\widehat{\varphi }$ is given from the trace equation ${\mathbf{T}}_{\mu }^{\mu }=0$ as [21,22]

$$\begin{array}{c}\hfill -\frac{b_c}{4{\pi }^2}B{\partial }_{\eta }^4\widehat{\varphi }+6M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)=0,\end{array}$$

(8)

where $M_{\mathrm{P}}=1/\sqrt{8\pi G}$ is the reduced Planck mass. The first term is from the Riegert action (2) and the second term is from the Einstein–Hilbert action, while there is no contribution from the Weyl and matter actions. In addition, letting $\rho$ be an energy density of matter fields, we obtain an energy conservation equation from the time-time component equation ${\mathbf{T}}_{00}=0$ as [22]

$$\begin{array}{c}\hfill \frac{b_c}{8{\pi }^2}B\left(2{\partial }_{\eta }^3\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }-{\partial }_{\eta }^2\widehat{\varphi }{\partial }_{\eta }^2\widehat{\varphi }\right)-3M_{\mathrm{P}}^2{\mathrm{e}}^{2\varphi }{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }+{\mathrm{e}}^{4\widehat{\varphi }}\rho =0.\end{array}$$

(9)

Here we introduce the scale factor $a=e^{\widehat{\varphi }}$ and the proper time defined by $d\tau =ad\eta$. Further, introducing the Hubble variable $H={\partial }_{\tau }a/a$ and their derivatives $Y={\partial }_{\tau }H$ and $Z={\partial }_{\tau }^{2}H$, we can rewrite (8) as

$$\begin{array}{c}\hfill B\left({\partial }_{\tau }Z+7HZ+4Y^2+18H^{2}Y+6H^4\right)-3H_{\mathbf{D}}^2\left(Y+2H^2\right)=0.\end{array}$$

(10)

When the coupling constant vanishes and thus $B=1$, this equation has a stable inflationary solution with

$$\begin{array}{c}\hfill H=H_{\mathbf{D}},H_{\mathbf{D}}=\sqrt{\frac{8{\pi }^2}{b_c}}{M}_{\mathrm{P}}.\end{array}$$

(11)

The constant $H_{\mathbf{D}}$ has a value between the Planck mass $m_{\mathrm{pl}}$ and the reduced Planck mass $M_{\mathrm{P}}$ for the typical value of $b_c$ mentioned before. Therefore, it is also called the Planck scale. This solution expresses that the scale factor expands exponentially in the proper time as $a\left(\tau \right)\propto e^{H_{\mathbf{D}}\tau }$. Thus, the universe begins to grow rapidly from the Planck time ${\tau }_{\mathrm{P}}=1/H_{\mathbf{D}}$.

The energy conservation Equation (9) can be rewritten as

$$\begin{array}{c}\hfill B\left(2HZ-Y^2+6H^{2}Y+3H^4\right)-3H_{\mathbf{D}}^2{H}^2+\frac{8{\pi }^2}{b_c}\rho =0.\end{array}$$

(12)

Substituting the solution (11) into this equation, we find that $\rho =0$ is an initial value.

3.2. Spacetime Phase Transition

As the space begins to expand, the running coupling constant $\overline{t}$ (5) increases accordingly. It implies that the universe gradually deviates from conformally-flat spacetime, and thus inflation terminates at the dynamical energy scale ${\Lambda }_{\mathrm{QG}}$ where ${\overline{t}}^2$ diverges. At this time, the conformal gravity dynamics disappears completely and the classical spacetime phase emerges. In order to represent it as a time evolution, we approximate the running coupling constant by a time-dependent mean field [22,23]. That is, replace Q with the reciprocal of the proper time $\tau$ and express it as

$$\begin{array}{c}\hfill {\overline{t}}^2\left(\tau \right)={\left[{\beta }_{0}log({\tau }_{\Lambda }^2/{\tau }^2)\right]}^{-1},\end{array}$$

(13)

where ${\tau }_{\Lambda }=1/{\Lambda }_{\mathrm{QG}}$ is called the dynamical time. This running coupling constant increases logarithmically as inflation progresses, and diverges rapidly in the dynamical time. Then, the Weyl part of the effective action disappears in proportion to the reciprocal of ${\overline{t}}^2\left(\tau \right)$.

The effective action that describes infrared dynamics of the Riegert sector is also expressed by replacing $t^2$ in the factor $B$(3) with ${\overline{t}}^2\left(\tau \right)$. We denote the factor that has undergone such a replacement as $\overline{B}\left(\tau \right)$ and consider the equation where B in the equation of motion is simply replaced with $\overline{B}\left(\tau \right)$. Furthermore, in order to express the disappearance of the dynamics of the Riegert sector with divergence of ${\overline{t}}^2\left(\tau \right)$, we assume the following form, called the dynamical factor [22,23]:

$$\begin{array}{c}\hfill \overline{B}\left(\tau \right)={\left[1+\frac{{\gamma }_1}{\kappa }\frac{{\overline{t}}^2\left(\tau \right)}{4\pi }\right]}^{-\kappa },\end{array}$$

(14)

where ${\gamma }_1,\kappa >0$. The constants ${\gamma }_1$, $\kappa$, and also ${\beta }_0$ are here treated as phenomenological parameters that adjust the inflationary scenario. In this paper, these dynamical-factor parameters are set to be ${\beta }_0=0.4$, ${\gamma }_1=0.1$, and $\kappa =1$.

In this way, the spacetime phase transition is expressed as a process in which all the fourth-derivative gravitational actions responsible for the conformal gravity dynamics vanish at the dynamical time ${\tau }_{\Lambda }$. The gravitational energy disappeared at that time is transferred to the energy of matter fields, causing the big bang. The energy conservation Equation (12) with B replaced by $\overline{B}$ expresses such a process, that is, the matter energy $\rho$, which is initially zero, is generated when $\overline{B}$ disappears as the running coupling constant increases. The change occurs rapidly near the phase transition, and $\rho \left({\tau }_{\Lambda }\right)=3M_{\mathrm{P}}^2{H}^2\left({\tau }_{\Lambda }\right)$ is obtained.

The number of e-foldings of inflation defined by ${\mathcal{N}}_e=log[a\left({\tau }_{\Lambda }\right)/a\left({\tau }_{\mathrm{P}}\right)]$ is roughly given by a ratio of the two physical scales $H_{\mathbf{D}}$ and ${\Lambda }_{\mathrm{QG}}$. Here, we take the ratio to be

$$\begin{array}{c}\hfill N=\frac{H_{\mathbf{D}}}{{\Lambda }_{\mathrm{QG}}}=60.\end{array}$$

(15)

In practice, the number of e-foldings ${\mathcal{N}}_e$ is slightly larger than this value, depending on the dynamical-factor parameters.

Since $b_c$ is set to 7 here, the dynamical energy scale is determined from (15) to be

$$\begin{array}{c}\hfill {\Lambda }_{\mathrm{QG}}=1.3\times {10}^{17}\mathrm{GeV}\end{array}$$

for $M_{\mathrm{P}}=2.4\times {10}^{18}$ GeV.

To clarify the physical meaning of the two mass scales, we introduce the corresponding comoving scales:

$$\begin{array}{c}\hfill \lambda =a\left({\tau }_i\right){\Lambda }_{\mathrm{QG}},m=a\left({\tau }_i\right)H_{\mathbf{D}},\end{array}$$

(16)

where $a\left({\tau }_i\right)$ is the scale factor at a time ${\tau }_i=1/E_i(E_i\gg H_{\mathbf{D}})$ far before the Planck time, that is, before the scale factor starts to grow. The current scale factor is normalized to 1 here. Letting $\lambda$ be a scale that explains the sharp falloff in low multipole components of the CMB angular power spectrum, its value is close to the Hubble constant $H_0$ and thus $a\left({\tau }_i\right)\sim {10}^{-59}$ is derived. This implies that most of the universe we see today was originally created within the correlation length ${\xi }_{\Lambda }$.

Furthermore, if the big bang occurred near ${10}^{17}$ GeV and the present universe was created there, then the universe expanded about ${10}^{29}$ times after the big bang. From this, it is required to expand ${10}^{30}$ times during the inflationary period. Although ${\Lambda }_{\mathrm{QG}}$ is not the only parameter determining the number of e-foldings, it turns out that this scale is a relatively tight one when considering the scenario of cosmic evolution. In this paper, we adopt $m=0.02{\mathrm{Mpc}}^{-1}$ as a reference value, and then $\lambda (=m/N)=0.00033{\mathrm{Mpc}}^{-1}$12.

A numerical solution of the equation of motion for the homogeneous component is shown in Figure 313. The horizontal axis is a dimensionless time defined by $t=H_{\mathbf{D}}\tau$ (not to be confused with the coupling constant). The Hubble variable H is obtained by solving the homogeneous Equation (10) with B replaced by $\overline{B}$. The matter energy density $\rho$ is then obtained by substituting the result of H into (12). See Appendix B for details of the computation method.

We will now consider how the low-energy spacetime after the phase transition is described. In QCD, when the gluon action given by the second derivative disappears, the gluon dynamics disappear completely. Meanwhile, in quantum gravity, even if the fourth-derivative gravitational action disappears, the second-derivative Einstein–Hilbert action remains, which gives the kinetic term of the gravitational field at low energies. Thus, Einstein’s theory of gravity appears as a low-energy effective theory. It is then no longer appropriate to consider separately the conformal-factor field and the traceless tensor field, and a composite field in which they are tightly coupled becomes a new variable of the gravitational field. The actual low-energy effective theory will be described by derivative expansions in the gravitational field with the Einstein–Hilbert action as the lowest order. According to this, the Hubble variable increases a little further after the phase transition, gains the e-foldings, and then oscillates and approaches the Friedmann solution rapidly, as shown in Figure 3 (see Ref. [22] for details).

3.3. Nonlinear Evolution Equation of Fluctuations

Here, we consider equations of motion for fluctuations around the inflationary solution. In fact, we should discuss the evolution of the physical scalars defined by composite fields, such as the quantum Ricci scalar. However, due to the fact that the scalar fluctuations reduce in amplitude during inflation, as previously mentioned, it turns out that it is sufficient to consider the evolution of the fluctuations, because the information we ultimately want to know is their reduced values at the phase transition, which gives the initial conditions of the Friedmann universe.

Now we derive a system of nonlinear evolution equations for the fluctuations. The conformal invariance indicates that the scalar fluctuation $\phi$ (7) dominates in the early universe beyond the Planck scale. In order to describe the process that the scalar fluctuation gradually deforms with time, we also take into account a scalar component of the traceless tensor field, defined by

$$\begin{array}{c}\hfill h_{00}=h,h_{ij}=\frac{1}{3}{\delta }_{ij}h,\end{array}$$

(17)

which is initially vanishing.

The system of linear evolution equations for the fluctuations has already been derived, and it has been shown that the fluctuations gradually reduce during inflation [22,23]. However, since $\phi$ is still large in the initial stage, we cannot accurately estimate the fluctuation spectrum without incorporating nonlinear terms involving multipoint functions of $\phi$ there. The exponential factor of $\phi$ existing in the Einstein–Hilbert action provides such terms. On the other hand, h is treated in the linear approximation, assuming that it is not very large, which is because h is initially small. It then grows, but reduces again before the coupling constant t becomes large. Thus, nonlinear terms are limited to first order in h.

In addition, scalar fluctuations contain nonlinear terms combining other tensor modes into scalars. These modes may also contribute to the equation of motion as nonlinear terms. However, since the initial smallness of these tensor modes is almost maintained [22], the nonlinear terms containing them will contribute little to the evolution of the scalar fluctuations.

First, consider the trace part of the equations of motion. The Riegert sector is given by ${\mathbf{T}}_{\mu }^{\mathrm{R}\mu }=-(b_c/16{\pi }^2)B(4{\overline{∆}}_4\varphi +{\overline{E}}_4)$. Expanding it to the first order of the traceless tensor field yields (A2). We expand this expression according to (7) and (17) and decompose it into the linear part of $\phi$ and h and the nonlinear part of $o\left(\phi h\right)$. Writing it as ${\mathbf{T}}_{\mu }^{\mathrm{R}\mu }={\mathbf{T}}_{\mu }^{\mathrm{R}\mu }{{|}_L+{\mathbf{T}}_{\mu }^{\mathrm{R}\mu }|}_{NL}$, we obtain

$$\begin{array}{ccc}\hfill {\mathbf{T}}_{\mu }^{\mathrm{R}\mu }{|}_L& =& \frac{b_c}{8{\pi }^2}B\{-2{\partial }_{\eta }^4\phi +4{\partial }_{\eta }^2|{\partial }^2\phi -2|{\partial }^4\phi -4{\partial }_{\eta }^4\widehat{\varphi }h-8{\partial }_{\eta }^3\widehat{\varphi }{\partial }_{\eta }h-\frac{16}{3}{\partial }_{\eta }^2\widehat{\varphi }{\partial }_{\eta }^{2}h\hfill \\ & & +\frac{32}{9}{\partial }_{\eta }^2\widehat{\varphi }|{\partial }^{2}h-\frac{4}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }^{3}h+\frac{20}{9}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }|{\partial }^{2}h-\frac{1}{3}{\partial }_{\eta }^{4}h+\frac{2}{9}{\partial }_{\eta }^2|{\partial }^{2}h+\frac{1}{9}|{\partial }^{4}h\}\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\mathbf{T}}_{\mu }^{\mathrm{R}\mu }{|}_{NL}& =& \frac{b_c}{8{\pi }^2}B\{-4h{\partial }_{\eta }^4\phi +\frac{8}{3}h{\partial }_{\eta }^2|{\partial }^2\phi +\frac{4}{3}h|{\partial }^4\phi -8{\partial }_{\eta }h{\partial }_{\eta }^3\phi +\frac{8}{3}{\partial }^{i}h{\partial }_{\eta }^2{\partial }_i\phi \hfill \\ & & +\frac{8}{3}{\partial }_{\eta }h{\partial }_{\eta }|{\partial }^2\phi +\frac{8}{3}{\partial }^{i}h{\partial }_i|{\partial }^2\phi -\frac{16}{3}{\partial }_{\eta }^{2}h{\partial }_{\eta }^2\phi +\frac{32}{9}|{\partial }^{2}h{\partial }_{\eta }^2\phi +\frac{16}{9}|{\partial }^{2}h|{\partial }^2\phi \hfill \\ & & -\frac{4}{3}{\partial }_{\eta }^{3}h{\partial }_{\eta }\phi +\frac{20}{9}{\partial }_{\eta }|{\partial }^{2}h{\partial }_{\eta }\phi -\frac{4}{3}{\partial }_{\eta }^2{\partial }^{i}h{\partial }_i\phi +\frac{4}{9}{\partial }^i|{\partial }^{2}h{\partial }_i\phi \}.\hfill \end{array}$$

On the other hand, the trace part of the Weyl sector disappears in the first order of h, resulting in ${\mathbf{T}}_{\lambda }^{\mathrm{W}\lambda }{|}_L=0$.

The trace part of the Einstein–Hilbert sector is given by ${\mathbf{T}}_{\mu }^{\mathrm{EH}\mu }=M_{\mathrm{P}}^2\sqrt{-g}R=M_{\mathrm{P}}^2{e}^{2\varphi }(\overline{R}-6{\overline{\nabla }}^2\varphi -6{\overline{\nabla }}^{\mu }\varphi {\overline{\nabla }}_{\mu }\varphi )$. Expanding it to the first order of the traceless tensor field, we get (A5). Expanding this further yields terms up to infinite order in $\phi$ due to the presence of the exponential factor $e^{2\varphi }$. Expressing it as ${\mathbf{T}}_{\mu \nu }^{\mathrm{EH}}={\sum }_{n=0}^{\infty }{\mathbf{T}}_{\left(n\right)\mu \nu }^{\mathrm{EH}}$, $n=0$ is a homogeneous component that gives the second term of (8) and (9), and ${\mathbf{T}}_{\left(n\right)\mu \nu }^{\mathrm{EH}}$ consists of the $o\left({\phi }^n\right)$ and $o\left({\phi }^{n-1}h\right)$ terms. The linear part is given by

$$\begin{array}{ccc}\hfill {\mathbf{T}}_{\left(1\right)\mu }^{\mathrm{EH}\mu }& =& M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\{6{\partial }_{\eta }^2\phi -6|{\partial }^2\phi +12{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\phi +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)\phi \hfill \\ & & +{\partial }_{\eta }^{2}h+\frac{1}{3}|{\partial }^{2}h+6{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h+6\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h\}\hfill \end{array}$$

and the nonlinear part composed of $o\left({\phi }^2\right)$ and $o\left(\phi h\right)$ is

$$\begin{array}{ccc}\hfill {\mathbf{T}}_{\left(2\right)\mu }^{\mathrm{EH}\mu }& =& M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\{12\phi {\partial }_{\eta }^2\phi -12\phi |{\partial }^2\phi +6{\partial }_{\eta }\phi {\partial }_{\eta }\phi -6{\partial }^i\phi {\partial }_i\phi +24{\partial }_{\eta }\widehat{\varphi }\phi {\partial }_{\eta }\phi \hfill \\ & & +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right){\phi }^2+2\phi {\partial }_{\eta }^{2}h+\frac{2}{3}\phi |{\partial }^{2}h+6{\partial }_{\eta }^2\phi h+2|{\partial }^2\phi h\hfill \\ & & +6{\partial }_{\eta }\phi {\partial }_{\eta }h+2{\partial }^i\phi {\partial }_{i}h+12{\partial }_{\eta }\widehat{\varphi }\phi {\partial }_{\eta }h+12{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\phi h\hfill \\ & & +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)\phi h\}.\hfill \end{array}$$

The case of $n\ge 3$ can also be easily derived.

Since there are two field variables, the trace equation alone is incomplete. Therefore, we consider another equation that combines the energy-momentum tensor such as

$$\begin{array}{c}\hfill \frac{1}{|{\partial }^2}\left({\mathbf{T}}_i^i-3\frac{{\partial }^i{\partial }^j}{|{\partial }^2}{\mathbf{T}}_{ij}\right)=0.\end{array}$$

(18)

Unlike the trace equation, there is a contribution from the Weyl sector, which is derived from the linear part of the energy-momentum tensor given in (A3). On the other hand, there is no contribution from the matter sector in this case either.

Since there is $1/t^2$ in front of the Weyl part, this part dominates at the early stage when $t^2$ is small. Multiplying the whole equation of motion by $t^2$ to remove it, we can then see that the contribution from the Riegert and Einstein–Hilbert sectors becomes $o\left(t^2\right)$ so that they can be treated in linear in the early stage. The contribution to (18) from the Riegert sector is thus given by linear part of $\phi$ and h obtained by substituting (7) into (A1). The contribution from the Einstein–Hilbert sector is also derived from (A4).

Furthermore, since the fluctuations decrease as the coupling constant t increases, Equation (18) can remain the linear approximation from beginning to end. This equation is expressed by the second derivative of the fields as

$$\begin{array}{ccc}& & \frac{16}{3}\frac{1}{t^2}\left({\partial }_{\eta }^{2}h-\frac{1}{3}|{\partial }^{2}h\right)+\frac{b_c}{8{\pi }^2}B\{\frac{4}{3}{\partial }_{\eta }^2\phi +\frac{8}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\phi +\left(8{\partial }_{\eta }^2\widehat{\varphi }-\frac{4}{3}|{\partial }^2\right)\phi \hfill \\ & & +\frac{2}{9}{\partial }_{\eta }^{2}h+\frac{4}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h+\left(\frac{20}{9}{\partial }_{\eta }^2\widehat{\varphi }-\frac{16}{9}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }+\frac{2}{27}|{\partial }^2\right)h\}\hfill \\ & & +M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\left(-4\phi +\frac{2}{3}h\right)=0.\hfill \end{array}$$

(19)

Thus, the system of evolution equations can be expressed as a coupled differential equation of this linear equation and the nonlinear trace equation ${\mathbf{T}}_{\mu }^{\mathrm{R}\mu }{{|}_L+{\mathbf{T}}_{\mu }^{\mathrm{R}\mu }|}_{NL}+{\sum }_{n=1}^{\infty }{\mathbf{T}}_{\left(n\right)\mu }^{\mathrm{EH}\mu }=0$.

However, it is not yet in a form that can be solved numerically. It is necessary to rewrite the ${\partial }_{\eta }^4\widehat{\varphi }h$ term in the ${\mathbf{T}}_{\mu }^{\mathrm{R}\mu }{|}_L$ part to terms with lower time-derivatives of $\widehat{\varphi }$ using the homogeneous Equation (8). In addition, noting that the trace equation has the structure

$$\begin{array}{ccc}& & -\frac{b_c}{4{\pi }^2}B\left({\partial }_{\eta }^4\phi -2{\partial }_{\eta }^2|{\partial }^2\phi +|{\partial }^4\phi \right)+M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\{6{\partial }_{\eta }^2\phi -6|{\partial }^2\phi \hfill \\ & & +12{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\phi +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)\phi \}=o\left(h\right),\hfill \end{array}$$

the $h{\partial }_{\eta }^4\phi$ term in the ${\mathbf{T}}_{\mu }^{\mathrm{R}\mu }{|}_{NL}$ part is rewritten to terms with lower time-derivatives of $\phi$ by ignoring the $o\left(h^2\right)$ term. In addition, we introduce a field variable

$$\begin{array}{c}\hfill \Phi =\phi +\frac{1}{6}h\end{array}$$

and treat $\Phi$ and h as independent variables. The fourth time-derivative of h in the linear part is then absorbed into a fourth time-derivative of $\Phi$ so that the time-derivative of h is given by up to at most the third one. Here, $\Phi$ is the gravitational potential when the line-element is expressed within linear as $d{s}^2=a^2[-(1+2\Psi )d{\eta }^2+(1+2\Phi )d{\mathbf{x}}^2]$, where the other one is given by $\Psi =\phi -h/2$.

The nonlinear terms that become effective in the early stage of the inflationary universe can be incorporated in this way. On the other hand, the nonlinear effect arising in infrared conformal gravity dynamics is represented by the running coupling constant. As mentioned earlier, it is expressed by adopting a type of the mean-field approximation that rewrites the coupling constant squared to the time-dependent average ${\overline{t}}^2\left(\tau \right)$(13) and also B to the time-dependent dynamical factor $\overline{B}\left(\tau \right)$(14). While this seems like a radical model, no detailed information of strong coupling dynamics other than the disappearance of conformal gravity dynamics will be needed, because the size of the fluctuations considered here becomes much larger than the correlation length ${\xi }_{\Lambda }$ near the phase transition point where the coupling constant increases.

The system of nonlinear evolution equations derived finally is given by

$$\begin{array}{ccc}& & \frac{b_c}{8{\pi }^2}\overline{B}\left(\tau \right)\{-2{\partial }_{\eta }^4\Phi +4{\partial }_{\eta }^2|{\partial }^2\Phi -2|{\partial }^4\Phi -\frac{4}{9}{\partial }_{\eta }^2|{\partial }^{2}h+\frac{4}{9}|{\partial }^{4}h-8{\partial }_{\eta }^3\widehat{\varphi }{\partial }_{\eta }h\hfill \\ & & -\frac{16}{3}{\partial }_{\eta }^2\widehat{\varphi }{\partial }_{\eta }^{2}h-\frac{4}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }^{3}h+\frac{32}{9}{\partial }_{\eta }^2\widehat{\varphi }|{\partial }^{2}h+\frac{20}{9}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }|{\partial }^{2}h\hfill \\ & & -\frac{16}{3}h{\partial }_{\eta }^2|{\partial }^2\Phi +\frac{16}{3}h|{\partial }^4\Phi -8{\partial }_{\eta }h{\partial }_{\eta }^3\Phi +\frac{8}{3}{\partial }^{i}h{\partial }_{\eta }^2{\partial }_i\Phi +\frac{8}{3}{\partial }_{\eta }h{\partial }_{\eta }|{\partial }^2\Phi \hfill \\ & & +\frac{8}{3}{\partial }^{i}h{\partial }_i|{\partial }^2\Phi -\frac{16}{3}{\partial }_{\eta }^{2}h{\partial }_{\eta }^2\Phi +\frac{32}{9}|{\partial }^{2}h{\partial }_{\eta }^2\Phi +\frac{16}{9}|{\partial }^{2}h|{\partial }^2\Phi -\frac{4}{3}{\partial }_{\eta }^{3}h{\partial }_{\eta }\Phi \hfill \\ & & +\frac{20}{9}{\partial }_{\eta }|{\partial }^{2}h{\partial }_{\eta }\Phi -\frac{4}{3}{\partial }_{\eta }^2{\partial }^{i}h{\partial }_i\Phi +\frac{4}{9}{\partial }^i|{\partial }^{2}h{\partial }_i\Phi \}\hfill \\ & & +M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\{6{\partial }_{\eta }^2\Phi +12{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\Phi -6|{\partial }^2\Phi +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)\Phi \hfill \\ & & +4{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h+\frac{4}{3}|{\partial }^{2}h-8\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h\hfill \\ & & +12\Phi {\partial }_{\eta }^2\Phi +6{\partial }_{\eta }\Phi {\partial }_{\eta }\Phi +24{\partial }_{\eta }\widehat{\varphi }\Phi {\partial }_{\eta }\Phi +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right){\Phi }^2\hfill \\ & & -12\Phi |{\partial }^2\Phi -6{\partial }^i\Phi {\partial }_i\Phi -8h{\partial }_{\eta }^2\Phi +4{\partial }_{\eta }h{\partial }_{\eta }\Phi -16{\partial }_{\eta }\widehat{\varphi }h{\partial }_{\eta }\Phi +8{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h\Phi \hfill \\ & & -16\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h\Phi +16h|{\partial }^2\Phi +4{\partial }^{i}h{\partial }_i\Phi +\frac{8}{3}|{\partial }^{2}h\Phi \}=0\hfill \end{array}$$

(20)

and

$$\begin{array}{ccc}& & \frac{16}{3}\frac{1}{{\overline{t}}^2\left(\tau \right)}\left({\partial }_{\eta }^{2}h-\frac{1}{3}|{\partial }^{2}h\right)+\frac{b_c}{8{\pi }^2}\overline{B}\left(\tau \right)\{\frac{4}{3}{\partial }_{\eta }^2\Phi -\frac{4}{3}|{\partial }^2\Phi +\frac{8}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\Phi +8{\partial }_{\eta }^2\widehat{\varphi }\Phi \hfill \\ & & +\frac{8}{27}|{\partial }^{2}h+\frac{8}{9}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h+\left(\frac{8}{9}{\partial }_{\eta }^2\widehat{\varphi }-\frac{16}{9}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h\}+M_{\mathrm{P}}{e}^{2\widehat{\varphi }}\left(-4\Phi +\frac{4}{3}h\right)\hfill \\ & & =0.\hfill \end{array}$$

(21)

The second Equation (21) plays an important role as a constraint for connecting the inflation and Friedmann phases. When the running coupling constant $\overline{t}$ is small, h is small, indicating that $\phi$ is dominant. Conversely, when $\overline{t}$ diverges, $1/{\overline{t}}^2$ the dynamical factor $\overline{B}$ vanishes, and the last Einstein term becomes dominant, this results in a transition to the Friedmann universe satisfying $\Phi =h/3$. In terms of the gravitational potentials, it expresses that the initial fluctuation of $\Phi =\Psi$ changes to $\Phi =-\Psi$.

Equations (20) and (21) will be solved simultaneously. At that time, it is necessary that ${\partial }_{\eta }^{3}h$ in (20) is rewritten to the third time-derivative of $\Phi$ and other lower time-derivative terms using an equation obtained by time differentiating (21) as

$$\begin{array}{ccc}\hfill {\partial }_{\eta }^{3}h& =& \frac{1}{3}{\partial }_{\eta }|{\partial }^{2}h+\frac{1}{{\overline{t}}^2\left(\tau \right)}{\partial }_{\eta }{\overline{t}}^2\left(\tau \right)\left({\partial }_{\eta }^{2}h-\frac{1}{3}|{\partial }^{2}h\right)\hfill \\ & & -\frac{b_c}{32{\pi }^2}{\overline{t}}^2\left(\tau \right)\left\{\overline{B}\left(\tau \right)\right[{\partial }_{\eta }^3\Phi +2{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }^2\Phi +8{\partial }_{\eta }^2\widehat{\varphi }{\partial }_{\eta }\Phi +6{\partial }_{\eta }^3\widehat{\varphi }\Phi -{\partial }_{\eta }|{\partial }^2\Phi \hfill \\ & & +\frac{2}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }^{2}h+\frac{4}{3}\left({\partial }_{\eta }^2\widehat{\varphi }-{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right){\partial }_{\eta }h+\frac{2}{3}\left({\partial }_{\eta }^3\widehat{\varphi }-4{\partial }_{\eta }^2\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h+\frac{2}{9}{\partial }_{\eta }|{\partial }^{2}h]\hfill \\ & & +{\partial }_{\eta }\overline{B}\left(\tau \right)[{\partial }_{\eta }^2\Phi +2{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\Phi +6{\partial }_{\eta }^2\widehat{\varphi }\Phi -|{\partial }^2\Phi +\frac{2}{3}{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h\hfill \\ & & +\frac{2}{3}\left({\partial }_{\eta }^2\widehat{\varphi }-2{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h+\frac{2}{9}|{\partial }^{2}h\left]\right\}\hfill \\ & & +M_{\mathrm{P}}^2{\overline{t}}^2\left(\tau \right)e^{2\widehat{\varphi }}\left\{\frac{3}{4}{\partial }_{\eta }\Phi +\frac{3}{2}{\partial }_{\eta }\widehat{\varphi }\Phi -\frac{1}{4}{\partial }_{\eta }h-\frac{1}{2}{\partial }_{\eta }\widehat{\varphi }h\right\}.\hfill \end{array}$$

(22)

4. How to Handle Nonlinear Terms

The fluctuation spectrum is expressed using Fourier transform in the three-dimensional comoving coordinate space. Letting $f\left(\mathbf{x}\right)$ be a dimensionless field variable, here we define its Fourier transform so that the transformed quantity is also dimensionless as

$$\begin{array}{c}\hfill f\left(\mathbf{x}\right)=\int \frac{d^3\mathbf{k}}{{\left(2\pi \right)}^3}\frac{1}{k^3}f\left(\mathbf{k}\right)e^{i\mathbf{k}·\mathbf{x}},\end{array}$$

where $k=\left|\mathbf{k}\right|$. Its dimensionless power spectrum ${\langle \left|f\left(\mathbf{k}\right)\right|}^2\rangle /2{\pi }^2$ is defined by $\langle f\left(\mathbf{k}\right)f\left({\mathbf{k}}^{\prime }\right)\rangle =\langle |f\left(\mathbf{k}\right){|}^2\rangle k^3{\left(2\pi \right)}^3{\delta }^3(\mathbf{k}+{\mathbf{k}}^{\prime })$. The scale invariance is expressed as this power spectrum being constant.

Fourier transform of a dimensionless field product $f\left(\mathbf{x}\right)g\left(\mathbf{x}\right)$ is also defined by

$$\begin{array}{ccc}\hfill \left(fg\right)\left(\mathbf{k}\right)& =& k^3\int d^3\mathbf{x}f\left(\mathbf{x}\right)g\left(\mathbf{x}\right)e^{-i\mathbf{k}·\mathbf{x}}\hfill \\ & =& k^3\int d^3\mathbf{x}\int \frac{d^3\mathbf{p}}{{\left(2\pi \right)}^3}\frac{d^3\mathbf{q}}{{\left(2\pi \right)}^3}\frac{1}{p^3{q}^3}f\left(\mathbf{p}\right)g\left(\mathbf{q}\right)e^{i(\mathbf{p}+\mathbf{q}-\mathbf{k})·\mathbf{x}}.\hfill \end{array}$$

(23)

Here, we rewrite the exponential part using the expansion formula with the spherical Bessel function,

$$\begin{array}{c}\hfill e^{i\mathbf{k}·\mathbf{x}}=4\pi \sum _{l=0}^{\infty }\sum _{m=-l}^l{i}^l{j}_l\left(kx\right)Y_{lm}^{\ast }\left({\Omega }_k\right)Y_{lm}\left({\Omega }_x\right),\end{array}$$

where $\Omega$ denotes angular components such that $d{\mathbf{x}}^3=x^{2}dxd{\Omega }_x$ and the spherical harmonics are normalized as $\int d\Omega Y_{lm}^{\ast }(\Omega )Y_{l^{\prime }{m}^{\prime }}(\Omega )={\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}$.

Henceforth, assuming that isotropic components of each function in the linear and nonlinear terms contribute predominantly, let

$$\begin{array}{c}\hfill f\left(\mathbf{p}\right)=f\left(p\right),g\left(\mathbf{q}\right)=g\left(q\right).\end{array}$$

(24)

In this case, the angular integration in (23) can be easily performed and we obtain

$$\begin{array}{c}\hfill \left(fg\right)\left(k\right)=\frac{8}{{\left(2\pi \right)}^3}{k}^3\int \frac{dp}{p}\frac{dq}{q}f\left(p\right)g\left(q\right)\int dx{x}^2{j}_0\left(kx\right)j_0\left(px\right)j_0\left(qx\right),\end{array}$$

where $j_0\left(kx\right)=sin\left(kx\right)/kx$. The x-integral yields

$$\begin{array}{c}\hfill {\int }_0^{\infty }dx{x}^2{j}_0\left(kx\right)j_0\left(px\right)j_0\left(qx\right)=\left\{\begin{array}{cc}\frac{\pi }{4}\frac{1}{kpq}\hfill & for|p-q|<k<p+q\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

and thus we get

$$\begin{array}{c}\hfill \left(fg\right)\left(k\right)=\frac{k^2}{{\left(2\pi \right)}^2}\int \frac{dp}{p^2}\frac{dq}{q^2}\theta (k-|p-q\left|\right)\theta (p+q-k)f\left(p\right)g\left(q\right),\end{array}$$

(25)

where $\theta$ is the Heaviside step function.

The Fourier transform of a product of functions with spatial derivatives can be calculated similarly. Here it suffices for us to consider the following type:

$$\begin{array}{ccc}& & \left({\partial }_{i}f{\partial }^{i}g\right)\left(k\right)=k^3\int d^3\mathbf{x}{\partial }_{i}f\left(\mathbf{x}\right){\partial }^{i}g\left(\mathbf{x}\right)e^{-i\mathbf{k}·\mathbf{x}}\hfill \\ & & =k^3\int d^3\mathbf{x}\int \frac{d^3\mathbf{p}}{{\left(2\pi \right)}^3}\frac{d^3\mathbf{q}}{{\left(2\pi \right)}^3}\frac{-\mathbf{p}·\mathbf{q}}{p^3{q}^3}f\left(\mathbf{p}\right)g\left(\mathbf{q}\right)e^{i(\mathbf{p}+\mathbf{q}-\mathbf{k})·\mathbf{x}}.\hfill \end{array}$$

Under the isotropic condition (24), integrating the angular components of the momentum $\mathbf{p}$ and $\mathbf{q}$ can be easily performed by expressing the vector $\mathbf{p}=(p_x,p_y,p_z)$ in terms of spherical harmonics as

$$\begin{array}{ccc}\hfill p_x& =& \sqrt{\frac{2\pi }{3}}p\left\{Y_{1-1}\left({\Omega }_p\right)-Y_{11}\left({\Omega }_p\right)\right\},\hfill \\ \hfill p_y& =& i\sqrt{\frac{2\pi }{3}}p\left\{Y_{1-1}\left({\Omega }_p\right)+Y_{11}\left({\Omega }_p\right)\right\},\hfill \\ \hfill p_z& =& \sqrt{\frac{4\pi }{3}}p{Y}_{10}\left({\Omega }_p\right).\hfill \end{array}$$

Performing the integration of the remaining angular component ${\Omega }_x$ using

$$\begin{array}{c}\hfill \int d{\Omega }_x{Y}_{l_1{m}_1}^{\ast }\left({\Omega }_x\right)Y_{l_2{m}_2}\left({\Omega }_x\right)Y_{l_3{m}_3}\left({\Omega }_x\right)=\sqrt{\frac{(2l_2+1)(2l_3+1)}{4\pi (2l_1+1)}}{C}_{l_{2}0l_{3}0}^{l_{1}0}{C}_{l_2{m}_2{l}_3{m}_3}^{l_1{m}_1}\end{array}$$

yields

$$\begin{array}{ccc}\hfill \left({\partial }_{i}f{\partial }^{i}g\right)\left(k\right)& =& \frac{8}{{\left(2\pi \right)}^3}{k}^3\int dpdqf\left(p\right)g\left(q\right)\int dx{x}^2{j}_0\left(kx\right)j_1\left(px\right)j_1\left(qx\right)\hfill \\ & =& \frac{k^2}{{\left(2\pi \right)}^2}\int \frac{dp}{p^2}\frac{dq}{q^2}\theta (k-|p-q\left|\right)\theta (p+q-k)\hfill \\ & & \times \frac{1}{2}(p^2+q^2-k^2)f\left(p\right)g\left(q\right).\hfill \end{array}$$

(26)

The lower bound of the comoving momentum integral is given by the physical comoving scale $\lambda$(16), while there is no upper bound on the integration. Using (25) and (26), Fourier transforms of the second-order nonlinear terms in the trace Equation (20) are given by

$$\begin{array}{ccc}& & \frac{b_c}{8{\pi }^2}\overline{B}\frac{k^2}{{\left(2\pi \right)}^2}{\int }_{\lambda }^{\infty }\frac{dp}{p^2}\frac{dq}{q^2}\theta (k-|p-q\left|\right)\theta (p+q-k)\{-8{\partial }_{\eta }h\left(p\right){\partial }_{\eta }^3\Phi \left(q\right)\hfill \\ & & -\frac{16}{3}{\partial }_{\eta }^{2}h\left(p\right){\partial }_{\eta }^2\Phi \left(q\right)-\frac{4}{3}{\partial }_{\eta }^{3}h\left(p\right){\partial }_{\eta }\Phi \left(q\right)+\left(-\frac{20}{9}{p}^2+\frac{20}{3}{q}^2-\frac{4}{3}{k}^2\right)h\left(p\right){\partial }_{\eta }^2\Phi \left(q\right)\hfill \\ & & +\left(-\frac{20}{9}{p}^2-\frac{8}{3}{q}^2\right){\partial }_{\eta }h\left(p\right){\partial }_{\eta }\Phi \left(q\right)+\left(-\frac{2}{3}{p}^2-\frac{2}{3}{q}^2+\frac{2}{3}{k}^2\right){\partial }_{\eta }^{2}h\left(p\right)\Phi \left(q\right)\hfill \\ & & +\left[-\frac{2}{9}{p}^4+4q^4+\frac{2}{9}{p}^2{q}^2+\left(\frac{2}{9}{p}^2+\frac{4}{3}{q}^2\right)k^2\right]h\left(p\right)\Phi \left(q\right)\}\hfill \end{array}$$

(27)

and

$$\begin{array}{ccc}& & M_{\mathrm{P}}^2{e}^{2\widehat{\varphi }}\frac{k^2}{{\left(2\pi \right)}^2}{\int }_{\lambda }^{\infty }\frac{dp}{p^2}\frac{dq}{q^2}\theta (k-|p-q\left|\right)\theta (p+q-k)\{6\Phi \left(p\right){\partial }_{\eta }^2\Phi \left(q\right)\hfill \\ & & +6{\partial }_{\eta }^2\Phi \left(p\right)\Phi \left(q\right)+6{\partial }_{\eta }\Phi \left(p\right){\partial }_{\eta }\Phi \left(q\right)+12{\partial }_{\eta }\widehat{\varphi }\Phi \left(p\right){\partial }_{\eta }\Phi \left(q\right)+12{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\Phi \left(p\right)\Phi \left(q\right)\hfill \\ & & +12\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)\Phi \left(p\right)\Phi \left(q\right)+3\left(p^2+q^2+k^2\right)\Phi \left(p\right)\Phi \left(q\right)-8h\left(p\right){\partial }_{\eta }^2\Phi \left(q\right)\hfill \\ & & +4{\partial }_{\eta }h\left(p\right){\partial }_{\eta }\Phi \left(q\right)-16{\partial }_{\eta }\widehat{\varphi }h\left(p\right){\partial }_{\eta }\Phi \left(q\right)+8{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }h\left(p\right)\Phi \left(q\right)\hfill \\ & & -16\left({\partial }_{\eta }^2\widehat{\varphi }+{\partial }_{\eta }\widehat{\varphi }{\partial }_{\eta }\widehat{\varphi }\right)h\left(p\right)\Phi \left(q\right)-\left(\frac{2}{3}{p}^2+14q^2+2k^2\right)h\left(p\right)\Phi \left(q\right)\}.\hfill \end{array}$$

(28)

Here, the momentum integral in (28) has no problem on the upper bound due to the strength of its negative power. On the other hand, the absence of problems in (27) derived from the fourth-derivative Riegert action will be guaranteed by the property that the traceless tensor field becomes smaller in the high-momentum region14.

5. Simplification of Nonlinear Terms and Numerical Evaluation

In order to numerically solve the above system of nonlinear evolution equations, we have to describe the momentum integrals by sums and rewrite it to a multiple simultaneous differential equation in which many fields with different momentum are linked. However, it is not easy to actually solve such a multiple differential equation. Here we further simplify the equations and examine them to evaluate how the nonlinear term specifically affects patterns of the spectrum. At that time, we do not care too much about the overall amplitude accuracy.

First of all, since we cannot solve the evolution equation if two time variables are mixed, we rewrite the conformal time $\eta$ to the proper time $\tau$. Using the scale factor a, the Hubble variable H, and its time derivatives $Y={\partial }_{\tau }H$ and $Z={\partial }_{\tau }^{2}H$, the differential operators with respect to $\eta$ are rewritten as

$$\begin{array}{ccc}\hfill {\partial }_{\eta }& =& a{\partial }_{\tau },\hfill \\ \hfill {\partial }_{\eta }^2& =& a^2{D}_{\tau }^2,D_{\tau }^2={\partial }_{\tau }^2+H{\partial }_{\tau },\hfill \\ \hfill {\partial }_{\eta }^3& =& a^3{D}_{\tau }^3,D_{\tau }^3={\partial }_{\tau }^3+3H{\partial }_{\tau }^2+\left(Y+2H^2\right){\partial }_{\tau },\hfill \\ \hfill {\partial }_{\eta }^4& =& a^4{D}_{\tau }^4,D_{\tau }^4={\partial }_{\tau }^4+6H{\partial }_{\tau }^3+\left(4Y+11H^2\right){\partial }_{\tau }^2+\left(Z+7HY+6H^3\right){\partial }_{\tau }.\hfill \end{array}$$

The inflationary variables are also rewritten as

$$\begin{array}{ccc}\hfill {\partial }_{\eta }\widehat{\varphi }& =& aH,{\partial }_{\eta }^2\widehat{\varphi }=a^2\left(Y+H^2\right),{\partial }_{\eta }^3\widehat{\varphi }=a^3\left(Z+4HY+2H^3\right),\hfill \\ \hfill {\partial }_{\eta }^4\widehat{\varphi }& =& a^4\left({\partial }_{\tau }Z+7HZ+4Y^2+18H^{2}Y+6H^4\right),\hfill \end{array}$$

where the last ${\partial }_{\eta }^4\widehat{\varphi }$ was used when deriving (10).

We decompose the scale factor as $a\left(t\right)=a\left(t_i\right)\times \overline{a}\left(t\right)$, where $t(=H_{\mathbf{D}}\tau )$ is the dimensionless time as introduced before and $\overline{a}\left(t_i\right)=1$. Then, a dimensionless quantity obtained by dividing the physical momentum squared $k^2/a^2$ by $H_{\mathbf{D}}^2$ can be expressed using the comoving Planck mass $m$(16) as

$$\begin{array}{c}\hfill P_a(t,k)=\frac{k^2}{H_{\mathbf{D}}^2{a}^2\left(\tau \right)}=\frac{k^2}{m^2{\overline{a}}^2\left(t\right)}.\end{array}$$

(29)

Similarly, we introduce a dimensionless differential operator defined by $D_t^n=D_{\tau }^n/H_{\mathbf{D}}^n={\partial }_t^n+\cdots (n=2,3,4)$, in which the Hubble variables are replaced with normalized $\overline{H}=H/H_{\mathbf{D}}$, $\overline{Y}=Y/H_{\mathbf{D}}^2$, and $\overline{Z}=Z/H_{\mathbf{D}}^3$.

Next, we simplify the nonlinear terms as follows. First, applying an approximation that the field correlation vanishes instantaneously when distance exceeds the correlation length ${\xi }_{\Lambda }$, let the field $f\left(k\right)$ be a variable with values such as a step function for $k\ge \lambda$. We then focus on the following function that arises in the Fourier transforms of the field products:

$$\begin{array}{c}\hfill \frac{1}{{\left(2\pi \right)}^2}\frac{k^2}{p^2{q}^2}\theta (k-|p-q\left|\right)\theta (p+q-k).\end{array}$$

(30)

When displaying this function in the physical range $p,q\ge \lambda$, there are two sharp peaks at $(p,q)=(\lambda ,k)$ and $(k,\lambda )$. Therefore, we extract these neighborhoods and evaluate them. Supposing that the field variable $f\left(k\right)$ is the one averaged over a narrow range $∆$ around the momentum k, and considering only the contribution from the narrow region near the peak where $dpdq$ is ${∆}^2$, we simplify (25) as follows:

$$\begin{array}{c}\hfill \left(fg\right)\left(k\right)=\frac{1}{{\left(2\pi \right)}^2}\frac{{∆}^2}{{\lambda }^2}\left[f\left(\lambda \right)g\left(k\right)+f\left(k\right)g\left(\lambda \right)\right].\end{array}$$

(31)

We also simplify (26) in the same way. The evolution equation derived in this way is written down in Appendix B, in which $T{R}_{\left[\lambda \right]}\left(k\right)=0$(A7) is a simplified version of the trace Equation (20) while $SP\left(k\right)=0$(A8) is the constraint Equation (21).

Hence, we have to solve a system of nonlinear evolution equations consisting of the four equations

$$\begin{array}{c}\hfill T{R}_{\left[\lambda \right]}\left(k\right)=0,SP\left(k\right)=0,T{R}_{\left[\lambda \right]}\left(\lambda \right)=0,SP\left(\lambda \right)=0,\end{array}$$

(32)

simultaneously, together with the homogeneous Equation (10). Here note that the last two $T{R}_{\left[\lambda \right]}\left(\lambda \right)=0$ and $SP\left(\lambda \right)=0$ form a system of equations representing the nonlinear evolution of the field with momentum $\lambda$, and can be solved with only these two equations. In this way, although originally we have to solve a multi-line coupled equation system consisting of many fields with different momentums, it is now reduced to a two-line system consisting only of the fields with momenta k and $\lambda$, which can be solved relatively easily.

Let us solve this simplified evolution equation to see how the nonlinear terms work. Here $∆$ is introduced as a phenomenological parameter that effectively expresses the strength of the nonlinear term, and its value is assumed to be approximately $∆/2\pi \sim \lambda$ because the coefficient of (31) becomes about 1.

In solving the nonlinear evolution Equation (32), we have to pay attention to the role of $SP=0$. As mentioned near the end of Section 3, this equation acts like a constraint condition connecting the initial stage where the running coupling constant $\overline{t}$ is small until the moment of the spacetime phase transition where it diverges. It indicates that the universe at ${\overline{t}}^2\to 0$ before inflation is in a state where only the fluctuation of $\Phi =\phi$ with $h=0$ exists. On the other hand, at the dynamical time $t_{\Lambda }$ when ${\overline{t}}^2$ diverges, the fluctuation changes to that satisfying $\Phi =h/3$.

To solve the evolution equation while preserving the constraint $SP=0$, we need to solve it as a boundary value problem. The initial time $t_i$ is set to ${10}^{-3}$, well before the Planck time, and the amplitude of $\Phi$ at that time is given by the square root of its scale-invariant power spectrum $\langle |\phi (t_i,k){|}^2\rangle ={\pi }^2/b_c$ obtained by Fourier transform of the two-point correlation function (4) [23]. Thus the initial value of $\Phi$ is set to

$$\begin{array}{c}\hfill \Phi (t_i,k)=\frac{\pi }{\sqrt{b_c}}\end{array}$$

(33)

and ${\partial }_t\Phi (t_i,k)={\partial }_t^2\Phi (t_i,k)={\partial }_t^3\Phi (t_i,k)=0$ for $k\ge \lambda$. In order to set the boundary condition on the field h, we rewrite it as $h=3(\Phi -X)$ and impose the boundary condition on X. Since the initial value of h is zero, $X(t_i,k)$ is equal to (33), while X should vanish at the dynamical time $t_{\Lambda }(=60)$, thus we set

$$\begin{array}{c}\hfill X(t_i,k)=\frac{\pi }{\sqrt{b_c}},X(t_{\Lambda },k)=0.\end{array}$$

(34)

To evaluate time evolution of the field with momentum k, we have to solve (32) with the boundary conditions (33) and (34) at k, together with those where k is replaced with $\lambda$. Further comments on how to solve it are summarized in Appendix B.

The fluctuation that was scale-invariant before inflation deforms and eventually reduces with time. Figure 4 shows the results of time evolution of the fluctuations solved in this way for various $∆$, in which one of the results (black) is already displayed in Figure 1 and Figure 2 together with an enlarged view near the phase transition. The line of $k=\lambda$ on the right end of each figure is the result of solving only the latter two equations in (32), while each k line is obtained by solving the four equations simultaneously.

Figure 5 plots the value of $\Phi$ at the time of the phase transition, and the square of it will give the primordial spectrum that is an initial condition of the subsequent time evolution. The right normalized panel shows that the nonlinear term has effects of flattening the spectrum even in the region $k>m$. In the previous linear equation results [22], the spectrum was greatly deformed at $k>m$, but it can be seen that the nonlinear effect suppresses the deformation up to the high momentum region, exceeding 3 times m.

Here, we examined the evolution equation simplified by introducing the parameter $∆$, but in fact it will be necessary to solve a system of evolution equations consisting of a very large number of the fields obtained by describing the momentum integrals of (27) and (28) by sums. In addition, although only nonlinear terms up to second order in the fields are considered here, it is believed that higher-order nonlinear terms contribute in higher momentum regions. In fact, here the linear approximation exhibits large undulations in the spectrum and some strange behavior, which suggests that the nonlinear effects given by the exponential function become essential in that region.

Why, then, do the nonlinear terms with the exponential factor play an important role in maintaining the scale invariance despite the presence of mass scale? The Einstein–Hilbert action having such a factor is one of the conformal fields in terms of conformal field theory. It is generally believed that correlation functions in the theory with it as a potential term exhibit a power-law behavior. It may be considered that this behavior is manifested.

Thus, we can expect that the scale invariance before inflation is maintained up to the relatively-high momentum region due to the nonlinear effects, and is conveyed to the present CMB. However, to what extent will such an effect be maintained? Could unexpected structures emerge in regions of even higher momentum [48]? The results here are the first step toward clarifying these matters.