{'fraction_non_alphanumeric': 0.03898763463552658, 'fraction_numerical': 0.016070284043833732, 'mean_word_length': 4.016074649102172, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 37, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Reheating and thermalization in a simple scalar model

We consider a simple model for the Universe reheating, which consists of a single self--interacting scalar field in Minkowskian space--time. Making use of the existence of an additional small parameter proportional to the amplitude of the initial spatially homogeneous field oscillations, we show that the behavior of the field can be found reliably. We describe the evolution of the system from the homogeneous oscillations to the moment when thermalization is completed. We compare our results with the Hartree--Fock approximation and argue that some properties found for this model may be the common features of realistic theories.


I. INTRODUCTION
The problem of decay of homogeneous oscillations of the inflaton field after inflation is one of the most important issues of inflationary cosmology. This process is responsible for the reheating of the Universe, in particular, the reheating temperature depends on the details of the decay of homogeneous oscillations. Recently, this problem has attracted considerable interest, and new investigations have revealed a rather complicated dynamics of the inflaton field. In the old approach, which is called "elementary theory" [1], the decay of the spatially homogeneous oscillations of the inflaton field is described by a friction term in the equation of motion, which has a coefficient proportional to the decay rate of the inflaton quanta. This picture leads to a slow, exponential decay of the amplitude of oscillations. It is not until recently that it has been discovered that in most important cases this picture is not correct.
The reason for that is usually related to the existence of the resonance bands in the spectrum of the scalars interacting with the inflaton field, which contains modes that are unstable on the background of homogeneous oscillations of inflaton field. The importance of parametric resonance was first realized by Kofman, Linde and Starobinsky [2]. The modes affected by this phenomenon may be of the inflaton field itself or belongs to any scalar coupled to the inflaton.
As the importance of the parametric resonance has been understood, a new scenario has be proposed [2] (see also [3]). According to the latter, the reheating process occurs in two steps. At the first step, dubbed "preheating", the amplitudes of the quantum modes inside the resonance bands increase exponentially, giving rise at some time moment to the creation of a large number of quanta of inflaton or other scalar fields. Due to the exponential character of the growth of the resonance modes, the creation of the particles by parametric resonance is an explosion-type process, where most particles are emitted at the very last moments of the preheating epoch. After this explosion, most energy is carried by the created particles, and the old theory of reheating can be applied to the thermalization of the products of explosion and to the decay of the remnant of the homogeneous oscillations.
Despite the fact that there is a growing number of works aimed at extracting physical consequences from the scenario advocated in Ref. [2] (see Refs. [4][5][6]), the latter deserves more detail theoretical elaboration in order to take into account important physical effects that may change many details of the evolution of the Universe after the explosion. Let us mention only two of these effects. The first is the phase coherence between the modes created by parametric resonance: the wave functions of particles created in this way have almost the same phase, which is correlated with the phase of the homogeneous oscillations.
This phase correlation, presumably, may have certain consequences, in the same manner as the coherent nature of the homogeneous oscillations lead to the parametric resonance. The second important effect is the scattering between the particles, which cannot be neglected even in theories with small coupling. In fact, according to the standard picture, the particles created by preheating have very large occupation numbers, typically of the order of the inversed coupling. These large occupation numbers lead to strong Bose enhancement of the scattering processes, and the system is effectively strongly coupled. This implies, in particular, that after the explosion, the relaxation time is of the same order as the period of oscillation, and the effect of particle collision cannot be neglected. Existing calculations [2,7] mostly rely on the Hartree-Fock approximations, where the effect of particle collision is not taken into account.
Closely related to the decay of homogeneous oscillation is the thermalization among products of parametric resonance particle production. It has been argued that since the particles produced by parametric resonance has energy much smaller than the equilibrium temperature, the thermalization takes a long time interval as the particles have to go through may collisions in order to overcome the wide separation of energy scales. However, no theory of thermalization has been developed, as well as no reliable estimation of the thermalization time has been made so far.
That the Universe after the preheating explosion is effectively strongly-coupled as mentioned above makes the study of the dynamics of the Universe quite difficult, and it has been suggested that numerical methods may be used to simulate processes occurring after inflation [8]. Numerical simulations of this type has been carried out for the simplest φ 4 theory in Ref. [9], but in order to understand the qualitative features in the evolution of the Universe after inflation, as well as to check the validity of the Hartree-Fock approximation, it is desirable to have a simple model where this evolution can be followed in the whole time interval from the beginning up to the thermalization in a reliable way.
In this paper we consider one of such models. We neglect the expansion of the Universe, and consider the theory of one self-coupled scalar filed without symmetry breaking. The initial amplitude of the homogeneous oscillations is taken to be small enough for the latter to be almost sinusoidal. More precisely, the amplitude is supposed to be A 0 m/ √ λ, where m is the boson mass and λ is the coupling, and A 0 is some small dimensionless parameter, A 0 ≪ 1.
Despite the fact that the effect of parametric resonance is suppressed by a power of A 0 , we will show that the analogue of the elementary theory of reheating is totally inapplicable when A 0 ≫ √ λ. In particular, this theory is invalid in the regime that we consider throughout the paper, where A 0 is small, but fixed, and λ → 0. By considering the regime where the elementary theory breaks down, we let our model possess the main feature of the theories where parametric resonance is important. That the evolution of the Universe can be found reliably in our regime is the consequence of the existence of an additional small parameter A 0 in our model.
Let us briefly outline our results. At first, when the field makes homogeneous oscillations, the modes inside the resonance bands are enhanced and their amplitudes increase exponentially with time. The modes that are enhanced most are those around the center of the resonance band near 2m. At some time moment, when the energy carried by the resonance modes is still small, the back-reaction of these modes to the homogeneous background and their self-interaction become important and the exponential growth of the modes terminates. The dynamics of the field becomes rather complicated, but during some time interval it can be treated by the Hartree-Fock approximation. We see that the coherent nature of the modes created by parametric resonance leads to a period of time where some portion of the energy is transfered back and forth between the homogeneous oscillations and the resonance modes. Therefore, the decay rate is reduced by the coherence effect. After some time interval the coherence is destroyed, and the system goes to another regime where created modes are uncorrelated. In this regime, there is still coherent enhancement of modes by parametric resonance, but now the width of the resonance band is much smaller than that of the particle spectrum. Moreover, the resonance band moves towards the region of smaller energies, so each mode with a given momentum is enhanced during only a finite time interval, after than they become uncorrelated and can be considered as ordinary particles.
By this mechanism, the long-run effect of parametric resonance is non-explosive production of particles with energy approximately 2m. The rate of particle production slowly changes with time, but not by many orders of magnitude.
After a period of particle emission, when most energy is still carried by the homogeneous mode, the density of created particles is already enough for the non-coherent scattering between the particles becomes important and the Hartree-Fock approximation breaks down.
To describe the system in this regime, we make use of the Boltzmann kinetic equation. We find that the production of particles via parametric resonance occurs concurrently with their dissipation to other modes by scattering, which leads to the excitation of modes with larger and larger energies. It is very interesting that the dynamics of the scalar field is similar to that of turbulent systems, where energy is transfered from long-distance to short-distance modes. In particular, the Kolmogorov power spectrum is formed after some time interval among the particles with nonzero momenta.
During the whole period describe above, the amplitude of homogeneous oscillations gradually decreases. When the energy carried by the homogeneous mode becomes much smaller than the total energy, the thermalization epoch begins. The thermalization can be divided into two steps where the distribution function shows two different self-similar behaviors, depending on whether the Bose condensation is effective compared to condensate evapora-tion. Let us only note that the typical energy of particles increases continuously, while the typical occupation number decreases. The time required to achieve thermal equilibrium is estimated to be t ∼ A 1/6 0 λ −7/4 m −1 . Therefore, the results of our paper show that the dynamics of the scalar field after the explosion is much more complicated than predicted by the Hartree-Fock approximation, and there is interplay between various processes, including particle creation, scattering, and Bose condensation, which must be accurately taken into account in order to find out the actual behavior of the system. Generally, we find that the scattering slows down the decay of the homogeneous oscillations. Our investigation suggests that features discovered in our model may be common in some realistic theories.
The paper is organized as follows. In Sect.2 we describe the model, and show that the elementary theory of reheating is not valid in a particular regime in this model. We consider the parametric resonance, in particular the strongest resonance band that is located at energies closed to 2m. We show that the linear regime, where back-reaction of the created particle to the background can be neglected, breaks down when the energy carried by produced particles is still small. Sect.3 is devoted to the Hartree-Fock approximation, where particles are considered as moving on a self-consistent background. We derived the Hartree-Fock equation for the amplitude of the wave functions of the particles, and show that the phase coherence of the latter slow down the decay of the homogeneous oscillations. We find analytically the rate of particle creation in the long run. In Sect.4 we take into account the scattering among produced particles and describe the evolution of the system in the regime where both particle creation and scattering occurs at the same speed. We show that the part of our system that contains particles with nonzero momenta behaves like a typical turbulent system. In particular, we find the Kolmogorov index of turbulence for our case. In Sect.5 we discuss the behavior of the system after the amplitude of homogeneous oscillations become small is considered. The thermalization is described, and the time required for it is also estimated. Finally, Sect.6 contains concluding remarks.

II. THE MODEL
Let us consider the theory of one scalar field in flat Minkowskian space-time, where the coupling constant λ is assumed to be small, λ ≪ 1. For convenience we take the mass of the boson to be 1. We suppose that at the beginning the field φ makes spatially homogeneous oscillations. This can be done by adding a source term Jφ, where J is a constant at t < 0, to the Lagrangian (1) and turn off J at t = 0. We assume that the amplitude of the oscillation is much smaller than 1/ √ λ, so that it is almost sinusoidal. To the leading approximation, the time dependence φ is of the form where A 0 ≪ 1. We will be interested in the process of energy transfer from the homogeneous mode to the modes with nonzero momentum. The regime that we are interested in is small, but fixed, A 0 , and λ → 0. Let us first show that in this regime, the elementary theory of reheating is not applicable.

A. Elementary theory
The homogeneous oscillation (2) can be considered as a condensate of particles with zero momentum. The particle density in the condensate that corresponds to the oscillations (2) is In the elementary theory of reheating, the decay of the homogeneous mode (which will also be called "condensate evaporation") is viewed as a perturbative process which usually occurs via the decay of inflaton quanta in the condensate to other particles. In our case where the theory contains only one scalar, the φ-quanta are stable, but there are still processes that creates φ-quanta with nonzero momentum from those in condensate. The first perturbative process of this type is the 4 → 2 process, in which 4 φ-bosons at rest become to 2 particles with energy E = 2 each. The amplitude of this 4 → 2 process is not vanishing and proportional to λ 2 . The rate of the 4 → 2 can be obtained by multiplying the squared amplitude to the fourth power of n 0 (since it involves 4 initial particles), and the elementary theory in our case predicts The time scale of the decay of the condensate is However this naive calculation is incorrect. Though usually the breakdown of the elementary theory is related to the phenomenon of parametric resonance, it can be easily understood also in the perturbative framework. To do that, let us recall that in our calculations we assume implicitly that the two final particles of the 4 → 2 process are emitted as in vacuum.
However, the rate of this process is enhanced if there are already particles with energy E = 2. In fact, provided the occupation number at the energy level E = 2 is n, then the 4 → 2 probability is enhanced by the Bose factor of (1 + n) 2 . There is also the time-reversal transition where two particles with energy E = 2 go back to 4 particles at rest, which has a relative rate of n 2 . If one combines the two processes one finds that the overall rate of transition from the condensate is (1 + n) 2 − n 2 = 1 + 2n times larger than predicted by the elementary theory, which implies that the latter is correct only when the occupation number n is much smaller than 1.
Let us estimate the occupation number at the energy level E = 2. At the beginning this occupation number is 0, but it increases due to the radiation by 4 → 2 transitions.
The particles with energy E = 2 would be accumulated in large number if they remain untouched at the same energy level after they have been created, however in reality these particles scatter on those in the condensate and have a finite mean free time of order (where σ is the typical cross section in particle collisions, σ ∼ λ 2 ). Note that at the values of A 0 we are considering, A 0 ≪ 1, the mean free time t 0 is much smaller than the time scale in eq.(5). Each particle created from the condensate lives at the energy level E = 2 for a time interval t 0 , therefore the maximal density of particles with energy E = 2 that can be achieved is Though this density is much smaller than n 0 , the occupation number would be infinite if all these particles have exactly the same energy. However, due to the Heisenberg principle, the energy of particles has an uncertainty of order t −1 0 since they have a finite lifetime.
The particle distribution in the momentum space, so, can be imagined as a spherical shell around E = 2 with the thickness of t −1 0 . The occupation number can be found by dividing the particle density (6) to the phase space that the particles occupy, Recall that the condition for the elementary theory of reheating to be valid is n ≪ 1, one finds that the latter can be applied only when A 0 ≪ √ λ, or, in other words, the amplitude of the homogeneous oscillations must be much smaller than 1. Since in the regime considered in this paper (A 0 = fixed, λ → 0) this condition is violated, the elementary theory is not appropriate for our purposes.

B. Parametric resonance
Since the elementary theory of reheating breaks down, let us turn to the parametric resonance. For further calculations, we need a more accurate expression for the solution to the field equation than eq.(2). At small A 0 , the following formula for the solution can be derived, where ω = 1+ is the frequency of the oscillations. As noted above, the oscillations is almost a cosine: the higher harmonics come with smaller coefficients than the leading one.
If the system under consideration is classical, and the spatial homogeneity is exact, the oscillations (7) would continue infinitely long. However, this is not true in the quantum case.
The reason for the behavior of the quantum system to be different is that the homogeneous oscillation (7) is unstable: some vacuum quantum fluctuations become larger with time. To see this instability let us quantize the scalar field on the background (7). For this end we decompose the file operator into classical and quantum parts, φ = φ 0 +φ. The part of the Lagrangian that is quadratic onφ has the form, The quantization is completely analogous to that of the free scalar field on the vacuum background. As the result, one representsφ via the creation and annihilation operators where [a k , a † k ′ ] = (2π) 3 2ω k δ(k−k ′ ). If there is no homogeneous oscillations, one has f (k, t) = e −iω k t , where ω k = √ k 2 + 1, and eq.(9) reduces to the usual formula for quantized scalar field on the vacuum background. On the background of φ 0 , f (k, t) is a solution to the mode equation To find the initial condition for eq.(10), let us note that at t < 0 the source term Jφ is present in the Lagrangian, so φ 0 is constant, and the Lagrangian (8) is just that of a scalar field on the vacuum background. Neglecting the small correction to the mass, one has We will call f (k, t) the wave function of particles with momentum k (in fact it gives the time dependent of the wave function). Eq.(10) has the same form as the equation of motion of an oscillator whose frequency is a periodic function of time, ω 2 (t) = ω 2 k + λ 2 φ 2 0 (t). Such systems are well known in classical mechanics to exhibits the phenomenon of parametric resonance, which means that there are values of k where any solution to f (k, t) to eq.(11) becomes arbitrarily large in at least one of the limits t → −∞ or t → +∞ [10]. In the case when the dependence of the frequency on time is weak, as in our problem (A 0 ≪ 1), the solution to eq.(10) can be found analytically, and the parametric resonance occurs when ω k is close to 1, 2, 3, etc. (the higher ω k , the weaker is the parametric resonance) . To demonstrate the technique which will be extensively explored in this paper, let us discuss the first resonance near ω k ≈ 1.

C. First resonance band
Let us follow the technique of [10] and find out the first resonance band near ω k ≈ 1, where the parametric resonance is strongest. The result in this subsection will not be used in other part of the papers, sine the resonance band lies in the region of unphysical values of ω k . However, the techniques used here to find the first resonance band can be applied, with minor modifications, to the second resonance and the Hartree-Fock approximation, so they are worth detail consideration.
First, let us substitute eq.(7) to eq.(10). For our purpose, only the leading term in eq. (7) is important. One obtains, The idea is to find f (k, t) in the form of a function with oscillatory behavior with the same frequency ω as φ 0 , but the amplitude of the oscillations is supposed to be a slowly varying function of time. So we look for the solution in the form where α(t) varies in a time scale much larger than 1, and dots stay for terms of higher harmonics. Substituting eq.(13) to eq.(12), and neglecting terms with second derivatives of a and b with time, as well as those with higher harmonics, we obtain Since α and α * vary much slower than e −iωt and e iωt , eq.(14) can be satisfied only if both the coefficients of e −iωt and e iωt vanish. Thus, one obtains the following equations, The solution to this set of equation can be searched in the form α(t) = α 0 e s k t , α * (t) = α * 0 e s k t . One finds Note that if one considers only the values of ω k so that ω k − ω ∼ A 2 0 , then s k ∼ A 2 0 ≪ 1, and our assumption that α varies on a time scale much larger than 1 is justified (in our case this time scale is A −2 0 ). Depending on whether |ω 2 k − ω 2 + A 2 0 | is larger or smaller than A 2 0 /2, s k is purely imaginary or real. If s k is purely imaginary, the solutions cannot be arbitrary large. The parametric resonance occurs only when s k is real, for example, the solution corresponding to positive s k can be arbitrary large if t is large enough. From eq.(16) one finds that the parametric resonance occurs when ω k lies in a finite window, However, it is easy to note that since ω 2 = 1 + A 2 0 2 , the upper limit of the range (17) coincides with 1. So, eq.(17) corresponds to unphysical values of ω k , ω k < 1, and there is no parametric resonance at physical ω k .
So far, our calculations are performed only to the first order on A 2 0 and do not exclude the possibility that the upper limit of the resonance band is differ from 1 by a value of O(A 4 0 ), so there is still a very small region of physical ω k where there is parametric resonance. However, the coincidence of the upper limit of the band and the boson mass can be shown to be hold exactly, to any order of A 0 . In fact, at ω k = 1, the mode equation (10) possesses a solution f (0, t) =φ 0 (t) that is a periodic function of t with the same period as φ 0 . This fact shows that at k = 0 one has s k = 0 exactly.

D. Second resonance band
Let us turn to the second resonance at ω k ≈ 2. Now one must take into account the second term in the expansion (7) and the linearized equation for f (k, t) is In analogy with the case of first resonance band, we will seek the solution to eq.(18) in the where we assume that α, β and γ are slowly varying functions that vary in a typical time scale of A −4 0 . We also expect that the coefficients of non-leading harmonics, β and γ, are suppressed by a factor of A 2 0 as compared to the leading harmonics, α. Substituting the ansatz (19) to eq.(18) and set the coefficients of e −2iωt , e −4iωt and 1 to zero, one obtains the following equations, respectively, From the last two equations one finds β = The first equation now can be written as 12 . This equation, except from the difference in numerical coefficients, has the same form as that of the first resonance (eq.(15). As is this case, there are two solutions proportional to e s k t where The resonance band is then Since ω c is approximately 2, the entire resonance band lies entirely in the physical region.
Some remarks are in order: 1. The center of the resonance band is located at ω k = ω c ≈ 2 + A 2 0 /4, i.e. is displaced by A 2 0 /4 from 2.
2. The width of the resonance band is of order A 4 0 , and is much smaller than the distance from the center of the resonance band to 2. The rate of exponential enhancement s k is also of order A 4 0 . There may be resonance bands at higher ω k , however we will ignore them since they are weaker than the one at ω k ≈ 2.
and since the energy carried by the modes with momentum k is proportional to |f (k, t)| 2 , the energy distribution has the Gaussian form centering at ω k = ω c and with width of order t −1 . We will call the regime when eq.(24) is valid the linear regime, since the modes obey the linearized equation on the background φ 0 . Now let us find out the region of validity of this regime. One may expect that the linear regime is valid once the energy of the inhomogeneous modes is much smaller than that of the homogeneous one. However in fact the linear regime breaks down when the energy carried by the inhomogeneous modes is still much smaller than the total energy. The reason for this is that the resonance band in our case has a small width proportional to A 4 0 , and a small change of the amplitude of homogeneous oscillations can drive the resonance band to a new location that does not overlap with the original resonance band. Let us estimate the upper limit on the energy carried by non-homogeneous modes for the linear approximation to be valid. Due to the energy conservation, when a portion of energy is transfered to non-homogeneous modes, the amplitude of the homogeneous oscillation decreases, and let us denote this as A 0 → A 0 − δA 0 . The center of the resonance band is located at k = 2 + A 2 0 4 , and we required that this values is shifted by a value much smaller than the width of the resonance band. This leads to the condition Since δA 0 /A 0 is proportional to the fraction of energy carried by inhomogeneous modes, the latter must be smaller than A 2 0 of the total, for the linear regime to be valid, Taking into account eq.(24), one finds that the linear regime breaks down at After t = T 0 , the back-reaction of the created modes onto the homogeneous oscillations cannot be neglected. Note that at t ∼ T 0 , the width of the Gaussian distribution, eq.(24), is A 4 0 log −1/2 (1/λA 2 0 ), i.e. much smaller than that of the resonance band. Before turning to Hartree-Fock approximation to go further than the linear regime, let us make an important remark on the phase coherence of the wave functions f (k, t) of particles created during the linear regime. This coherence can be seen in the simple case of ω = ω c .
In this case s k = A 4 0 96 , and from eq.(22) one sees that the exponentially increasing solution corresponds to α having the phase π/4, i.e. α = α 0 e s k t e iπ/4 where α 0 is a real constant.
From eq.(19) one finds, Thus, except from the slowly varying exponential part, the wave function f (k, t) has a welldefined phase shift, namely π/4, with φ 2 0 (t). It is worth noting that this phase shift is the same for all vectors k pointing to all possible directions that have ω k = ω c . If one does not require that ω k = ω c ,, the phase shift is different, but since most excited particles lies in a narrow central part of the resonance band, this phase shift is approximately constant for all modes excited during the linear regime.

A. Hartree-Fock equations
In the previous section we have seen that the linear regime is violated rather early when the energy transfered to the inhomogeneous modes is still small. After the linear regime ceases to be valid, one must take into account back-reaction of the inhomogeneous modes to the homogeneous oscillations, in particular the decrease of the amplitude of the latter due to the energy lost, and also the interaction between the modes. In this section, we will apply the Hartree-Fock approximation to this problem (the Hartree-Fock approximation was discussed earlier in connection with the problem of reheating in Ref. [7,2]). The region when this approximation is justified will be discussed as well.
The Hartree-Fock equations has the form, where φ 0 andφ are the classical and quantum parts of φ, respectively. Recalling eq.(9), the mean value ofφ 2 is where f is the mode function that is still to be found. It is convenient to decompose f into real and imaginary parts, where F 1 and F 2 are both real. Eqs.(26) can be rewritten in the form The initial conditions for these equations are formulated at t = 0. For φ 0 these conditions while for F 1,2 one has from eq.(11) In analogy with the technique used in previous section, we will look for the solution in the following form, where Ω(t) = dt ω(t) and the functions ω(t), A(t), B(t), α i (k, t), β i (k, t), and γ i (k, t) are assumed to be slowly varying functions of t. Note that ω(t) can chosen so that A(t) is always real, in this case ω(t) has the meaning of the time-dependent frequency of the φ 0 oscillations. Other functions, i.e. B, α i , β i are, in general, complex. The technique used here is similar to that applied to the case of first and second resonance bands. The calculation is rather cumbersome, but straightforward (see Appendix A). For our purpose, we need only the equations determining the amplitudes of the main harmonics of φ 0 and F i (k, t), i.e. A(t) and α i (k, t). The equation where i = 1, 2, and while A(t) can be determined from α i (k, t) by energy conservation, Note that when I and C are still much smaller than A 4 0 , eq.(32) reduces to eq.(20), and the linear regime is restored. It is useful to rescale ω k − ω c , t and α to expel the small parameter A 0 from eq.(32). Making use of the following change of variables eq.(32) becomes the following simple equation for a(τ, κ) where now The initial conditions for a i (κ, t) can be found from eqs.(30), (31) and (33), It is also useful to write eq.(34) in real notation. Denoting the real and imaginary parts of a i as u i and −v i , respectively, where u i and v i are real, the equations for u i and v i arė where The initial conditions for u 1,2 , v 1,2 are Let us briefly consider the limit when log(1/λ) is very large. More precisely, we consider the limit log(1/λ) → ∞, while t − T 0 = fixed (T 0 is the time moment when the linear regime breaks down, eq.(25)). In this limit, only modes close to the center of the resonance band are important, so in eq.(34) one can replace κ with 0. Then their solution has the form is arbitrary function of κ, which will be normalized for convenience, dκ f 2 (κ) = 1, and U i and V i do not depend on κ and satisfy the following equations, Giving the initial conditions, this set of equation can be solved numerically. The typical time dependence of the total energy of inhomogeneous modes, which is proportional to , is presented in Fig.1. As one can see from this figure, the energy of carrying by inhomogeneous modes first grows exponentially, but at some time moment reach a maximum and then drops back to 0. Note that at the maximum the total energy of inhomogeneous modes is of order A 2 0 of the total energy. If one follows eqs.(38) further, the energy will reach some very small minimum and then the process would repeat again, however this occurs outside the regime we are considering (t − T 0 = fixed).
The reason why the energy of inhomogeneous modes goes back to the homogeneous one after some time moment is that the behavior of the system in this limit is the same as of the system of two coupled oscillators. In fact, imagine that the particle spectrum consists of only modes with momenta k and −k, then if one repeats the whole calculation we have made one finds again eq.(38). The behavior shown in Fig.1 is typical for the system of two non-linearly coupled oscillators, where the energy is concentrated at first in one oscillator and the frequency of the second one lies in a resonance band on the background of the first.
C. Finite log(1/λ). Rate of particle production The regime discussed above (log(1/λ) → ∞, t − T 0 = fixed) is not, however, interesting from the physical point of view. In fact, even when the coupling λ is small, log(1/λ) is typically not very large. On other hand, we are interested in the behavior of the system up to thermalization, so t − T 0 may be very large. In these cases, the effect of the finite width of the distribution of resonance modes is important, and the system cannot be reduced to that of two coupled oscillators.
Since now the equation involves an infinite number of degrees of freedom, it seems hard to find the exact solution analytically. We have made numerical simulations of the solution.
In Fig.2 the time dependence of the total energy of inhomogeneous modes (I) is presented at λ = 10 −12 . As seen from this figure, at the beginning the energy of inhomogeneous modes is very small by it increases exponentially. At some time moment, however, it reaches some maximum, and the energy decreases. It does not reach 0 however, since log(1/λ) is finite.
The energy makes some more oscillations, and subsequent evolution is rather complicated, but shows a clear tendency to grow.
Fortunately, it is possible to find analytically the asymptotics of I in the limit of large τ , if one makes some reasonable assumptions on the behavior of I(τ ) and C(τ ). Let us suggest that I(τ ) and C(τ ) are smooth, slowly varying functions of τ at large τ . We assume that the imaginary part of C is very small and can be neglected, and ReC < 1. The last assumption is that I(τ ) + C(τ ) and I(τ ) − C(τ ) are monotone, increasing functions of τ . To find the asymptotics of I(τ ) and C(τ ), we will try to solve eq.(34) at each κ and demand latter that the conditions (35) are satisfied, so that our calculations are self-consistent. For the reason that will be clear in further discussions, we are interested only in large negative values of κ.
First assume that I and C are constant. Then eq.(34) has the same form as eq.(20), except that the center of the resonance band is now located at κ = I and its width is 1 − C. to be the time-dependent rate of exponential growth, and if s k (τ ) is slowly changing the enhancement factor during the whole period from τ = τ 1 to τ = τ 2 is proportional to exp dτ s κ (τ ) = exp π 2 · (1 − C) 2 I whereİ = ∂I/∂τ we have assumed thatİ and C can be considered as constant when τ varies from τ 1 to τ 2 . Since this enhancement factor must be of order λ −1/2 for a i to increase from √ λ to 1, one obtains the following equation, which is valid with the logarithmic precision.
3. τ > τ 2 , when κ is again outside the resonance band. Now u and v are again oscillatory and satisfy the equationu Making use of the fact that I and C are slowly varying functions, the solution to eqs. (40) has the form where and λ i (τ ) are slowly varying functions, and ϕ i (0) are some constants that will not be important in our further calculations. The time dependence of λ(τ ) can be found, making use of the adiabatic invariance technique [10]. The results reads where λ i (κ) does not depend on τ . So, we find the solution to eq(40) for given functions I(τ ) and C(τ ), Now we have to satisfy the self-consistency conditions, eqs.(37). At large τ one has I ≫ (1 − C) and both integrals are dominated by the region −I + (1 − C) < κ < 0. The first self-consistency condition reads, Now let us note that as cos(ϕ i ) and sin(ϕ i ) are rapidly oscillating functions, one can replace cos 2 (ϕ i ) and sin 2 (ϕ i ) by 1/2 and obtains Since this equation must be valid for any large value of τ , or equivalently, for any large negative value of −I + (1 − C), one finds λ 2 1 (κ) + λ 2 2 (κ) = 1 for every κ, |κ| ≫ 1. The second condition then reads, The integral in the r.h.s. can be evaluated to the logarithmic accuracy, and as the result one obtains the following equation, Note that both the integrals in (43) and (42) So, at large τ , I(τ ) grows almost linearly, except from a logarithmic correction. In further calculations we will be interested mostly in the order of magnitude of physical quantities, and in all estimations we will neglect the logarithms (i.e. we consider log λ −1 and log A −1 0 as quantities of order 1). Recall that I(τ ) is proportional to the density of particles created by parametric resonance, and restore all power of A 0 from eq.(33), one can sees that eq.(45) means that at t ≫ A −4 0 the rate of particle creation is of order λ −1 A 8 0 . Note that this rate is larger than that predicted by the elementary theory, eq.(4), by a factor of 1/λ (so one haṡ n 0 ∼ λ 2 n 4 0 instead of λ 4 n 4 0 . To summarize, we find that in the long run the resonance band moves toward lower energies (I increases with t and the center of the resonance band is located at κ = −I).
Due to this movement, particles with lower and lower energies are created, which in its turn causes the resonance band to move to smaller energies. The width of the resonance band decreases during this process (as log −1 τ ). It is worth noting, however, that the width of the particle spectrum, though much larger than that of the initial resonance band, i.e. A 4 0 , is still much smaller than 1, so all created particles have energy close to 2. The excitation of a given mode by parametric resonance is a coherent process, but one the parametric enhancement terminates, the phase of the mode become rapidly uncorrelated with that of the homogeneous oscillations and of the other modes. At t ≫ A −4 0 , most of the modes that have been excited lie outside the resonance band and can be considered as ordinary particles.
It is useful to find the distribution function of the particles produced by this mechanism. This can be done by comparing eq.(27) with the one giving the fluctuation φ 2 for the state characterized by the distribution n k , If n k ≫ 1 (as will be obtained below), one finds where the average in the r.h.s. can be taken over time or over a finite volume in the kspace. Now consider some mode that the resonance band has gone by. At this mode one has The occupation number is then Therefore, at given time moment t ≫ A −4 0 , the particle distribution has the form of a thin spherical shell. The outer sphere of this shell is fixed at ω k = ω c , while the inner sphere moves inward when t increases. The occupation number of modes inside the shell is of order 1/λ. Now let us estimate when the number of particle with energy E = 2 when their scattering cannot be neglected. Consider the scattering process where two particles with energy E = 2 scatter to one particle in the condensate (i.e. with energy 1) and another particle with energy E = 3. This process is in fact the most probable since it involves one particle in the condensate, and thus has a large Bose-enhancement factor associated with the latter. The criterion for the scattering to be important is that the rate of production of particles with energy E = 3 is comparable with those with energy E = 2. Let us suppose that the density of the particle with energy E = 2 is n. Then the number of particle with energy E = 3 created by the scattering process in an unit volume during an unit time interval is σn 2 n 0 , where σ is the cross section and n 0 is the density of the condensate. Substituting σ ∼ λ 2 and n 0 ∼ A 2 0 /λ, one finds that this rate is λA 2 0 n 2 . Compare with the rate of creation of particles with energy E = 2, which is of order λ −1 A 8 0 , one finds that the scattering becomes important when n ∼ λ −1 A 3 0 , which occurs at the time moment t ∼ A −5 0 . Note that at t ∼ A −5 0 the total energy of inhomogeneous modes is λ −1 A 3 0 , which is still much smaller than the total energy (λ −1 A 2 0 ).

IV. PARTICLE SCATTERING AND FORMATION OF KOLMOGOROV TURBULENCE SPECTRUM
In the previous section we have found that the Hartree-Fock approximation breaks down at t ∼ A −5 0 , when most of the energy is still carried by the homogeneous mode. In this section we will consider the behavior of the system after the moment t ∼ A −5 0 up to the time when the the homogeneous mode begins to lost a considerable portion of its energy. Since the energy flux from the homogeneous mode is λ −1 A 8 0 and the total energy is of order λ −1 A 2 0 , the latter time moment is A −6 0 . Thus, in this section we consider A −5 0 ≪ t ≪ A −6 0 . The system can be characterized by the distribution function n k (t), while the homogeneous oscillations correspond to the delta-functional part of n k at k = 0.
Let us first write down the generic Boltzmann equation (the Boltzmann equation for the φ 4 field theory has been studied in a different setting, see Refs. [12,13]). Since most of our discussion will be based on order-of-magnitude estimations, we will neglect all numerical factors. One writes, In fact, one should take into account the particle creation by parametric resonance, however let us first neglects this effect. Since the occupation numbers are typically of order 1/λ which are much larger than 1, one can neglect the terms 1 compared to n. The Boltzmann equation now becomes dn k dt ∝ λ 2 dk 1 dk 2 dk 3 ω k ω k 1 ω k 2 ω k 3 (n k 1 n k 2 n k 3 + n k 1 n k 2 n k − n k n k 3 n k 1 − n k n k 3 n k 2 ) Let us make two remarks on eq.(46). First, since we have n k ∼ λ −1 , the time scale set by the Boltzmann equation does not depend parametrically on λ (but may depend on A 0 ).
Second, if at some t the mode k is unoccupied, n k = 0, the particles will be created at this k only when there is a process k 1 , k 2 → k 3 , k where the other three momenta participating in the process have nonzero occupation numbers. Since the initial condition for the Boltzmann equation is that where all particles have energy E = 1 or 2, approximately, one can expect that the particle distribution will have peaks at E = 3, 4... i.e. at integer values of E when t ∼ A −5 0 . In fact, if three particles participating in a scattering have integer energies, the energy of the fourth is also integer due to energy conservation. This "peak-like" structure of the spectrum will not preserve infinitely long: since the particles produced by parametric resonance have not exactly the same energy, the width of the peaks will become broader and broader, and at some time moment comparable to 1, however this happens at some time scale much larger than A −5 0 . Therefore, we take the distribution function in the form where we have introduced the notation k l = √ l 2 − 1, which is the momentum corresponding to the energy ω k = l. Note that n l has the meaning of the density of particles with energy l. Before substituting eq.(47) to eq.(46), let us also notice that the scattering processes can be divided into 2 types: those without participation of condensate particles (i.e. all k 1 , k 2 , k 3 are nonzero), and those with one condensate particle participating in the scattering (one of k 1 , k 2 , k 3 is 0). In order to avoid the overlap with the Hartree-Fock effects, it is necessary to exclude from consideration the scattering where two of the participating particles have zero momenta. It is natural to expect that when the density of non-condensate particles is still small, the processes of the second type (which will be called "three-particle" processes) dominate over those of the first type ("four-particle" transitions). The criterion for the 4-particle scattering to be negligible compared to the 3-particles ones can be easily found. In fact, the 3-particle processes are enhanced by the Bose factor of n 0 due to the participation of a condensate particle in the scattering, while the 4-particle processes contain instead of this factor an additional momentum integral dk n k /ω k . So, the criterion for the 4-particle collision to be negligible is Running ahead, let us state that this condition is satisfied in the whole evolution up to the moment when thermalization is completed.
At first, one can consider the density of particles in the condensate n 0 as a constant.
Substituting eq.(47) to the Boltzmann equation, and throwing away terms corresponding to four-particle collisions, after some calculations (see Appendix B) the following equation is obtained, n l 1 n l 2 k l 1 k l 2 l 1 l 2 + l 1 −l 2 =l−1 n l 1 n l 2 k l 1 k l 2 l 1 l 2 − n 2 l k 2 l l 2 − 2n l k l l l 1 <l n l k l l   + jδ l2 (49) where we have included into the r.h.s. the term jδ l2 to take into account the production of particles with energy closed to E = 2 that we have found in the Hartree-Fock approximation.
Therefore, the Boltzmann equation is now written in terms of a discrete set of numbers n l which is proportional to the density of particles with energy closed to l. Let us emphasize that in eq.(49) we do not write down the numerical constants. For estimation, it is sufficient to take a constant j ∼ λ −1 A 8 0 , although in fact we found in Sect.3 that the production rate is slowly changing.
As seen from eq.(49), due to the scattering, particles with energy E = 3 are created by colliding two particles with energy E = 2. Then particles with energy E = 4 and 5 come to existence in the collision of particles with energy E = 2 and 3, etc. Therefore, particles with higher and higher energies are created, despite the fact that the source term in eq.(49) produces particles with E = 2 only.
The situation when energy is transfered from the long-distance modes to short-distance ones is the common situation in turbulent systems (see, for instance, [11]). Hence one can make use of the vast intuition elaborated in the theory of turbulence in order to predict the long-term behavior of the system described by eq.(49). In particular, in the turbulent systems, if one has a constant flow of energy into the modes at some momentum scale during an infinite time interval, a stationary distribution is formed. To find this distribution in our case, one looks for those which has vanishing collision integral, While eq.(50) seems hard to solve analytically, the asymptotics of n l at large l can be found.
In analogy with the Kolmogorov spectrum of turbulence, the stationary solution can be expected to have the form n l ∼ l η at large l. For convenience, it is useful to introduce n(E) = n l 4πl 2 , l ≈ E which (at large E) has the meaning of the average occupation number n k at energy ω k = E (one can imagine n(E) as the distribution function after one smears off the peaks of the spectrum). At large l, the sums in eq.(49) can be transformed to integrals, and one obtains, Substituting the ansatz n(E) ∼ E η to this equation, one obtains the following equation on the Kolmogorov index η for n(E) to be stationary, The l.h.s. is finite when η lies in the region −3 < η < −1/2. It can be verified that eq.(51) has exactly two solutions, η = −1 and η = −3/2 (accidentally, these values coincide with the Kolmogorov indices in the case of non-relativistic Bose gas with condensate, see Ref. [13]). The physical meaning of the first solution is trivial: it corresponds to the thermal distribution n k ∼ ω −1 k , and is not appropriate for our purposes. The second solution is a non-trivial one, n k ∼ ω −3/2 k , and describes the particle distribution developed in our system after an infinite time interval.
In fact, this stationary spectrum can be established only if one have the energy flux to the mode with E = 2 during an infinitely long time period (this can be seen from the fact that the stationary spectrum n k ∼ ω −3/2 k has infinite energy). In the real situation when this time interval is finite, one expects that the Kolmogorov spectrum is cut off in the ultraviolet.
In other words, at finite t one has n k ∼ ω −3/2 k at k ≪ k 0 (t) and n k ≈ 0 at k ≫ k 0 (t). The cutoff k 0 is expected to increase with t. Note that for this type of spectrum both energy and particle number is concentrated at the cutoff, k ∼ k 0 .
That the Kolmogorov index is equal −3/2 can be understood in a very intuitive way.
Suppose that at some time moment the particle distribution has the form of n k ∼ ω η k , with a cutoff at k 0 . If η is the correct Kolmogorov index, the cutoff k 0 increases with time, while the occupation numbers at k ≪ k 0 do not change. The time scale during which the cutoff is, say, doubled, is of the same order as the time scale set by the Boltzmann equation, which can be estimated to be proportional to k 0 /n(k 0 ) ∼ k 1−η 0 . During this time interval, an amount of energy, also proportional to k 1−η 0 is fed into the system, which is spent to the creation of particles with energies between k 0 and 2k 0 . The total energy of these particles is n(k 0 )k 4 0 ∼ k 4+η 0 . So, we obtain the equation 1 − η = 4 + η, which gives the Kolmogorov The power law n k ∼ ω −3/2 k is valid only at k ≪ k 0 . Let us consider the particle spectrum at the momentum scale of the cutoff, k ∼ k 0 . In analogy with the situation in turbulence systems [11][12][13], one can expect that the distribution has the following self-similar behavior, where c 1 and c 2 are some constants chosen so that the function f (ξ) varies on the typical values of ξ of order 1 and is itself of order 1. We will try to find p and q but will not attempt to obtain the form of the function f (ξ). Substituting this ansatz to the Boltzmann equation, and comparing the time dependence of the left and right hand sides, one finds that the latter can be satisfied only if the following relation between p and q holds, On the other hand, one expects that when k ≪ k 0 eq.(52) must reproduce the Kolmogorov distribution. This fact implies that the asymptotics of the function f (ξ) at small ξ is f (ξ) ∼ ξ −3/2 . Furthermore, in the same region k ≪ k 0 , the spectrum must be stationary, i.e. n k does not depend on t. Taking into account these two facts, one obtains one more equation on p and q, Solving eqs. (53) and (54), we find p = 2/5, q = 3/5. In particular the cutoff k 0 grows with time as k 0 ∼ t 2/5 . The total energy of the system behaves like c 1 c −4 2 t 4p−q ∼ c 1 c −4 2 t, which agrees with the fact that the energy is fed into the system with a constant rate. Since the rate of creation of particles with energy close to 2 is On the other hand, at t ∼ A −5 0 , the typical momentum of particles is 1, so c 2 = A −2 0 . Therefore, c 1 = λ −1 , and eq.(52) looks like the following, Let us emphasize that in our calculations are limited to the range A −5 0 ≪ t ≪ A −6 0 . It is straightforward to verify the condition for the 4-particle collisions to be negligible: the l.h.s. of eq.(48) is of order λ −1 A 2 0 while the r.h.s. is always smaller than A 14/5 0 in the time interval we are considering, A −5 0 ≪ t ≪ A −6 0 . Another note is in order. The total kinetic energy of the produced particles grow almost linearly with time, but since the average energy of the particles also grows as t 2/5 , the number of particles grows only as t 3/5 . This means that the larger is t, the slower is the growth of the total number of particles outside the condensate. This slowdown occurs despite the fact that particles with energy E = 2 are produced with approximately constant rate. There is in fact no contradiction since the dominant collision process in our case is the 2 → 2 scattering where two particles with nonzero momenta collide and become one particles in the condensate and one particle with nonzero momentum: each such act effectively reduces the number of particles with nonzero momentum by 1. At large t most particles produced at the energy level 2 effectively go back to the condensate by this mechanism, transferring their kinetic energy to the particles that are already outside the condensate. The condensate looses 4 particles in each 4 → 2 transition but gains back 2 particles from the scattering, therefore the actual rate of condensate evaporation is smaller than predicted by the Hartree-Fock approximation by a factor of 2.
For further reference let us note that at the end of the epoch that we are considering, . Since now the typical energy is much larger than 1, one can neglect the peak structure of the spectrum and deal with the average occupation number n(E) only, regardless how effective the processes broadening the peaks are.

V. THERMALIZATION
A. Early step of thermalization: When t become much larger than A −6 0 , the energy carried by the condensate becomes much smaller than the total. Therefore, one can consider the regime t ≫ A −6 0 as the epoch of thermalization of particles produced by condensate evaporation. Though the condensate now carries only a small fraction of the energy, it still plays an important role in the scattering process. In fact, we will find that the 3-particle collisions are still dominant. Let us first find out the time dependence of the condensate density. To do that let us first assume that the the dominant mechanism of the decay of the condensate density is the coherent particle creation via parametric resonance. Then the decay of n 0 is describes by the equatioṅ n 0 = λ 3 n 4 0 , from which one finds The Boltzmann equation, where only 3-particle collisions are taken into account, now has the form It is natural to assume that the particle distribution still shows a self-similar behavior (52) in this regime with a some values of the parameters p, q, c 1 and c 2 . Substituting the ansatz (52) into eq.(56) one obtains the following relation between p and q, On the other hand, the energy should be conserved (since most energy has been emitted from There are two assumptions that have been used in order to find eq.(57) that we now have to check. The first is that particles in the condensates are lost mainly due to the parametric resonance condensate evaporation. To find when this assumption is justified, let us mention that the concurrent process that changes the density of the condensate is the particle collision: as noted above, the dominant scattering is supposed to be those where 2 non-condensate particles collide and becomes one non-condensate particle and a condensate one and the number of particle in the condensate is increased by 1. Therefore, there is a flux of particles coming from the non-condensate modes back to the condensate, which can be calculated by taking the time derivative of the total number of non-condensate particles.
According to eq.(57), the total number of particles with nonzero momentum is 0 /t 2/15 , so the flux of particles to the condensate is λ −1 A 8/5 0 /t 17/15 . On the other hand, the rate of particle creation by parametric resonance is n 4 0 ∼ λ −1 t −4/3 . Comparing the two quantities, one sees that the former process can be neglected when t ≪ A −8 0 . The Bose condensation becomes comparable to condensate evaporation at t ∼ A −8 0 . So, the region of validity of eq.(57) is A −6 0 ≪ t ≪ A −8 0 . Another requirement that must be verified is eq.(48). The l.h.s. is λ −1 t −1/3 while the r.h.s. is of order λ −1 A 6/5 0 /t 4/15 . It is easy to see that the l.h.s. is larger than the r.h.s. everywhere in the region A −6 0 ≪ t ≪ A −8 0 . According to eq.(57), at the end of the regime, t ∼ A −8 0 , the typical energy of particles is A −2/3 0 , while the typical occupation number is λ −1 A 14/3 0 . The density of the condensate is At t ∼ A −8 0 , the flux of particles coming to the condensate from the hard modes becomes comparable to the lost of particles due to the parametric resonance. One can expect, therefore, that the decay of n 0 is now slower than t −1/3 . We make the following assumptions in the regime t ≫ A −8 0 : 1. n 0 ∼ t −r , where r is some unknown constant.
2. The particle distribution outside the condensate has the self-similar form (52) with some p and q.

The loss of particles by condensate evaporation is compensated almost exactly by
Bose condensation. The change of n 0 with time is due to a small discrepancy between the rates of the two processes.
4. The four-particle collisions can be neglected compared to three-particle ones.
Let us find p, q and r. First, since the four-particle scattering can be neglected, the Boltzmann equation has the same form as eq.(56), where the factor of t −1/3 in the r.h.s. is replaced by t −r , so p + q + r = 1 (58) As before, the energy conservation implies Finally, the rate of particle creation by parametric resonance is n 4 0 ∼ t −4r . On other hand, the number of particles outside the condensate is t −p , so the flux of non-condensate particles coming back to the condensate is t −(1+p) . According to our third assumption, the following relation should hold.
As anticipated, n 0 decays slower than t −1/3 . Now one can verify our starting assumptions.
First,ṅ 0 ∼ λ −1 A 8/21 0 /t 9/7 is much smaller than the rate of particle creation by parametric Therefore, particle loss is much smaller than the rate of evaporation by parametric resonance, which implies that the latter is compensated almost exactly by Bose condensation. That the four-particle collisions can be neglected compared with three-particle ones can be checked easily: the l.h.s. of eq.(48) is A 8/21 0 /t 2/7 , while the r.h.s. is A 22/21 0 /t 2/7 . Since A 0 ≪ 1, this condition is always satisfied. Therefore, all our starting assumptions are justified a posteriori. It is interesting to note that the fourparticle scattering is negligible in the whole history of the evolution of the system: even when the density of the condensate is small, the dominant relaxation processes are those involving one condensate particle. The remnant of homogeneous oscillations, so, play an important role in the thermalization.
The thermal equilibrium is achieved when the typical occupation number is of order 1.
This occurs at the the time moment Naturally, the smaller the coupling, the larger the time required for thermalization. The residual particle density in the condensate at these times is For comparison, the density of particles outside the condensate is A 3/2 0 λ −3/4 , which is much larger than n 0 . Naturally, the condensate density will decrease further and tend to 0.

VI. CONCLUSION
In this paper we have followed the evolution of the scalar field from the initial homogeneous oscillations to thermalization. We found that this evolutions is a rather complex process, occurring through following regimes: 0 : exponential enhancement of resonance modes. This is the feature that makes our model similar to the ones considered in Ref. [2].

A −4
0 ≪ t ≪ A −5 0 : particles with energy ω k ≈ 2 are produced with a roughly constant rate.

A −5
0 ≪ t ≪ A −6 0 : the typical energy of non-condensate particles is k 0 ∼ A 2 0 t 2/5 , and the typical occupation number is n k ∼ λ −1 t −3/5 . Most energy is still concentrated in the condensate, whose density remains approximately constant, : the condensate density decays as n 0 ∼ λ −1 t −1/3 . The typical energy and occupation number of non-condensate particles are k 0 ∼ A 2/5 0 t 2/15 and n k ∼ The thermalization finishes at t ∼ A 1/6 λ −7/4 . We see that in our model the Hartree-Fock approximation works only in the regime 1 our calculations give the rate of particle loss from the condensate twice smaller than predicted by the Hartree-Fock approximation. Finally, at t ≫ A −8 0 the scattering changes the decay law of the condensate, so instead of n 0 ∼ t −1/3 one has n 0 ∼ t −2/7 . Though all calculations in this paper are based on the existence of the additional small parameter A 0 , one may expect that some features similar to those discovered in our model also occur in the some of the realistic theories. As one of them, we speculate that the exponential enhancement of the parametric resonance modes is terminated by the backreaction of produced modes and the coherent nature of the latter may slow down the energy loss of the homogeneous oscillations. It is also natural to argue that the thermalization process takes a long time, during which some types of self-similar behavior of the distribution functions may have place in some models. Hopefully, the techniques used in this paper may prove helpful in the understanding of the dynamics of the inflaton in these theories, which would provide more accurate physical predictions than those provided by the Hartree-Fock approximation.
It is particularly interesting to compare our analytical results with the numerical simulations in the massless φ 4 theory performed in Ref. [9]. It is worth noting that though the massless φ 4 theory lacks any additional small parameter beside λ, the width of the resonance band, as well as the value of the maximal index of exponential enhancement s k are numerically small in this model. Therefore, one can expect that many qualitative features found in our model are also valid for the massless φ 4 theory. Though late thermalization epochs are not considered in Ref. [9], the behavior of the field show many similarity with our model in earlier regimes. In particular, the fluctuations of the field comming from the non-condensate modes are relatively small at the end of the linear regime. The peaks in the spectrum are also found to be present during some time interval and they are shown to be caused by scattering processes. From the analogy with our model, the bounces found in the time evolution of the amplitude of homogeneous oscillations can be attributed to the coherence of resonance modes, not to Bose condensation as argued in Ref. [9].
The author thanks P. Arnold and L. Yaffe for stimulating discussions, and L. Kofman and A. Linde for explaining many points of Ref. [2]. I thank S.Yu. Khlebnikov and I.I. Tkachev for numerous comments on the manuscript. This works is supported, in part, by the U.S. Department of Energy under grant #DE-FG03-96ER40956. In this appendix we will show how eq.(32) can be derived. For convenience, we will to work with complex notations. In analogy with the technique used in Sect.2 to find the first and second resonance band, we will try to find the solution to the Hartree-Fock equation in the following form, φ 0 = 1 √ λ Ae −iΩ + Be −3iΩ + h.c. + · · · F i (k) = 1 √ λ α i (k)e −2iΩ + β i (k)e −4iΩ + δ i (k)e −6iΩ + h.c. + γ i (k) + · · · (A1) where Ω = dt ω(t) and ω(t) is the time-dependent frequency of the oscillations that is chosen so that A is real.
The order of magnitude of the parameters in eq.(A1) are expected to be as follows (the consistence of these relations can be verified a posteriori), Let us substitute the ansatz (A1) to the Hartree-Fock equations, eqs. (28) and (29). First let us consider eq.(29). To the order of A 6 , the l.h.s. is expanded into terms proportional to 1, e ±2iΩ , e ±4iΩ and e ±6iΩ . Each coefficient of these exponents must be equal 0, which gives us 7 equations. Those corresponding to e ±6iΩ are the equations that determine d k . Since d k do not emerge in other equations, these equations can be ignored. It turns out that we will need the coefficient of 1 and e ±4iΩ only with the accuracy of O(A 4 ) (not O(A 6 )). To this order, terms proportional to 1 read 4γ i (k) + A 2 2 (α i (k) + α * i (k)) = 0 while the equations corresponding to e ±4iΩ are −12β i (k) + A 2 2 α i (k) = 0 and its complex conjugation. These equations relate β i (k) and γ i (k) with α i (k), The term e −2iΩ implies the following equation − 4iα i (k) + −4ω 2 + ω 2 k + A 2 + I α i (k) + (AB + C) α * i (k) + A 2 2 (β i (k) + γ i (k)) = 0 (A3) where we introduce the notations I = dq (2π) 3 2ω q (|α 1 (q)| 2 + |α 2 (q)| 2 ) C = 1 2 dq (2π) 3 2ω q (α 2 1 (q) + α 2 2 (q)) Eq.(A3) can be rewritten in the following form after taking into account eqs.(A2) Now let us turn to eq.(28). To the order of O(A 7 ), both sides can be expanded in terms proportional to e ±iΩ , e ±3iΩ , e ±5iΩ and e ±7iΩ . We will not need the equation corresponding to e ±5iΩ and e ±7iΩ . Moreover we will need the equation coming from the term e −3iΩ with the accuracy of O(A 5 ) only, which has the form, (−9ω 2 + 1 + A 2 )B + A 3 6 + AC = 0 form which one obtains, The equation coming from the terms containing e −iΩ is, Together with eqs.(A2), (A5), and the condition that A is real, one obtain the formula for