A quantum heating as an alternative of reheating

To model a realistic situation for the beginning we consider massive real scalar $\phi^4$ theory in a (1+1)-dimensional asymptotically static Minkowski spacetime with an intermediate stage of expansion. To have an analytic headway we assume that scalars have a big mass. At past and future infinities of the background we have flat Minkowski regions which are joint by the inflationary expansion region. We use the tree-level Keldysh propagator in the theory in question to calculate the expectation value of the stress-energy tensor which is, thus, due to the excitations of the zero-point fluctuations. Then we show that even for large mass, if the de Sitter expansion stage is long enough, the quantum loop corrections to the expectation value of the stress-energy tensor are not negligible in comparison with the tree-level contribution. That is revealed itself via the excitation of the higher-point fluctuations of the exact modes: During the expansion stage a non-zero particle number density for the exact modes is generated. This density is not Plankian and serves as a quench which leads to a thermalization in the out Minkowski stage.


I. INTRODUCTION
Standard reheating [1] after inflation [2], [3], [4], [5] demands the presence of inflaton field. One of the important ingredients of this picture is the belief that during the rapid expansion stage the matter in the unverse can only cool down. This is based on the common wisdom that quantum effects provide only very tiny contributions to such quantities as the stress-energy fluxes.
In fact, the so called one loop expectation value of the stress-energy tensor, which is calculated with the use of the tree-level two-point Hadamard function or Keldysh propagator usually does not provide strong contributions to the particle number density after inflation. The rapid expansion always wins. The non-trivial value of the stress-energy tensor, if any in such circumstances, is due to the excitation of the zero-point fluctuations, which are the only ones contributing to the tree-level two-point function in the ground state. At the same time, it is believed that quantum loop corrections to the expectation value under discussion provide only UV renormalization to physical quantities. However, this observation is based on the intuition gained in stationary vacuum quantum field theory loop calculations, which in generic situations are not applicable in such non-stationary situations as rapid inflationary expansion.
There is already a vast literature showing that quantum loop corrections during inflation can grow in time and become comparable with tree-level contributions. See e.g. [6]- [57] for incomplete list of references. There are different sorts of secular effects. There are such effects which are specific only for the de Sitter space-time massless non-conformal scalars and gravitons (see e.g. [24]). The other secular effects are present in generic non-stationary situations and appear when the time separation between arguments of the two-point function is growing (see e.g. introduction of [56] for the review). The latter type of secular effects leads just to a mass renornalization or reveals a quasi-particle instability, if any. However, there are also such effects which may provide dramatic particle production [48]- [57] within comoving volume. They appear in the two-point functions when both of their arguments are taken to the future infinity, while the time separation between them is held fixed. Of course the intensity of the particle production depends on their mass and on their initial density. However, even for very big mass and for vanishing initial density it may be comparatively relevant, if the expansion stage is long enough, as we show in the present note. It goes without saying that for light fields there can be even more dramatic effects [56].
It is known in condensed matter theory that in non-stationary situations generally loop corrections can be strong [58], [59] for various reasons. Furthermore this fact has been observed in de Sitter space [57], [56] in strong electric field backgrounds [61], [62], in the case of loop corrections to the Hawking radiation [63] and in the case of moving mirrors [64], [65]. In all, even fundamental quantum fields in non-stationary situations reveal a similar behaviour to quantum fields in a medium, as in condensed matter theory.
The secular effect of interest for us in the present paper is due to the excitation of the higher levels (on top of the zero-point fluctuations) of the exact modes in background gravitational field. This is due to the non-stationarity of the inflationary expansion and due to the presence of interactions in the quantum field theory [57] (see also [56]). In non-Gaussian (selfinteracting) theories in non-stationary situations initial ground state of the quantum field theory is changing in time that may reveal itself in dramatic changes of correlation functions, which are not observed in the proper vacuum quantum field theory measurements and calculations. In particular, the secular effects under discussion lead to a growth of the physical particle number density in time and may result in its non-zero value at the final stage of the expansion. That can be true even for very massive fields, if the expansion stage is long enough. As a result in this note we point out that heating can be achieved even without inflaton field, for any physical origin of the cosmological expansion and even for very massive fields.
For simplicity in this note we consider a model situation of the 2D φ 4 quantum field theory in the gravitational background as follows. At past and future infinities of the background we have flat Minkowski regions which are joint by a long inflationary expansion stage.
The paper is organized as follows. In the section 2 we define the gravitational background and describe the properties of the free mode functions. In the section 3 we perform the standard calculation of the stress-energy expectation value with the use of the treelevel two-point Hadamard or Keldysh function. In the section 4 we calculate lowest loop corrections to the correlation functions and show that some of them grow with time. In the section 5 we perform the resummation of the leading corrections from all loops during the expansion stage and find the created particle density at the moment of exit from the expansion. In the section 6 we derive the kinetic equation in the final Minkowski stage and show that the created particle density serves as a quench for the thermalization. The section 7 contains conclusions.

II. SETUP OF THE PROBLEM
The action of the theory that we consider in this paper is where the metric is that of (1+1)-dimensional Robertson-Walker spacetime: where C(η) is the so called Conformal factor.
In particular we consider a (1+1)-dimensional asymptotically static universe with Minkowskian in-and out-regions and with a de Sitter phase of expansion, with the following metric: This metric describes three phases (see fig. 1): 1. First -in-Minkowski stage

de Sitter expansion stage
3. Second -out-Minkowski stage We assume that there is the infinite flat region with the metric which is glued to the last stage. Physically we can think that at some moment of time a cosmological constant was somehow created resulting in a de Sitter expansion that stops in a while, leaving the space flat again. The origin of this cosmological constant is irrelevant for the discussion in the present paper. For the de Sitter expansion stage the scale factor has the exponential form a(t) ≈ e 2Ht .
In our case the Hubble expansion rate is H ≈ 1 T .

A. Modes in the expansion stage
During the expansion stage the equations of motion for the free modes in the theory under consideration (λ = 0) are as follows: One can represent the modes which solve this equation as φ k (η, x) = g k (η)e ∓ikx and the equation of motion for their temporal part, g k (η), is then If one assumes the ansatz g k (η) = |η| 1/2 h(kη), then we obtain the Bessel equation for h(kη) with the index ν = iµ = i m 2 T 2 − 1 4 . Hence, one can write that [60]: where H (1) and H (2) are Hankel functions of the first and second kind and C 1 ,C 2 are some complex coefficients that we will fix by gluing solutions at η = −T and by the normalization conditions.
Hankel functions with such an index behave as follows: for some complex constants A ± .
If the field is heavy, m > 1/2T , in de Sitter space it is said to belong to the principal series, while the light field, m < 1/2T , belongs to the complementary series. Modes of the principal series oscillate and decay to zero as η 1/2±iµ when η → 0, while for the light fields from the complementary series modes homogeneously decay to zero as η 1/2± when η → 0. In the present paper we work under the assumption of the heavy fields, m 1/2T , so that ν ≡ iµ is purely imaginary. It is understood in this case that we are talking about the analogue of principal and complementary series of de Sitter space as they are representations of the de Sitter isometry group that is not present in our case.
So for this reason from now on we will just refer to heavy and light fields.
The quantum field can be expanded in the usual way: Annihilation, a k , and creation, a † k , operators obey the standard commutation relations as a corollary of the fact that g k (η) obeys the Klein-Gordon equation.
Because of the metric that we are considering here, the free Hamiltonian depends on time. Hence, the system under consideration is in a non-stationary state and one has to apply the Schwinger-Keldysh (aka in-in, aka Closed Time Path) diagrammatic technique.
An introduction to it can be found in [58], [59]. Every particle in this formalism is described by the matrix propagator: for which the entries are the Keldysh D K , Retarded D R and Advanced D A propagators: In this paper we consider only spatially homogeneous states and due to the spatial homogeneity of the background it is convenient to perform the Fourier transformation of all quantities along the spatial direction: Here n(k) ≡ a † k a k , κ(k) ≡ a k a −k and κ * (k) ≡ a + k a + −k are the number density and anomalous quantum average for the exact modes with respect to the state under consideration. In the free field theory n(k) and κ(k) do not change in time and, in particular, are vanishing if one starts in the vacuum state, a k |in = 0, for all k. However, if one turns on interactions, then generally these quantities reappear and depend on time [57]. Furthermore, if the anomalous quantum average vanishes, then n(k) g k g * k acquires the meaning of the particle number density per comoving volume, given that in this case the free Hamiltonian of the theory is diagonal. This will be important for the physical interpretation of the kinetic equation approach that we will present later.

B. Gluing of the modes
The approximate behaviour of the modes in the entire spacetime under consideration is as follows: where ω in (k) = √ k 2 + m 2 and ω out (k) = k 2 + m 2 T 2 2 . This can be seen from the approximate form of the Klein-Gordon equation in the corresponding parts of the space-time.
To estimate the coefficients A, B, C, D we need to glue the solutions and their first derivatives at η = −T and η = − . We start with the gluing at η = −T : For the large momenta, when k|η| µ for all η ∈ [−T, − ] (note that in this case both kT µ and k µ) we have that: and from this system of equations we find the coefficients: where the index b indicates that the coefficients are valid in the large momenta approximation.
Since we are interested in the very heavy case mT 1, for high momentum kT µ = (mT ) 2 − 1 4 ≈ mT so that k m and hence ω in (k) ≈ k. So we have For low momenta, when k|η| µ for all η ∈ [−T, − ] (note that in this case both kT µ and kT µ), we have at η ≈ −T that: Solving this system of equations, we find the coefficients For later convenience, let us write the coefficients as follows: where Since we are working in the approximation that mT 1 for small momentum kT µ = m 2 T 2 − 1 4 ≈ mT and this leads to k m so that in the leading approximation For intermediate momenta, when k µ kT , we can distinguish between two regions: k|η| > µ and k|η| < µ. Hence, at η = − , the modes behave as H iµ ∝ (k|η|) ±iµ , while . For the intermediate momenta in the gluing at We continue with the gluing of the modes at η = − : We can then calculate the coefficients C k and D k from eq. (14) in complete analogy to the previous case. On this side, for large momenta, k|η| µ for all η ∈ [−T, − ], we have that: From this system of equations we find the coefficients: When mT 1 and k , kT µ we can use the approximations k µ ≈ mT , hence, k mT so that ω out (k) = k 2 + mT ≈ k as well as ω in (k) ≈ k. In this case we can For low momenta, k|η| µ for all η ∈ [−T, − ], we have that: Solving the system we find the coefficients When mT 1 and kT, k µ we can use the approximations k µ ≈ mT , hence, k mT so that ω out (k) = k 2 + mT mT as well as ω in (k) ≈ m. In this case we can approximate the coefficients C s k and D s k as Finally, in the intermediate region of momenta, for the η = − gluing, we find the following relations: which in the limit we are working with, may be further approximated as follows: We will need to perform the gluing of the modes with the use of the de Sitter out Jost modes in the section 6 to look at the out Minkowski region. The reason for that will be explained below.

III. TREE-LEVEL EXPECTATION VALUE OF THE STRESS-ENERGY TENSOR
We now evaluate the tree-level flux as this is usually associated with particle creation in such a situation as we discuss here. The expectation values of the stress-energy tensor components can be calculated with the use of the well known relation where T µν [u k , u * k ] denotes the bilinear expression for the energy momentum tensor in terms of the mode functions u k [66]. To calculate the particle flux we need to use the in-modes in the out region: because we are interested in the value of in| T µν |in as t → +∞, where |in is the ground state with respect to the in-modes. The latter we consider as the initial state of the problem under consideration.
We can then calculate the integral defining in| T µν |in with the use of the asymptotic expressions for the coefficients C k and D k that we have found in the previous section.
To do that we separate the integration into three regions of low, intermediate and high momenta: where the regions of integration, Ω i , are defined as follows The flux is then given by the expectation value of the non-diagonal components of the stress-energy tensor The last equality follows from the fact that the integrands are odd functions in symmetric integration intervals. In principle one could obtain non-zero fluxes from the momenta For the other components of the stress-energy tensor we find: In the last two expressions we obtain the standard UV divergences as in flat space-time, which are coming from the Ω 3 region, as it should be. By the normal ordering we can cancel these contributions. The rest are the expectation values that we are looking for. In the next section we will calculate loop contributions to the two-point correlation functions, which will substantially correct the expressions found in this section.

IV. TWO-LOOP CONTRIBUTIONS IN DE SITTER EXPANSION STAGE
In this section we discuss the leading infrared loop contributions to the retarded, advanced, Keldysh propagators and vertexes in the limit as η 1,2 → − and T → +∞. See, for example, [67] for the Schwinger-Keldysh diagrammatic technique in φ 4 scalar field theory in cosmology.
It is straightforward to see that one loop bubble diagram corrections to the two-point functions lead just to a mass renormalisation (perhaps time dependent). They generally do not grow as T / → ∞. Also it is straightforward to see that retarded and advanced propagators do not receive growing with T / corrections at two-loop sunset diagram order [57] (see also [59] for a related general discussion). Furthermore, it can be shown that for massive fields vertexes also do not receive growing corrections [57]. (Note that this is not true for very light fields [56].) Thus, we continue with the two-loop sunset diagrams for the Keldysh propagator.
Leading contributions to the Keldysh propagator at 2-loop order D K 2 (k|η 1 , η 2 ), in the limit T → +∞, are contained in the following expression [57]: with where F (k|η 3 , η 4 ) = dq 1 2π and g k (η) is defined in eq. (14), while C(η) is in eq. (3). The subscript 2 of D K , n, κ and κ * indicates that we are discussing here only the second loop corrections to these quantities.
The above expressions appear from the following 2-loop correction to the Keldysh propagator: if we neglect the difference between η 1 , η 2 and η = √ η 1 η 2 → − and set the initial time to be −T . Such a neglection can be done if we keep only the largest contribution to n, κ and κ * as T → +∞ [57].
The Schwinger-Keldysh diagrams defining eq. (39) are as follows: if subleading terms are neglected. Using this expression we can find from eq. (38) that n 2 is as follows where q 3 = |k − q 1 − q 2 |. The largest contribution to this expression comes from the region where q 1,2,3 k [57]. (Note that this is true only for the massive fields [56].) To evaluate the last integral we make the following change of the integration variables: Then, it becomes equal to The integration over the internal momenta, q 1,2 , can be separated into four regions: The first one is when |q 1,2 | µ √ η 3 η 4 , the second region is when |q 1,2 | µ √ η 3 η 4 , the third region is when |q 1 | µ √ η 3 η 4 and |q 2 | µ √ η 3 η 4 and, finally, the fourth region is when |q 1 | µ √ η 3 η 4 and |q 2 | µ √ η 3 η 4 . As a result, For low internal momenta part of the integral, |q 1,2 | µ √ η 3 η 4 , we have that: and the corresponding contribution to n 2 is as follows: where σ 1 (mT ) denotes the integral factor as a function of T and the logarithmic behaviour follows from the fact that the integrand of du/u does not depend on u. This is the consequence of the fact that we have an approximate scaling symmetry, η → aη and p → p/a, in the expansion stage, which is just a part of the de Sitter isometry group.
For large internal momenta part of the integral, |q 1,2 | µ √ η 3 η 4 , we have that: because of the approximate form of the coefficient A b q , B b q in eq. (18). The corresponding contribution to n 2 is as follows: Hence, the region of high internal momenta also provides the logarithmic growth.
In the other two regions of integration ( we have a very similar situation, with a resulting growth of the form λ 2 m 4 log T σ 3,4 (mT ), with some σ 3 = σ 4 . Thus, for low external momenta, k|η| µ, there is a logarithmic growth of n 2 (k) over the entire region of internal momenta: where σ small (mT ) = 4 i=1 σ i (mT ) and in σ i (mT ) we include T-dependent prefactors. For the case of κ 2 (k) we have the same situation. The only difference is in the Tdependent prefactors. In this case we have that: As a result: Hence, with the same manipulations as for n 2 we can just state that with a factor Γ small (mT ) which follows from the previous equation.

Now let us look at the intermediate external momentum region, k µ and kT µ.
In this case we have that: Hence, and as a result, Similarly to the above discussion one can show that the product under this integral depends only on v. Hence, the calculations are essentially the same as we did for small external momenta. The only difference is in the argument of the logarithm and in the T-dependent prefactors. Hence, Similarly for κ 2 (k) we use that: and, as a result, Here again as before, the product under the integral depends only on v. Hence, even for intermediate momenta we have logarithmically growing corrections both for n and κ, κ * .
At the same time, on general physical grounds it should be obvious that for high external momenta k|η| µ we do not have any secular growth because the corresponding modes are not sensitive to the curvature of the space-time and behave as if they are in flat space. It is straightforward to show the latter fact explicitly (the situation is similar to that in flat space [57]). Hence, below we concentrate on the small and intermediate parts of the external momenta, k, regions.

V. RESUMMATION OF LEADING CONTRIBUTIONS FROM ALL LOOPS IN THE EXPANSION STAGE
We see that even if λ 2 is small λ 2 log ( /T ) and λ 2 log (k /µ) ≈ λ 2 log (k /mT ) can become large as T / → +∞, hence quantum loop corrections are not suppressed in comparison with classical tree level contributions to the propagators. That evidently may strongly affect particle production. This means that we need to sum at least the leading corrections from all loops, and for that we have to use the system of Dyson-Schwinger (DS) equations [57]. That is the mathematical part of the story. As a result, we have to deal only with the single equation for the Keldysh propagator.
Furthermore, because we are working here with the heavy fields, mT 1, the modes are oscillatory functions for both small and large momenta, i.e. they oscillate even when the system exits from the expansion stage. This results in the standard separation of scales between the time dependence of the modes g k (η) and that of particle density n k (η) and anomalous quantum average κ k (η) and κ * k (η) [58], [59]. Namely n, κ and κ * can be considered as very slow functions of time in comparison with g k (η). This fact also allows to simplify the DS equation for the Keldysh propagator. (Note that this is not the case for the light fields [56].) Thus, we assume that the quantum system under consideration had started its evolution at the vacuum state at the beginning of the de Sitter expansion stage, i.e. n and κ were remaining zero before η = −T . Hence, the tree-level or initial value of the Keldysh propagator D K 0 is given by eq. (12) with n k = 0, κ k = 0 at η = −T . Our goal is to find the values of D K and of n k , κ k at the exit from the expansion stage, i.e. at η = − . Then, the relevant part of the system of DS equation takes the following form [57]: where D R,A 0 are the initial (tree-level) values of the retarded and advanced propagators and D K is the exact Keldysh propagator whose form will be specified in a moment.
This equation is covariant under Bogoliubov rotations between different modes. For all the family of the modes which are related to each other via Bogolubov rotations [69], [68], the following ansatz solves this equation [57]: Compare it to the eq. (12).
Above to simplify equations we have adopted a bit different notations. Namely we denote the exact n and κ, κ * taken at a time η as n k (η) and κ k (η). Note also that due to the approximate scale invariance, which is present in the expansion stage as a part of the de Sitter isometry (and which is violated only in the normalization factors of the modes), we can assume that n k (η) ≈ n(kη) ≡ n kη [57]. up the initial quench for the reheating stage in the future Minkowski region, which will be discussed in the next section. At the same time, the generation of non-zero anomalous averages κ and κ * just means that the initial ground state of the theory at past infinity is not anymore its ground state at future infinity.
Furthermore, the fact that the DS equation is covariant under Bogoliubov rotations between different modes and that for out-modes this equation has the solution with vanishing κ and κ * is an evidence that the proper ground state of the theory is the out-vacuum [57].
All this means that to find a solution of the DS equation we have to perform the Bogoliubov rotation from in-to the out-modes within the equation (50), plague there the ansatz (51) with κ = κ * = 0, assume that n is a slow function in comparison with g k , take correspondingly g k to be the out-modes and, finally, pick up on the right hand side the leading contribution in the IR limit in question [57]. The latter extraction is similar to that which have been done above for the second-loop correction. As the result of these manipulations we find the integro-differential equation: × (1 + n kη )(1 + n l 1 )(1 + n l 2 )(1 + n |l 1 +l 2 | ) − n kη n l 1 n l 2 n |l 1 +l 2 | , where l i ≡ q i η and we neglected k in comparison with q i , i = 1, 2, 3, as that is the region where from the largest contributions are following [57]. That is similar to the situation with the two-loop corrections.
For the low external momenta the out Jost modes behave as This means that the function φ in eq. (52) has the following form: The obtained eq. (52) is valid for both low and intermediate external momenta.
First thing to observe is that the kinetic type equation (52) does not possess Planckian distribution as its solution, because there is no energy conservation in time dependent backgrounds.
If we assume that the initial value of n after the rotation to out Jost modes is very small 1 , we can use the approximations adopted in [57] to find the stationary solution of eq. (52). If n 1 the latter equation reduces to [57]: where Γ 1 and Γ 2 are decay and production rates respectively, given by Hence, for all |k| < µ/ we have flat stationary distribution: 1 It should be much smaller than one, which is the case for the vanishing initial values of n, κ and κ * , if the mass m is very big. Note that the initial values of n, κ and κ * are taken to be vanishing for in-modes, but after the Bogoliubov rotation to the out-modes neither initial value of n nor that of κ and κ * does vanish any more. However, for large mT they are much smaller than one. In such circumstances it was shown in [56] that κ and κ * evolve to zero and n solves the eq (52).
We can estimate the order of magnitude of this ratio if µ ≈ mT 1 to verify the consistency of the result using the fact that the integrals are saturated at l ∼ µ, hence, neglecting the contributions from the remaining integration region [57]. It is easy to see then that the so approximated integrand functions are just limits of products of Bessel functions, so we have the following estimate, in the limit T → +∞: This is the initial value for the dynamics in the future Minkowski region, which we discuss in the next section.

VI. KINETIC EQUATION IN THE MINKOWSKI OUT REGION
Thus, at the exit from the expansion stage for all such momenta that |k| < µ/ we have the flat particle number density (56), which is small, but non-zero. It is parametricaly small because we are dealing with large mass mT 1. This distribution serves as an initial quench for the dynamics in the outer Minkowski region. We expect an eventual thermalisation, because the latter region is flat and particle kinetics there is the same as in flat space. But we would like to see it from the first principles. That is what we will do in this section.
Above to calculate the tree-level expectation value of the stress-energy tensor we have glued the in-modes across η = − to find their behavior in the out Minkowski region.
However, in the previous section we have observed that due to self-interactions the theory under consideration relaxes to the out-vacuum. Hence at the exit from the expansion stage the single quasi-particle states are represented by out-modes. As a result to proceed we have to glue the de Sitter out-modes rather than in-modes across η = − . Namely, to find the form of the modes in the out Minkowskian region we have to glue the solution in the out flat region with the out Jost modes from the de Sitter stage.
For small momenta, k|η| µ, the gluing conditions give the following system of equa- Solving this system we find the coefficients and with the usual approximations ω in ≈ m, ω out ≈ mT and µ ≈ mT we have the following so the modes in the out region behave as single waves where we also used the explicit form of the coefficient of the Jost function, F s k , calculated in the previous section in eq. (53).
Let us now apply the same gluing conditions for the modes of intermediate momenta As a result We have seen that the particle density created during the de Sitter expansion stage is not Planckian (57). Hence, our purpose now is to find a kinetic equation for n k (t) in the out Minkowski region using as initial value for n the one from (57). The point is that if the initial density is not Plankian, i.e. is not stationary for flat space kinetic equation, it should evolve in time towards the stationary Plankian distribution. To see that, essentially we have to perform similar calculations to those we already did to find the kinetic equation during the de Sitter expansion stage, but using this time the modes we found in this section above. Note that in general anomalous average κ k also grows along with n k , but because the modes under consideration represent single waves, it is not hard to see that κ k ≈ 0 [57].
Let us write the modes in this out region in the following way: As we will see this form will be good for both regions of momenta, considering |Q T, | 2 = 1 in both small and intermediate momenta case. From now on, for the sake of notational simplicity we will write just ω having in mind we are talking about ω out .
Computing the cosine time integrals we obtain terms of the form sin[∆ω(t+ )] ∆ω and in the limit of t → +∞ these are reduced to δ-functions ensuring energy conservation, as it should be the case in flat space. The only allowed process is then the scattering one between scalar particles, so the collision integral of the kinetic equation just contain the term with (1 + n k )n q 1 n q 2 (1 + n |k−q 1 −q 2 | ) − n k (1 + n q 1 )(1 + n q 2 )n |k−q 1 −q 2 | . Given that our initial density is small, n k 1, we can approximate the expression as n q 1 n q 2 −n k n |k−q 1 −q 2 | .
After these manipulations and approximations we are left with the standard Boltzmann equation and we can immediately see that the stationary solution is given by the equilibrium Boltzmann distribution n k ∝ e − ω(k) τ for a constant τ , resulting in thermalization. The value of this constant can be estimated from the total energy that was stored by the particles of the density (57) at the moment of exit from the expansion stage. I.e. we should have an equality µ/ −µ/ n ω out (k)dk ≈ τ , hence, τ ∼ mT e −3mT .

VII. CONCLUSIONS
In the section 3 we have calculated the stress-energy expectation value due to the zero point fluctuations. That is due to 1/2 term in eq. (12), i.e. when n, κ and κ * are vanishing. In the section 4 we have shown that quantum loop corrections to n, κ and κ * are not negligible, if the expansion stage is long enough. This signals the breakdown of the perturbation theory or semi-classical approximation and calls for the resummation of at least the leading corrections from all loops. That is done in the section 5. The generation of non-zero n, κ and κ * affects the picture observed in the section 3. At the exit from the expansion stage we have a non-zero physical particle number density even for very massive fields. The latter one is thermalized in the future Minkowski region, as is shown in the section 6.
In this note we have considered very massive fields to have an analytical headway.
For the massive fields the backreaction on the background gravitational field is negligible.
That is because the comoving density (56) remains finite, which means that the physical one gets diluted as the spatial volume increases during the expansion. As a result, the expectation value of the stress-energy tensor due to created particles has a negligible effect on the background geometry.
However, in [57] it is shown that even for massive fields the kinetic equation (52) has explosive solution with a dense enough initial state. The problem with such a situation is that for the latter value of the initial particle number density the matter energy density is comparable or even grater than the cosmological constant and, hence, the discovery of the explosive solution is not consistent.
However, for the light particles the situation is quite different, as is shown in [56]. In the latter case the explosive solution of the corresponding substitution of the kinetic equation appears even when the initial matter energy density is much smaller than the cosmological constant. Moreover, light particles are much more interesting from the phenomenological point of view. Hence, we plan to discuss light particles elsewhere. But from the computational point of view the consideration of the light particles is much more complicated than that of the massive particles. Apart from all it demands a consideration of the coupled system of the gravitational field and of the matter, because in such a case the backreaction on the gravitational background can be non-negligible [56].
We would like to acknowledge discussions with S.Alexeev and L.Astrakhantsev. The work of ETA was done under the partial support of the RFBR grant 16-02-01021 and of the state grant Goszadanie 3.9904.2017/8.9.