Evolution of Vacuum Fluctuations of an Ultra-Light Massive Scalar Field generated during and before Inflation

We consider an ultra-light scalar field with a mass comparable to (or lighter than) the Hubble parameter of the present universe, and calculate the time evolution of the energy-momentum tensor of the vacuum fluctuations generated during and before inflation until the late-time radiation-dominated and matter-dominated universe. The equation of state changes from $w=1/3$ in the early universe to $w=-1$ at present, and it can give a candidate for the dark energy that we observe today. It then oscillates between $w=-1$ and $1$ with the amplitude of the energy density decaying as $a^{-3}$. If the fluctuations are generated during ordinary inflation with the Hubble parameter $H_I \lesssim 10^{-5} M_{\rm Pl}$, where $M_{\rm Pl}$ is the reduced Planck scale, we need a very large e-folding number $N \gtrsim 10^{12}$ to explain the present dark energy of the order of $10^{-3} {\rm eV}$. If a Planckian universe with a large Hubble parameter $H_P \sim M_{\rm Pl}$ existed before the ordinary inflation, an e-folding number $N \sim 240$ of the Planckian inflation is sufficient.


Introduction
Our universe is well described by the spatially-flat ΛCDM model. According to the PLANCK 2013 results [1], only 5.1% of the energy density is attributed to known form of baryonic matter, while 26.8% is attributed to cold dark matter, and 68.3% to dark energy. Although its equation of state w(= p/ρ) = −1 seems like that of vacuum energy of quantum fields, there is no reasonable explanation for its magnitude, ρ DE = 3(M Pl H 0 ) 2 Ω Λ ∼ (2.2meV) 4 . Here, M Pl = (8πG N ) −1/2 ∼ 2.4 × 10 30 meV is the (reduced) Planck scale and H 0 ∼ 1.4 × 10 −30 meV is the current Hubble parameter. ρ DE is by far smaller than the expected magnitude of vacuum energy Λ 4 in a theory with an ultra-violet (UV) cutoff Λ. If we take Λ to be M Pl , ρ DE is smaller than Λ 4 by more than 120 orders of magnitude. This is the cosmological constant problem [2].
On the other hand, we may try to explain the dark energy as the Casimir energy in the current universe. In particular, vacuum energy for fluctuations of massless fields in de Sitter background with Hubble parameter H is of order H 4 , and has w = −1 [3,4]. Since our present universe is close to de Sitter space, one may wonder if dark energy can be explained as vacuum energy in de Sitter with the current Hubble parameter H 0 , but this does not seem to be plausible. Dark energy that we observe M 2 Pl H 2 0 is much larger than the expected contribution from a single field H 4 0 . However, H 0 is not the only dimensionful quantity that affects the renormalized energy-momentum tensor (EMT). To compute expectation values of fluctuations, we need to specify the vacuum state. This may depend on global properties of the geometry and the whole history of the universe, thus different scales might be introduced in the problem. There is by now strong evidence [5] that there has been a period of inflation with the Hubble parameter H I much larger than H 0 . It would be reasonable to take the vacuum to be the Bunch-Davies vacuum [3] for de Sitter space with H I . Fluctuations of massless scalar in de Sitter background is of order H I . Fluctuations are frozen (remain constant) outside the Hubble radius (see e.g. [6]), thus, infra-red (IR) modes could have a large value in the universe after inflation. In fact, these fluctuations are considered to be the origin of the fluctuations of cosmic microwave background (CMB) that is observed today [1,5].
In [7] and [8], the time evolution of the EMT is calculated for a minimally coupled massless scalar field. The fluctuations are generated during the inflationary universe and evolve until the late-time universe of the radiationdominated (RD) and the matter-dominated (MD) era. The equation of state w approaches w = 1/3 and w = 0 in the RD and MD periods respectively. The magnitude of the present energy density is of order H 2 I H 2 0 , and is still much smaller than that of the dark energy in our universe. The analysis is extended to a non-minimally coupled scalar field in [9].
In order to circumvent the smallness of energy density, we considered a double inflation model in [8]. We assumed that there was an inflation with a Hubble parameter H P of the order of the Planck scale M Pl (which should be natural in the Starobinsky's inflation [10], for example) before the usual inflation with H I started. We fix the initial condition of the fields in the Planckian inflation period by taking Bunch-Davies vacuum with the Hubble H P , and study the time evolution afterwards. In this case, IR mode is enhanced to H P , and the present value of vacuum energy becomes of order H 2 P H 2 0 . In order to make these large fluctuations consistent with the observed value of CMB fluctuations, the enhancement needs to be restricted in the far IR modes whose wave lengths are larger than the current Hubble radius H −1 0 . Since the enhanced modes are still out of the horizon now and the dominant contribution to the EMT is given by the spacial-derivative parts |∇φ| 2 , the equation of state is given by w = −1/3 instead of w = 1/3 or w = 0 as in the ordinary inflation.
In the present paper, we extend our previous analysis in [8] to an (ultra-light) massive scalar field. If the mass is smaller than the Hubble parameter m ≪ H I , the wave function receives large enhancement during the inflation as in the massless case. The low-momentum modes remain frozen until the mass becomes larger than the Hubble parameter in the late universe when the field starts to oscillate. Such ultra-light scalars have been studied extensively as candidates for the dark matter (see [11] and papers citing it). The scenario has attracted renewed interest after [12]. For m 10 −24 eV, it is indistinguishable from the standard cold dark matter [13]. For lighter mass 10 −32 eV m 10 −25.5 eV, its abundance is strongly constrained by the CMB and galaxy clustering data. If its mass is smaller than the current Hubble parameter H 0 ∼ 10 −33 eV, it can be a candidate for the dark energy at present [14,15]. If the light particle is an axion-like particle (ALP), the initial value of the field is set by the misalignment mechanism, so the amplitude of the energy density can be chosen by hand. A similar idea is given by the quintessence scenario [16], where the initial value of the field is also set by hand. In this paper, we investigate a possibility that the initial field value is dynamically determined by the fluctuations generated during the primordial inflation. The magnitude of the energy density is expected to be (m 2 /2)(H I /2π) 2 N where N is an e-folding number of the inflation. 1 For m H 0 , we may need a large N or/and a large H I to make it comparable to the current dark energy of the unvierse.
The purpose of the present paper is to investigate detailed behavior of the time evolution of the EMT of the vacuum fluctuations generated during the inflation until the late-time universe of the RD and MD eras.
The paper is organized as follows. In section 2, we solve the equation of motion of a massive scalar field in the history of the universe with the inflation period followed by the RD universe. In section 3, we calculate the EMT and study its time evolution. In section 4, we obtain the conditions that the vacuum fluctuation of the ultra-light scalar explains the dark energy at present. The last section 5 is devoted to conclusions and discussions. In Appendix A, we calculate the time evolution of the EMT by using the zero-momentum approximation. The time evolution in the MD period is given by using this approximation.

Massive scalars in expanding universe
Our universe is approximated by the Robertson-Walker spacetime with the inflation, RD (radiation-dominated) and MD (matter-dominated) periods. In this section, we focus on the first two stages of the universe and study the detailed behaviors of the wave function of a massive scalar field. 2 The metric is given by ds 2 = a(η) 2 dη 2 − (dx i ) 2 in terms of the conformal coordinates (η, x i ), and the scale factor a(η) is given by where η ini and η 1 denote the beginning and the end of the inflation, and η 2 the beginning of the RD period. The continuity conditions for the scale factor a and its derivative a ′ = ∂ η a at the boundary of the inflation and the RD periods give the relations (2. 2) The Hubble parameter H = a ′ /a 2 is given by different physical setting in which we suppose a period of pre-inflation with a large Hubble parameter HP ∼ M Pl before the ordinary inflation starts. In this case, we do not need an extremely large e-folding number to explain the dark energy in the present universe. 2 In [18], particle content and the degree of classicality were studied in a similar background (de Sitter in inflation, followed by RD, and late-time de Sitter).
The CMB fluctuations give a constraint H I < 3.6 × 10 −5 M Pl . We consider a minimally coupled massive scalar field φ with a mass m. Quantum fields are expanded as where the mode functions u k (η) with the comoving momentum k are the solutions of the equation of motion, ( + m 2 )u = 0, and are chosen to asymptote to positive-frequency modes in the remote past. A vacuum is then defined by a k |0 = 0. The vacuum |0 , which is an in-state, evolves as η increases, and if an adiabatic condition is broken, the state gets excited above an adiabatic ground state at each moment η.

Wave functions in the inflationary period
In the inflationary period, a solution of the wave equation is given by where H (1) ν is the Hankel function of the first kind and for m ≪ H I . Another solution is given by its complex conjugate. Using an asymptotic form of the Hankel function at |z| → ∞, For η → 0, on the other hand, the wave function can be approximated by using the expansion formula of the Hankel function H (2.13) It is valid when k|η| = k phy /H I < 1 is satisfied where the physical momentum in the inflationary period is given by k phy = k/a = kH I |η|. In the massless limit, the wave function u BD (η) = χ BD (η)/a Inf (η) becomes almost constant in time η, and behaves as k −3/2 . It is consistent with the η → 0 behavior of the exact massless wave function (2.14)

Wave functions in the RD period
In the RD period, the wave equation (2.5) is written as where we have rescaled the parameters η and k as Since the Hubble parameter is given by H = H I (η 1 /η) 2 , the new parameters x and q have the following physical meaning: where k phy = k/a RD is the physical momentum in the RD period. The equation (2.15) is nothing but the Schrödinger equation in an inverted harmonic oscillator with the Hamiltonian H = −∂ 2 x − x 2 /4 and the energy eigenvalue q 2 . The exact solution can be obtained by analytical continuation of the solution in the harmonic potential, and is given by where c = (2m) −1/4 . D n (z) is the parabolic cylinder function satisfying the Weber differential equation It has the following asymptotic behavior The exact solution (2.18) approaches the WKB-approximated one (2.24) for x → ∞, as will be seen below.
If the adiabaticity condition is satisfied, the exact wave function is approximated by the WKB wave function, The integral of ω(x) is given by where the x-independent terms are fixed by comparing with the asymptotic behavior of the exact wave function. By taking c = (2m) −1/4 , the WKB wave function (2.21) is normalized by the KG inner product with respect to η-derivative as (χ WKB , χ WKB ) η = 1. The adiabaticity condition is given by For q < x/2, namely k phy < m, the WKB wave function (2.21) is further approximated as Then the exact wave function (2.18) can be shown to asymptote to the WKB wave function by using (2.20). Note also that x 2 /4 = m/(2H) = mt where t is the physical time in the RD period, and the wave function oscillates as e −imt . For q > x/2 (i.e. k phy > m), the WKB wave function is reduced to the plane wave  Figure 1: The wave function in the RD period changes its behavior depending on the parameters x = 2m/H (time) and q = k phy / √ 2mH (momentum). The thick solid curve represents ǫ = 1; outside the curve, the WKB approximation is valid. The vertical line is x = √ 2, which corresponds to m = H due to (2.17); in the left to the line, m < H. The dashed curve is qx = k phy /H = 1; under the curve, k phy < H. The dot-dashed line is q = x/2, which corresponds to k phy = m; under the line, k phy < m.
Since qx = k phy /H = 2k phy t, the oscillation behavior e −2ik phy t is controlled by the momentum, not by the mass.
For small x, we can approximate the exact wave function χ EX of (2.18) by a power series of x. The parabolic cylinder function is expanded as where the first two coefficient functions are given by The rest of the coefficient functions c r (q 2 ) are related to them as, e.g.
If we take the highest powers of q 2 in each c r (q 2 )/c 0 (q 2 ) and c r (q 2 )/c 1 (q 2 ), the parabolic cylinder function is expanded as For small x, it can be approximated as Since the above expansion turns out to be a double expansion of x 4 and (qx) 2 , it is a good approximation for small q 2 at x < 1. At the same time, since we took the highest orders of q 2 in each c r (q 2 ), it should also be a good approximation for large q 2 . Indeed, we have confirmed numerically that the remaining terms in the square bracket in (2.29) decrease and oscillate as a function of qx. It will be related to the fact that the expansion of (qx) 2 is an alternative power series. Hence, (2.30) gives a good approximation for small x 1, irrespective of the magnitude of q.

Determination of the Bogoliubov coefficients
We now solve the wave equation throughout the inflationary and the RD periods by imposing the BD initial condition. The wave function can be written as where χ BD (η) and χ EX (η) are defined in (2.8) and (2.18). 3 The Bogoliubov coefficients A(k) and B(k) can be determined by imposing continuity of χ and ∂ η χ at the boundary of the inflationary and the RD periods, and are given by The Bologiubov coefficients A(k) and B(k) in the IR region k|η 1 | < 1 are calculated by using the wave function χ BD (η) in (2.13). Its derivative is written as where (m/H I ) ≪ 1 was used in the second equality, and (2.32) become where kη 2 < 1 and √ 2mη 2 ≪ 1 were used in the second equality.

Behaviors of the wave functions with BD initial condition
By using the Bogoliubov coefficients A(k) and B(k), the wave function u(η) in the RD period with the BD initial condition is given by where c 0 is defined in (2.27). It is one of the main results of the paper, and gives a starting point to calculate the EMT in the RD period. This expression is valid for IR modes with momentum k|η 1 | = kη 2 < 1. For UV modes with k|η 1 | > 1, the wave function is not enhanced like this. For early times x = √ 2mη = 2m/H 1, the exact solution χ EX (η) is approximated by (2.30), and the wave function u(η) is simplified as (2.37) The momentum is set to be k = 0. The wave function is frozen until x ∼ √ 2 (i.e. m ∼ H) and then starts oscillation.
Here we have used the identity which can be proved by using the property of the KG inner product (χ EX , χ EX ) η = 1 at η = 0. For kη < 1, u(η) becomes almost constant, i.e. frozen. For later times x √ 2, the WKB approximation becomes valid, and the wave function u(η) is written as The WKB approximation is also valid for the UV modes (outside of the semicircle in Figure 1) even at early times x < √ 2. As we saw previously, the WKB wave function has two different behaviors, (2.25) at q > x/2 (k phy > m), and (2.24) at q < x/2 (k phy < m). In the region q > x/2, it can be shown numerically that the x-independent phase in the square bracket in (2.39) vanishes for large q > 1 and (2.39) becomes identical with (2.37). On the contrary, for q < x/2, the WKB wave function (2.24) oscillates as χ WKB ∼ e −ix 2 /4 = e −imt .
In Figure 2, the behavior of the wave function of the exact solution in (2.36) is plotted for q = 0. As expected, it is almost constant (frozen) for x √ 2, and slowly decreases with an oscillation for x √ 2.

Evolution of energy-momentum tensor
In this section, we calculate the EMT (energy-momentum tensor) in the RD period. The physical results that are relevant to the application to the dark energy are summarized in section 4. The vacuum expectation value of the EMT is given by The first term in ρ and p is the time-derivative term, which gives the equation of state w = p/ρ = 1. The second term is the spacial-derivative term, and gives w = −1/3. The third term is a contribution of the mass term, and gives w = −1.
The integrals (3.1) and (3.2) are divergent in the UV region, and need subtraction of the UV divergences. After the subtraction is performed, we simply cut off the integral at the momentum k = 1/|η 1 |, since only those IR modes are enhanced during the inflationary period. (For more details, see Section 6 of ref. [8].) The conformal anomaly is given by the subtraction term, but its contribution to the EMT is of the order of H 4 and is negligibly small compared to the vacuum fluctuations generated during the inflation of the order of H 2 I H 2 . So we will not consider it in the present paper.

Massless case
Before investigating the massive scalar field, let us first summarize the evolution of the EMT in the massless case studied in [7,8]. The wave function with the BD initial condition is given by for k|η 1 | < 1. Then the energy and pressure densities become with f (y) = sin(y)/y. The k-integral is performed for momenta with which the wave function is enhanced during the inflation. Hence, if the inflation continues during η ∈ [η ini , η 1 ], the integral region of the comoving momentum is restricted There is no IR divergence in (3.4), and the IR cutoff does not play an important role.
In the RD period, some of the enhanced modes enter the horizon again. The modes with k < 1/η are still out of the horizon and are frozen. Then the time-derivative term, the first term in (3.4), vanishes and the spacial-derivative term, the second term in (3.4), gives On the other hand, the UV modes with k > 1/η have already entered the horizon and the wave functions are time dependent. In performing the k-integration over [1/η, 1/η 2 ], we estimate the oscillating integrals as Then the energy and the pressure densities contributed from the UV modes become Here we dropped higher order terms with respect to (kη) −1 . We also defined an e-folding during the RD period, where BB stands for the big bang, namely, a BB means the scale factor at the beginning of the RD period or the end of inflation. N RD represents the number of degrees of freedom that were enhanced during the inflation and have already entered the horizon. The total energy density is given by ρ = ρ IR + ρ UV , but the UV contribution (3.7) dominates the IR part (3.5) when N RD 1. Thus the equation of state of the vacuum fluctuations generated during the inflation approaches w = 1/3 soon after the RD period begins.

Massive case at early times
We now study the evolution of the EMT in the massive case. In terms of the variables x and q, (3.1) and (3.2) become 2m/H I correspond to the beginning and the end of the inflation period. The integration region represents the momentum region where the wave function is amplified.
We first consider the behavior of the EMT at early times x = 2m/H < 1. The wave function is approximated by (2.37), and the EMT becomes where f (y) = sin(y)/y. The kinetic terms, (∂ y f (y)) 2 + f (y) 2 in the square brackets, are the same as those of the massless case (3.4) except for an additional power of the momentum q depending on m. But since the kinetic terms are IR convergent, it does not affect the integrals much. Hence, the contribution from the kinetic terms to (3.11) is given by the same form as in the massless case (3.7). The contribution from the mass term, namely the term proportional to (x 2 /4q 2 )f (y) 2 , would be IR divergent if we used the massless wave function.
However, the additional power of the momentum, q 2 3 m H I 2 , which comes from the mass deformation of the wave function, reduces the IR divergence. In terms of k and η, the mass term contribution is written as Since the function f (kη) starts decreasing at k ∼ 1/η, it gives an effective UV cut off, and we have Here, we used m ≪ H I and defined an effective e-folding 14) where is an e-folding number during the inflation period and N RD is that of the RD period defined in (3.8). Thus N eff represents the number of modes which were enhanced during the inflation period, and are still outside the horizon and frozen at time η.
For sufficiently large N eff , (3.13) becomes 3(H I /2π) 4 . It is nothing but the thermal equilibrium energy at the de Sitter temperature T = H I /2π, and independent of m. On the contrary, (3.13) can be approximated by is also given by a simple physical argument that in the de-Sitter spacetime an ultra-light field experiences the Brownian motion at temperature T = H I /2π. The growth of fluctuation since the initial time in de Sitter space was studied in [19,20,21,22]. 4 The pressure density is p mass = −ρ mass and it gives a candidate for the dark energy. As will be shown in section 4, it can explain the present dark energy if N eff satisfies N eff 24π 2 (M Pl /H I ) 2 ∼ 10 12 for H I ∼ 10 −5 M Pl .

Massive case at late times
We then consider the behaviors of the EMT at late times x = 2m/H > √ 2. In this parameter region, the WKB approximation becomes valid, and we can use the wave function (2.39). In performing the momentum integration for the EMT, we divide the integration region into two; the UV region q > x/2 which corresponds to k phy > m, and the IR region q < x/2 corresponding to k phy < m.
In the UV region with q > x/2, the WKB wave function is reduced to the plane wave ∝ e −iqx in (2.25). As we discussed below eq.(2.39), the WKB wave function with the BD initial condition (2.39) becomes identical with (2.37). Then the EMT reduces to (3.11) except for the integration region, where the lower bound of the integration region is replaced by q = x/2. Since q > x/2 in the UV region, the terms proportional to x 2 /4q 2 in (3.11) can be neglected. Then the UV contribution to (3.11) becomes where we used ln 2 It is slightly different from N RD defined in (3.8), since the integration is cut off at q = x/2 (k phy = m), not at q = 1/x (k phy = H). However, the difference is negligible as far as we consider the region x ∼ √ 2, or m ∼ H.
We next consider the IR region, q < x/2 (i.e. k phy < m). In this region, the q-dependence of the WKB wave function (2.39) is given by u ∝ q −3/2+1/3(m/H I ) 2 . The other factors in (2.39) have mild q-dependences, and can be approximated by those at q = 0. 5 Then the energy density (3.9) can be evaluated as Here, we dropped the spacial-derivative terms proportional to q 2 in (3.9) since we are considering the IR region q < x/2. Note that the q-and x-dependences of the integrand have been separated in (3.19). The q-integration can be performed as where we used In this approximation, the time dependence of ρ IR is represented by the zero-momentum mode. In Appendix A, we give a simple derivation of the time evolution in the zero-momentum approximation. In Figures 3 and 4, we numerically calculate the time evolution of ρ without using such an approximation for the IR modes.
at m ∼ H, or x ∼ √ 2. The square brakcket in (3.19) is calculated as (3.23) When the higher order corrections to the WKB approximation (2.21) are included, the wave function in (3.19) is modeifeid as Then the second term with x −2 in (3.22) remains unchanged, but the third term with x −4 receives modifications. Thus, we will drop it. Therefore, (3.19) becomes The time dependence is dominantly given by the overall factor H 3/2 ∼ a −3 with an additional oscillation given by the square bracket term. Note that it is a nonanalytic function of m and the behaviour cannot be obtained from the massless theory by perturbation with respect to the mass. An interpolating solution between the early-time behavior (3.16) and the late-time behavior (3.25) can be obtained by the zero-momentum wave function given in Appendix A. Similarly, one can estimate the pressure density. From (3.10), one obtains almost the same equation as in (3.19), but with the second term in the square bracket having the minus sign. Accordingly, (3.22) is replaced by and we have Note that the EMTs, (3.25) and (3.27), at late times decrease as a −3 . The behaviour is consistent with our knowledge that once the scalar field φ starts oscillation in the quadratic potential mφ 2 /2, it behaves as an ordinary dust. However, the pressure is non-vanishing and the equation of state oscillates, Let us see how these behaviors are consistent with the conservation law of the EMT. The energy and the pressure densities behave as It can be easily shown that they satisfy the conservation law of the EMT with ∂ x a/a = 1/x up to the order x −4 . Hence, the energy density decreases as a −3 as if it were the pressureless dust though the pressure is nonvanishing but oscillating. To see the consistency at the order x −6 and higher, we need to include higher order corrections to the WKB approximation in (2.21). The ratio of the UV contribution (3.17) to the IR contribution (3.25) is given by As we saw at the end of the previous section, N RD ∼ 60, and a very large e-folding N eff 10 12 is necessary. Note also that m > H is satisfied for the late times. Therefore the r.h.s. in (3.31) is much smaller than one, and the IR contribution ρ IR gives a dominant contribution to the EMT compared to ρ UV .
In Figure 3, we numerically performed the q-integrations of (3.9) and (3.10) for the WKB wave-function (2.39) with (2.21), and plotted the time evolution of the energy and the pressure densities in the left. The equation of state is plotted in the right. The energy and pressure densities are normalized by The integration region is taken as q ∈ [q min , x/2] to obtain the IR contribution. The lower bound corresponds to as (3.25) and (3.27). On the other hand, if we take the IR cutoff larger, the approximation that the oscillation frequency of the integrand is independent of the momentum in the IR region is invalidiated. Then, by summing various modes with different momenta, the oscillating behavior in the time direction is incoherently averaged and is expected to be diminished. Indeed, in the lower figures, the magnitude of p, and accordingly, w decreases. However, even in these large q min , the oscillating behavior of w still remains. In the real setting, as we saw before, we need large N eff ∼ 10 12 , and thus small q min ∼ e −10 12 is required. Hence the equation of state w oscillates between w = −1 and 1 rather than diminishes.

Numerical results for the evolution of EMT
In order to see how the early and late times behaviors are smoothly connected, we evaluate the evolution of the EMT by using the exact solution (2.36) 6 instead  Figure 4: Similar plots to Figure 3, but using the exact wave function in (2.36). Note that the the horizontal axis runs from x = 0, unlike from x = 2 in Figure 3. The integration region over q is taken to be q ∈ [10 −5 , 10 −2 ]. At early times, the IR wave function is frozen and the IR part of the EMT behaves like the dark energy w = −1.
of using the WKB approximation. We insert the wave function (2.36) into the EMTs, (3.9) and (3.10), and perform numerical integration over q. In order to see the contributions from the IR modes, we set the integration region as q ∈ [10 −5 , 10 −2 ]. The results are shown in Figure 4, where ρ IR and p IR are normalized by the factor (3.32) as in Figure 3. For x > √ 2, the results of the exact solution in Figure 4 agree with those of the WKB wave function given in Figure 3. For x < 1, w approaches −1, which agrees with the previous results given in section 3.2. The early time behavior at x < 1 is smoothly connected with the late time behavior at x > √ 2. For technical reasons, the numerical integrations are performed in a restricted region of q > 10 −5 and q < 10 3 . In order to integrate over q ∈ [q min , q max ] where q min ∼ e −N eff < e −10 12 and q max = x −1 1 = (2m/H I ) −1/2 ∼ 10 30 , we use the analytical results based on the approximations of the wave function discussed in the previous sections.

Vacuum fluctuations as dark energy
We now investigate possibilities for the vacuum fluctuations of the ultra-light scalar to explain the dark energy at present. We give conditions for the mass and the e-folding number, and then discuss how the EMT evolves through the RD and MD periods. We consider two scenarios. In section 4.1, the ordinary inflation model is discussed. In section 4.2, a double inflation model is considered where we assume another inflation with a larger Hubble parameter H P ∼ M Pl before the ordinary inflation starts.

Ordinary inflation model
We first summarize how the EMT of an ultra-light scalar field evolves in the RD period. In the ordinary inflation model, the enhanced mode during the inflation with the largest momentum soon enter the horizon after the inflation ends. Hence the EMT is given by a sum of UV and IR contributions. At early times with m H, as discussed in section 3.2, the EMT is approximately given by ρ = ρ UV + ρ IR and p = p UV + p IR where The UV modes have already entered the horizon. N RD represents the number of UV degrees of freedom. Hence it is time dependent. For the UV modes, the kinetic terms in the EMT mainly contribute and w = 1/3 is obtained. The IR modes are still out of the horizon and frozen. Hence the mass term mainly contributes to the EMT, and we have w = −1.
If the mass of the scalar field is heavier than the Hubble parameter at the matter-radiation equality; m > H eq = 10 −28 eV, the condition m H becomes satisfied in the late RD period, and the EMT is described by the late time behaviors discussed in section 3.3. The coherent oscillation (the motion of the zero mode) gradually starts and the behavior (4.2) of the EMT is changed to (3.25) and (3.27): where s and c are defined in (3.23) and The amplitudes are proportional to √ mH 3/2 and decaying as H 3/2 ∝ a −3 . Hence it is like the dust but the pressure density is non-vanishing and oscillating. The UV contribution (4.1) remains unchanged. Now let us study the evolution of the EMT in the MD period. If the mass is lighter than H eq , m < 10 −28 eV, the coherent oscillation of the zero-mode of the scalar field has not yet started at the beginning of the MD period. In the early times of the MD period when the condition m H is satisfied, the EMT is again written as a sum of the UV and IR contributions, ρ = ρ UV + ρ IR . The UV part comes from the modes that have already entered the horizon, and is given by eq. (7.20) of ref. [8] where a eq and a are the scale factor at the matter-radiation equality and at the time in the MD period. N RD = ln(a eq /a BB ) is the e-folding number during the RD period and is constant in time. Compared to the IR contributions discussed below, this UV contribution ρ UV becomes negligible because of the factor (a eq /a). The IR part in the EMT has contributions from the kinetic term and the mass term, ρ IR = ρ IR,kin + ρ IR,mass . They are evaluated as ρ IR,kin was obtained in eq. (7.18) in ref. [8]. ρ IR,mass is the same as in (4.2), but N eff takes a slightly different value since it represents the number of degrees of freedom that are still out of the cosmological horizon, and is time dependent. Note that ρ IR,kin in (4.6) receives larger contributions from the modes with momenta k ∼ η −1 (i.e., k phy ∼ H). In contrast, ρ IR,mass in (4.7) has dominant contributions from the modes with much lower momenta. It is amusing that a single ultra-light scalar simultaneously contains the dark-energy-like component ρ IR,mass and the dark-matter-like component ρ IR,kin . At later times in the MD period when the condition m H is satisfied, the coherent oscillation starts and the IR contribution is changed. As shown in appendix A, if time evolution is represented by the zero-momentum approximation, the energy and pressure densities are given by (A.9) and (A.10): where s and c are defined in (A.11) and Note that the amplitude of the energy density decreases as H 2 ∝ a −3 . The interpolating solution between the early-time behavior (4.7) and the late-time behavior (4.8) is obtained in (A.9) and (A.10) in Appendix A.
A dark energy candidate is given by (4.7). We need three conditions for (4.7) to explain the dark energy in the present universe; where 0 denotes the present time. The first condition states that the present time corresponds to the early times before the coherent oscillation of the bosonic field starts. The second condition requires that the mass term contribution with w = −1 dominates over the kinetic term contribution with w = 0. The second condition gives a lower bound for the mass (or a lower bound for the effective e-folding). Combining these two conditions we have (C1, C2) : The third condtion is necessary if the observed magnitude of the present dark energy is given by (4.7). It can be written as (C3) : Inserting (C3) into (C1,C2), we have the following conditions for H I and N eff : In the ordinary inflation, we already have a constraint from the CMB observation that H I < 3.6 × 10 −5 M Pl . Then the second condition in (4.13) is already satisfied. The first one requires quite a large e-folding, N eff > 1.8 × 10 11 . (4.14) This may indicate that the observed universe with the Hubble radius 1/H 0 is embedded in a huge universe whose size is e 1.8×10 11 times larger.
In Figure 5, we show the time evolution of the EMT in the RD and MD periods. The upper figures plot the energy density ρ, divided by the critical value ρ cr , while the lower figures plot the equation of state w = p/ρ. We use m/H for the horizontal axis to denote time evolution. We used the following numerical values: the Planck scale M Pl = 2.4 × 10 30 meV, the present Hubble parameter H 0 = 1.4 × 10 −30 meV, and the z-factor at the matter-radiation equality z eq = 3.4 × 10 3 . Thus the Hubble parameter at the equality is given by H eq = H 0 z 3/2 eq = 2.8 × 10 25 meV. In drawing the figures, we chose the Hubble At early times, the UV contribution to the kinetic term ρ UV is dominant and gives the equation of state w = 1/3, while its magnitude is much smaller than the critical value. As time passes, the IR contribution to the mass term ρ IR,mass grows and dominates ρ UV at some time, when w changes to −1. For the parameters we chose, the transition between the above two behaviors occurs, accidentally, at around the same time of the matter-radiation equality. 8 As time passes further, the energy density approaches the critical value, which gives the present dark energy, and then w begins to oscillate. As the lower figures in Figure 5 show, the oscillating behavior already begins at present m/H = 1 and the equation of state is given by w ∼ −0.9. If we choose a smaller mass like m ∼ 0.1H 0 , the equation of state w = −1 can be realized at present, but we need a 100 times larger N eff. Such difference will be detected in the future observations.
It is also interesting to note that the energy density goes over the critical value when m > H, as shown in the upper figures in Figure 5. Indeed, the coefficient of (4.9) is 9/4 times larger than that of (4.7). Then the condition of (C3) in (4.10) necessarily leads to a situation that the energy density exceeds the critical density of the background universe in future. So we need take back reactions into account to extrapolate our analysis to obtain the behaviors of the future universe.
In [29], it was shown that the vacuum energy of a quantum field drives de Sitter expansion in the inflation period, by studying the back reaction in a selfconsistent way. It is also an interesting theoretical problem to study the latetime behavior to understand the fate of the universe. In this case, the interplay between the scale factor and the IR behavior of the wave function, with the Bunch-Davies initial condition, determines the dynamics of the universe selfconsistently.

Double inflation model
Let us now consider another possibility that the EMT of the ultra-light scalar field is enhanced due to the fluctuations created before the ordinary inflation period. As a simple example, we consider a cosmic model with two inflationary periods: the ordinary inflation with the Hubble parameter H I and a preinflation with a larger H P before the ordinary inflation. A similar model was studied to obtain a modified CMB spectrum (see e.g. [30]). It will be reasonable that such a period existed before the ordinary inflation, because in the very early time of the Planck scale, quantum gravity effects possibly generate a large-Hubble de-Sitter expansion. Some concrete examples are the Starobinski type of inflation [10] and the eternal inflation [31,32,33,34], where our universe is surrounded by the region with a larger Hubble parameter.
In addition to the Hubble parameter in the preinflation period H P , the double inflation model has another important parameter, i.e. the conformal time η * when the preinflation period ended and the ordinary inflation started. As studied in ref. [8], the wave function is enhanced to a larger amplitude of the order of H P during the preinflation period, and the enhanced modes are restricted within the momentum region k < 1/|η * |. If |η * | is larger than the current conformal time η 0 , all the enhanced modes are still outside of the current horizon, and do not affect the CMB data. In such a case, H P is not constrained by the CMB observation, and can be taken as large as the Planck scale M Pl .
We now study the time evolution of the EMT. We first consider the RD period. The EMT is given by a sum ρ = ρ UV +ρ IR . Here ρ UV is the contribution to ρ from the UV modes which are enhanced during the ordinary inflation, but not during the preinflation. Hence it is given by the same equation as (4.1), On the other hand, the IR part ρ IR is the contribution from the IR modes that are enhanced largely during the preinflation with H P . ρ IR is given by a sum of the mass and the kinetic terms ρ IR = ρ IR,mass + ρ IR,kin . Before the coherent oscillation starts, i.e., when m < H, the mass term in the EMT becomes ρ IR,mass = 1 8π 2 H 2 P m 2 N Preinf , p IR,mass = −ρ IR,mass , (4.16) where is an e-folding number during the preinflation period. Here, η ini and η * denote the conformal time at the beginning and the end of the preinflation period. The kinetic term in the EMT becomes (eq. (8.22) in ref. [8]), As mentioned above, we assume that all the enhanced modes are out of the current horizon. Then they are also out of horizon and frozen always in the past. Hence, the time-derivative term in the EMT vanishes, and the spatialderivative term gives (4.18) with w = −1/3.
In the MD period, the UV contribution becomes negligibly small as in (4.5 At late times (m H), the IR contribution becomes (4.8) and (4.9), with H I and N eff replaced by H P and N Preinf . The oscillating behavior of w is obtained again.
As in the ordinary inflation case, (4.19) gives a dark energy candidate. In order to explain the dark energy in the present universe, the three conditions (4.10) are required. The first and the second conditions give (C1, C2) : (4.24) Since the enhanced modes must be outside of the horizon in the present universe, |η * | > η 0 /(2π) needs to be satisfied, which gives (1/16)(η 0 /η * ) 2 < π 2 /4. Then, a large Hubble parameter as H P ∼ M pl satisfies the second condition. The first condition requires that the number of e-folding during the preinflation must satisfy N Preinf > 24π 2 ∼ 2.4 × 10 2 . Compared to the ordinary inflation, the e-folding number is not necessary to be so large as (4.14).  Figure 6: Similar plots to Figure 5, but for the vacuum fluctuations generated by the preinflation. The same numerical values are used as in Figure 5. For the additional parameters, we used H P = M Pl , |η * | = η 0 , and thus N Preinf = 2.4 × 10 2 . Compared to Figure 5, a new state with w = −1/3 appears during the transition from w = 1/3 to w = −1.
In Figure 6, we show the time evolution of the energy density ρ/ρ cr and the equation of state w = p/ρ in the RD and MD periods. We take the same parameters as in Figure 5. For the additional parameters, we use H P = M Pl and |η * | = η 0 . Then N Preinf = 2.4 × 10 2 is required by (4.23).
At early times, the UV contribution to the kinetic term ρ UV is dominant and gives w = 1/3, while its magnitude is much smaller than the critical value. As time passes, the IR contribution to the kinetic term ρ IR,kin dominates, and the era with w = −1/3 starts. As time passes further, the IR contribution to the mass term ρ IR,mass dominates, and the era with w = −1 starts. At later times, when m H, the era with oscillating w starts.
The existence of the intermediate stage with w = −1/3 depends on values of the free parameters we take. The transition time m/H from w = 1/3 to w = −1/3, i.e., when ρ IR,kin dominates ρ UV , can be shown to be proportional to m(H I η * /H P ) 2 , by comparing (4.15) and (4.18). On the other hand, the time m/H when ρ IR,mass dominates over ρ UV is proportional to mH I , by comparing (4.15) and (4.16) with (4.23). Hence, if we choose a larger η * and/or a smaller H P , the era with w = −1 stats before w = −1/3 might start; the intermediate stage with w = −1/3 does not arise.

Conclusions and discussions
In this paper, we calculated the time evolution of the energy-momentum tensor of an ultra-light scalar field with a mass m 10 −33 eV. In the case of axion-like particles, the initial condition is set by hand by the misalignment mechanism. We instead assume that the fluctuations generated during de Sitter expansion in the primordial inflation gave the initial condition of the amplitude of the vacuum energy. If the fluctuations are enhanced during the ordinary inflation with the Hubble parameter H I ∼ 10 −5 M Pl , a very large e-folding N eff ∼ 10 12 is necessary to explain the dark energy at present. But if we consider a cosmic history with another Planckian universe with a large Hubble parameter H P ∼ M Pl before the ordinary inflation, much smaller number of e-folding N eff ∼ 240 during the preinflation is sufficient. We furthermore calculated how the dark energy evolve in future. The amplitude decreases as a −3 where a is the scale factor and it is the same as the pressureless dust. However, the pressure does not vanish and the equation of state oscillates between w = −1 and 1 with a very large oscillation period 1/m. We showed that this behavior is consistent with the equation of motion. If the mass is a bit larger than the current Hubble, e.g. m = 10 −32 eV, such an oscillatory behavior may be detectable (see e.g. [35]). In most studies of the quintessence scenario, the classical equation for the zeromode is used to compare with the observational data. It is interesting to extend the analyses to include non-zero momentum modes discussed in this paper and to give more detailed observational constraints on the model parameters.
Another important issue not discussed in the present paper is an effect of interactions. In the de Sitter expanding universe, we often encounter with IR divergences ln(k|η|) when we calculate loop corrections of various quantities (see e.g. [36]). The IR divergences are related to the secular time growth of these quantities in |η| → 0 limit and in many cases they can be resummed. For example, the secular growth in a massless scalar with λφ 4 interactions can be cured by resumming the logarithmic factors so that the massless field acquires an effective mass m 2 eff = λ/2(H I /2π) 2 N where N is an e-folding number [37,38]. Then such an interaction generates the vacuum energy proportional to ρ ∼ λ/8(H I /2π) 4 N 2 . For an infinite N , it approaches an equilibrium value 3H 4 I /16π 2 which is independent of λ. The ratio of this energy density to the critical density of the present universe is given by Ω φ 4 ∝ λN 2 H 4 I /(384π 4 M 2 Pl H 2 0 ). If H I = 10 −5 M Pl and N = 10 2 , it becomes Ω φ 4 ∼ λ × 10 100 . Hence unless we take λ ∼ 10 −100 , it exceeds the critical energy density of the universe. We want to come back to this issue in future.

A Time evolution of the zero-momentum mode
In this appendix, we study zero-momentum approximation to obtain the time evolution of the EMT at late times in the MD period, (4.8). We also reproduce the late time behavior in the RD period, (3.25) and (3.27) 9 .
The wave equation, ( + m 2 )u = 0, in the Robertson-Walker metric is written as u ′′ + 2Hu ′ + k 2 + (ma) 2 u = 0 , We first see the freezing behavior of the low-momentum wave function at early times. Since at early times,u ∼ Hu and m < H, we can neglect the m 2 term in (A.3). Then it is easily solved as where F and G are arbitrary constants. The first term with F is time independent and represents the frozen wave function. In order to connect to the Bunch-Davies wave function in the inflation period, which is also time independent, one has to choose F (k) = 0 and G = 0. Note that once these coefficients are fixed, the solution (A.4) continues to be valid throughout the history, both in the RD and in the MD periods, as far as the condition k phy , m ≪ H is satisfied. If the scale factor is a ∝ t p , the Hubble parameter is given by H = p/t and the solution to the equation where F ′′ and G ′′ are arbitrary constants. Note that the solutions (A.5) and (A.7) are good approximations to (A.2) for low-momentum modes with k phy ≪ m, √ mH, irrespective of early times with mt < 1 (i.e., m < H) or late times with mt > 1 (m > H).
As we studied in section 4, we are interested in a situation that the wave functions continue to be frozen until the MD period and then start oscillation. Then the coefficients F ′′ and G ′′ in (A.7) can be determined by requiring that it is constant near t ∼ 0 and the amplitude is given by the wave function in the RD period (2.37), or equivalently by that in the inflation period (2.13). Hence the wave function in the MD period is given by where s = sin(2mt) = sin 4m 3H , c = cos(2mt) = cos 4m 3H . (A.11) These expressions are valid both at the early times m < H (mt < 1) and at the late times m > H (mt > 1) since the wave function (A.7), and thus (A. 8), are valid at both times. Indeed, they also give the early time behavior (4.7) by expanding the trigonometric functions with respect to m/H. We may consider another scenario where the frozen behavior turns to the oscillating one in the RD period, as analyzed in section 3. In this case, we can connect the wave function (A.5) in the RD period, with p = 1/2 and ν = 1/4, directly to the wave function in the inflation period (2.13), and obtain u = iH I √ 2 k where we used the expansion formula (2.12) of the Bessel function near the origin. The late time behavior is obtained by using the asymptotic form of the Bessel function J ν (z) = 2 πz 1 + O(z −2 ) cos z − 2ν + 1 4 π + O(z −1 ) sin z − 2ν + 1 4 π (A.13) at |z| ≫ 1. This agrees with the WKB wave function at low momenta (3.24). Then we can obtain the late time behavior of the EMT in the RD period, (3.25) and (3.27).