# Analytic and Numerical Study of Preheating Dynamics

###### Abstract

We analyze the phenomenon of preheating, i.e. explosive particle production due to parametric amplification of quantum fluctuations in the unbroken symmetry case, or spinodal instabilities in the broken symmetry phase, using the Minkowski space vector model in the large limit to study the non-perturbative issues involved. We give analytic results for weak couplings and times short compared to the time at which the fluctuations become of the same order as the tree level terms, as well as numerical results including the full backreaction. In the case where the symmetry is unbroken, the analytical results agree spectacularly well with the numerical ones in their common domain of validity. In the broken symmetry case, interesting situations, corresponding to slow roll initial conditions from the unstable minimum at the origin, give rise to a new and unexpected phenomenon: the dynamical relaxation of the vacuum energy. That is, particles are abundantly produced at the expense of the quantum vacuum energy while the zero mode comes back to almost its initial value. In both cases we obtain analytically and numerically the equation of state which in both cases can be written in terms of an effective polytropic index that interpolates between vacuum and radiation-like domination.

We find that simplified analysis based on harmonic behavior of the zero mode, giving rise to a Mathieu equation for the non-zero modes miss important physics. Furthermore, such analysis that do not include the full backreaction and do not conserve energy, results in unbound particle production. Our results rule out the possibility of symmetry restoration by non-equilibrium fluctuations in the cases relevant for new inflationnary scenarios. Finally, estimates of the reheating temperature are given, as well as a discussion of the inconsistency of a kinetic approach to thermalization when a non-perturbatively large number of particles are created.

###### pacs:

11.10.-z;11.15.Pg;98.80.Cq^{†}

^{†}preprint: PITT-96-; CMU-HEP-96-; LPTHE-96/32

## I Introduction

It has recently been realized[2, 3, 4] that as the zero momentum mode of a quantum field evolves, it can drive a large amplification of quantum fluctuations. This, in turn, gives rise to copious particle production for bosonic fields, creating quanta in a highly non-equilibrium distribution, radically changing the standard picture of reheating the post-inflationary universe[5, 6, 7]. This process has other possible applications, such as in understanding the hadronization stage of the quark-gluon plasma[8] as well as trying to understand out of equilibrium particle production in strong electromagnetic fields and in heavy ion collisions [9, 10, 11, 12].

The actual processes giving rise to preheating can be different depending on the potential for the scalar field involved as well as the initial conditions. For example, in new inflationary scenarios, where the inflaton field’s zero mode evolves down the flat portion of a potential admitting spontaneous symmetry breaking, particle production occurs due to the existence of unstable field modes which get amplified until the zero mode leaves the instability region. These are the instabilities that give rise to spinodal decomposition and phase separation. In contrast, if we start with chaotic initial conditions, so that the field has large initial amplitude, particles are created from the parametric amplification of the quantum fluctuations due to the oscillations of the zero mode.

In this paper we analyze the details of this so-called preheating process both analytically as well as numerically. Preheating is a non-perturbative process, with typically particles being produced, where is the self coupling of the field. Due to this fact, any attempts at analyzing the detailed dynamics of preheating must also be non-perturbative in nature. This leads us to consider the vector model in the large limit. This is a non-perturbative approximation that has many important features that justify its use: unlike the Hartree or mean-field approximation[4], it can be systematically improved in the expansion. It conserves energy, satisfies the Ward identities of the underlying symmetry, and again unlike the Hartree approximation it predicts the correct order of the transition in equilibrium.

This approximation has also been used in other non-equilibrium contexts[9, 10, 11, 12]. In this work, we consider this model in Minkowski space, saving the discussion of the effects of the expansion of the universe for later work.

Our findings are summarized as follows.

We are able to provide consistent non-perturbative analytic estimates of the non-equilibrium processes occurring during the preheating stage taking into account the exact evolution of the inflaton zero mode for large amplitudes when the quantum back-reaction due to the produced particles is negligible i.e. at early and intermediate times. We also compute the momentum space distribution of the particle number as well as the effective equation of state during this stage. Explicit expressions for the growth of quantum fluctuations, the preheating time scale, as well as the effective (time dependent) polytropic index defining the equation of state are given in sec. III.

We then go beyond the early/intermediate time regime and evolve the equations of motion numerically, taking into account back-reaction effects. (That is, the non-linear quantum field interaction). These results confirm the analytic results in their domain of validity and show how, when back-reaction effects are large enough to compete with tree level effects, dissipational effects arise in the zero mode. Energy conservation is guaranteed in the full backreaction problem, leading to the eventual shut-off of particle production. This is an important ingredient in the dynamics that determines the relevant time scales.

We also find a novel dynamical relaxation of the vacuum energy in this regime when the theory is in the broken phase. Namely, particles are produced at the expense of the quantum vacuum energy while the zero mode contributes very little. We find a radiation type equation of state for late times () despite the lack of thermal equilibrium.

Finally we discuss the calculation of the reheating temperature in a class of models, paying particular attention to when the kinetic approach to thermalization and equilibration is applicable.

There have been a number of papers (see refs.[2, 3, 13] -[14]) dedicated to the analysis of the preheating process where particle production and back-reaction are estimated in different approximations [30].

The layout of the paper is as follows. Section II presents the model, the evolution equations, the renormalization of the equations of motion and introduces the relevant definitions of particle number, energy and pressure and the details of their renormalization. The unbroken and broken symmetry cases are presented in detail and the differences in their treatment are clearly explained.

In sections III through V we present a detailed analytic and numerical treatment of both the unbroken and broken symmetry phases emphasizing the description of particle production, energy, pressure and the equation of state. In the broken symmetry case, when the inflaton zero mode begins very close to the top of the potential, we find that there is a novel phenomenon of relaxation of the vacuum energy that explicitly shows where the energy used to produced the particles comes from. We also discuss why the phenomenon of symmetry restoration at preheating, discussed by various authors[3, 14, 26, 27] is not seen to occur in the cases treated by us here and relevant for new inflationnary scenarios [30].

In section VI we provide estimates, under suitably specified assumptions, of the reheating temperature in the model as well as other models in which the inflaton couples to lighter scalars. In this section we argue that thermalization cannot be studied with a kinetic approach because of the non-perturbatively large occupation number of long-wavelength modes.

Finally, we summarize our results and discuss future avenues of study in the conclusions. We also include an appendix where we gather many important technical details on the evaluation of the Floquet mode functions and Floquet indices used in the main text.

## Ii Scalar Field Dynamics in the Large limit

As mentioned above, preheating is a non-perturbative phenomenon so that a non-perturbative treatment of the field theory is necessary. This leads us to consider the vector model in the large limit.

In this section we introduce this model, obtain the non-equilibrium evolution equations, the energy momentum tensor and analyze the issue of renormalization. We will then be poised to study each particular case in detail in the later sections.

The Lagrangian density is the following:

(1) |

for fixed in the large limit. Here is an vector, and represents the “pions”. In what follows, we will consider two different cases of the potential , with () or without () symmetry breaking.

We can decompose the field into its zero mode and fluctuations about the zero mode:

(2) |

The generating functional of real time non-equilibrium Green’s functions can be written in terms of a path integral along a complex contour in time, corresponding to forward and backward time evolution and at finite temperature a branch down the imaginary time axis. This requires doubling the number of fields which now carry a label corresponding to forward (), and backward () time evolution. The reader is referred to the literature for more details[15, 16]. This generating functional along the complex contour requires the Lagrangian density along the contour, which for zero temperature is given by[4]

(3) |

The tadpole condition will lead to the equations of motion as discussed in [4] and references therein.

A consistent and elegant version of the large limit for non-equilibrium problems can be obtained by introducing an auxiliary field and is presented very thoroughly in reference[10]. This formulation has the advantage that it can incorporate the corrections in a systematic fashion. Alternatively, the large limit can be implemented via a Hartree-like factorization[4] in which i) there are no cross correlations between the pions and sigma field and ii) the two point correlation functions of the pion field are diagonal in the space of the remaining unbroken symmetry group. To leading order in large both methods are completely equivalent and for simplicity of presentation we chose the factorization method.

The factorization of the non-linear terms in the Lagrangian is (again for both components):

(4) | |||||

(5) | |||||

(6) | |||||

(7) | |||||

(8) |

To obtain a large N limit, we define

(9) |

where the large N limit is implemented by the requirement that

(10) |

The leading contribution is obtained by neglecting the terms in the formal large limit. The resulting Lagrangian density is quadratic, with a linear term in :

where,

(12) | |||||

(13) | |||||

(14) |

Note that we have used spatial translational invariance to write

(15) |

The necessary (zero temperature ) non-equilibrium Green’s functions are constructed from the following ingredients

(16) | |||

(17) |

while the Heisenberg field operator can be written as

(18) |

where are the canonical destruction and annihilation operators and the quantization volume.

The evolution equations for the expectation value and the mode functions can be obtained by using the tadpole method[4] and are given by:

(19) | |||

(20) | |||

(21) | |||

(22) |

The initial state is chosen to be the vacuum for these modes, i.e. . The frequencies (i.e. ) will determine the initial state and will be discussed for each particular case below.

The fluctuations obey an independent equation, that does not enter in the dynamics of the evolution of the expectation value or the fields to this order and decouples in the leading order in the large limit[4].

It is clear from the above equations that the Ward identities of Goldstone’s theorem are fulfilled. Because , whenever vanishes for then and the “pions” are the Goldstone bosons. This observation will be important in the discussions of symmetry breaking in a later section.

Since in this approximation, the dynamics for the and fields decouple, and the dynamics of does not influence that of , the mode functions or , we will only concentrate on the solution for the fields. We note however, that if the dynamics is such that the asymptotic value of the masses for and the “pion” multiplet are different, and the original symmetry is broken down to the subgroup.

### ii.1 Renormalization of the Model

We briefly review the most relevant features of the renormalization program in the large limit that will be used frequently in our analysis. For more details the reader is referred to[10, 4, 12].

In this approximation, the Lagrangian is quadratic, and there are no counterterms. This implies that the equations for the mode functions must be finite. This requires that

(23) |

Defining

(24) |

The function is written as linear combinations of WKB solutions of the form

(25) |

with obeying a Riccati equation[12] and the coefficients are fixed by the initial conditions. After some algebra we find

(27) | |||||

(28) |

Using this asymptotic forms, we obtain[4, 12] the following renormalized quantities

(29) | |||

(30) |

and

(31) |

where we have introduced the (arbitrary) renormalization scale . Eqs. 23 and 31 lead to the renormalization conditions valid in the large limit.

At this point it is convenient to absorb a further finite renormalization in the definition of the mass and introduce the following quantities:

(32) | |||

(33) | |||

(34) | |||

(35) | |||

(36) | |||

(37) |

For simplicity in our numerical calculations later, we will chose the renormalization scale . The evolution equations are now written in terms of these dimensionless variables, in which dots now stand for derivatives with respect to .

### ii.2 Unbroken Symmetry

In this case , and in terms of the dimensionless variables introduced above we find the following equations of motion:

(38) | |||

(39) | |||

(40) | |||

(41) |

As mentioned above, the choice of determines the initial state. We will choose these such that at the quantum fluctuations are in the ground state of the oscillators at the initial time. Recalling that by definition , we choose the dimensionless frequencies to be

(42) |

### ii.3 Broken Symmetry

In the case of broken symmetry and the field equations in the limit become:

(45) | |||

(46) |

where is given in terms of the mode functions by the same expression of the previous case, (44). Now the choice of boundary conditions is more subtle. The situation of interest is when , corresponding to the situation where the expectation value rolls down the potential hill from the origin. The modes with are unstable and thus do not represent simple harmonic oscillator quantum states. Therefore one must chose a different set of boundary conditions for these modes. Our choice will be that corresponding to the ground state of an upright harmonic oscillator. This particular initial condition corresponds to a quench type of situation in which the initial state is evolved in time in an inverted parabolic potential (for early times ). Thus we shall use the following initial conditions for the mode functions:

(47) | |||

(48) | |||

(49) |

along with the initial conditions for the zero mode given by eq.(41).

### ii.4 Particle Number

Although the notion of particle number is ambiguous in a time dependent non-equilibrium situation, a suitable definition can be given with respect to some particular pointer state. We consider two particular definitions that are physically motivated and relevant as we will see later. The first is when we define particles with respect to the initial Fock vacuum state, while the second corresponds to defining particles with respect to the adiabatic vacuum state.

In the former case we write the spatial Fourier transform of the fluctuating field in (9) and its canonical momentum as

(50) | |||||

(51) |

with the time independent creation and annihilation operators, such that annihilates the initial Fock vacuum state. Using the initial conditions on the mode functions, the Heisenberg field operators are written as

(52) | |||||

(53) | |||||

(54) |

with the time evolution operator with the boundary condition . The Heisenberg operators are related to by a Bogoliubov (canonical) transformation (see reference[4] for details).

The particle number defined with respect to the initial Fock vacuum state is defined in term of the dimensionless variables introduced above as

(55) |

It is this definition of particle number that will be used for the numerical study.

In order to define the particle number with respect to the adiabatic vacuum state we note that the mode equations (39,46) are those of harmonic oscillators with time dependent squared frequencies

(56) |

with for the unbroken symmetry case and for the broken symmetry case, respectively. When the frequencies are real, the adiabatic modes can be introduced in the following manner:

(57) | |||||

where now is a canonical operator that destroys the adiabatic vacuum state, and is related to by a Bogoliubov transformation. This expansion diagonalizes the instantaneous Hamiltonian in terms of the canonical operators . The adiabatic particle number is

(59) |

As mentioned above, the adiabatic particle number can only be defined when the frequencies are real. Thus, in the broken symmetry state they can only be defined for wave-vectors larger than the maximum unstable wave-vector, . These adiabatic modes and the corresponding adiabatic particle number have been used previously within the non-equilibrium context[9, 10, 11] and will be very useful in the analysis of the energy below. Both definitions coincide at because . Notice that due to the fact that we are choosing zero initial temperature. (We considered a non-zero initial temperature in refs.[4, 12]).

### ii.5 Energy and Pressure

The energy momentum tensor for this theory is given by

(60) |

Using the large factorization (4-7, 9) we find the energy density operator to be

(61) | |||||

(62) |

Taking the expectation value in the initial state and the infinite volume limit, defining , and recalling that the tadpole condition requires that the expectation value of vanishes, we find the expectation value of the energy to be

(63) |

It is now straightforward to prove that this bare energy is conserved using the equations of motion (19)-(21). It is important to account for the last term when taking the time derivative because this term cancels a similar term in the time derivative of .

Since we consider translationally as well as rotationally invariant states, the expectation value of takes the fluid form

(64) |

We want to emphasize that the full evolution of the zero mode plus the back-reaction with quantum fluctuations conserves energy. Such is obviously not the case in most treatments of reheating in the literature in which back-reaction effects on the zero mode are neglected. Without energy conservation, the quantum fluctuations grow without bound. In cosmological scenarios energy is not conserved but its time dependence is not arbitrary; in a fixed space-time background metric it is determined by the covariant conservation of the energy momentum tensor. There again only a full account of the quantum back-reaction will maintain covariant conservation of the energy momentum tensor.

We can write the integral in eq.(63) as

(65) | |||

(66) |

where is a spatial upper momentum cutoff, taken to infinity after renormalization. In the broken symmetry case, is the contribution to the energy-momentum tensor from the unstable modes with negative squared frequencies, and is the adiabatic particle number given by eq.(59). For the unbroken symmetry case and .

This representation is particularly useful in dealing with renormalization of the energy. Since the energy is conserved, a subtraction at suffices to render it finite in terms of the renormalized coupling and mass. Using energy conservation and the renormalization conditions in the large limit, we find that the contribution is finite. This also follows from the asymptotic behaviors (28).

In terms of dimensionless quantities, the renormalized energy density is, after taking :

(67) | |||||

(68) | |||||

(69) |

where the lower sign and apply to the broken symmetry case while the upper sign and correspond to the unbroken symmetry case. The constant is chosen such that coincides with the classical energy for the zero mode. The quantity is identified as the effective (dimensionless) mass for the “pions”.

We find using the renormalized eqs. (38), (39), (44), (45) and (46), that the renormalized energy is indeed conserved both for unbroken and for broken symmetry.

The pressure is obtained from the spatial components of the energy momentum tensor (see eq.(64)) and we find the expectation value of the pressure density to be given by

(70) |

Using the large- behavior of the mode functions (28), we find that aside from the time independent divergence that is present also in the energy the pressure needs an extra subtraction compared with the energy. Such a term corresponds to an additive renormalization of the energy-momentum tensor of the form

(71) |

with a (divergent) constant[17]. Performing the integrals with a spatial ultraviolet cutoff, and in terms of the renormalization scale introduced before, we find

(72) |

In terms of dimensionless quantities and after subtracting a time independent quartic divergence, we finally find setting ,

(74) | |||||

(75) |

At this stage we can recognize why the effective potential is an irrelevant quantity to study the dynamics.

The sum of terms without in (67) for are identified with the effective potential in this approximation for a time independent . These arise from the “zero point” energy of the oscillators with time dependent frequency in (65).

In the broken symmetry case the term describes the dynamics of the spinodal instabilities[12] since the mode functions will grow in time. Ignoring these instabilities and setting as is done in a calculation of the effective potential results in an imaginary part. In the unbroken symmetry () case the sum of terms without give the effective potential in the large limit, but the term describes the profuse particle production via parametric amplification, the mode functions in the unstable bands give a contribution to this term that eventually becomes non-perturbatively large and comparable to the tree level terms as will be described in detail below. Clearly both in the broken and unbroken symmetry cases the effective potential misses all of the interesting non-perturbative dynamics, that is the exponential growth of quantum fluctuations and the ensuing particle production, either associated with unstable bands in the unbroken symmetry case or spinodal instabilities in the broken symmetry phase.

The expression for the renormalized energy density given by (67-69) differs from the effective potential in several fundamental aspects: i) it is always real as opposed to the effective potential that becomes complex in the spinodal region, ii) it accounts for particle production and time dependent phenomena.

The effective potential is a useless tool to study the dynamics precisely because it misses the profuse particle production associated with these dynamical, non-equilibrium and non-perturbative processes.

## Iii The Unbroken Symmetry Case

### iii.1 Analytic Results

In this section we turn to the analytic treatment of equations (38,39,44) in the unbroken symmetry case. Our approximations will only be valid in the weak coupling regime and for times small enough so that the quantum fluctuations, i.e. are not large compared to the “tree level” quantities. We will see that this encompasses the times in which most of the interesting physics occurs.

Since , the back-reaction term is expected to be small for small during an interval say . This time , to be determined below, determines the relevant time scale for preheating and will be called the preheating time.

During the interval of time in which the back-reaction term can be neglected, we can solve eq.(38) in terms of elliptic functions, with the result:

(77) | |||||

(78) |

where cn stands for the Jacobi cosine. Notice that has period , where is the complete elliptic integral of first kind. In addition we note that since

(79) |

if we neglect the back-reaction in the mode equations, the ‘potential’ is periodic with period . Inserting this form for in eq.(39) and neglecting yields

(80) |

This is the Lamé equation for a particular value of the coefficients that make it solvable in terms of Jacobi functions [18]. We summarize here the results for the mode functions. The derivations are given in the Appendix.

Since the coefficients of eq.(80) are periodic with period , the mode functions can be chosen to be quasi-periodic (Floquet type) with quasi-period .

(81) |

where the Floquet indices are independent of . In the allowed zones, is a real number and the functions are bounded with a constant maximum amplitude. In the forbidden zones has a non-zero imaginary part and the amplitude of the solutions either grows or decreases exponentially.

Obviously, the Floquet modes cannot obey in general the initial conditions given by (20) and the proper mode functions with these initial conditions will be obtained as linear combinations of the Floquet solutions. We normalize the Floquet solutions as

(82) |

We choose and as an independent set of solutions of the second order differential equation (80). It follows from eq.(81) that has as its Floquet index.

We can now express the modes with the proper boundary conditions (see 20) as the following linear combinations of and

(83) |

where is the Wronskian of the two Floquet solutions

(84) |

Eq.(80) corresponds to a Schrödinger-like equation with a one-zone potential[19]. We find two allowed bands and two forbidden bands. The allowed bands correspond to

(85) |

and the forbidden bands to

(86) |

The last forbidden band is for positive and hence will contribute to the exponential growth of the fluctuation function .

The mode functions can be written explicitly in terms of Jacobi -functions for each band. We find for the forbidden band,

(87) |

where is a function of in the forbidden band defined by

(88) |

and is the Jacobi zeta function [20]. It can be expanded in series as follows