8 we first discuss perturbative thermalization. Then we summarize the results of recent non-perturbative studies of inflaton decay. We begin with perturbative considera- tions. In a full thermal equilibrium the energy density ρ and the number density n of relativistic particles scale as: ρ ∼ T 4 and n ∼ T $^{3}$, where T is the temperature of the thermal bath. Thus, in full equilibrium the average particle energy is given by: 〈 E 〉 eq = ( ρ/n ), which obeys the scaling, 〈 E 〉 eq ∼ ρ 1 / 4 ∼ T . On the other hand, if the inflaton decays perturba- tively, then right after the inflaton decay has completed, the energy density of the universe is given by: particles increases very rapidly [52] 3 As a result, full thermal equilibrium will be established shortly after that. Based on the above analysis, one can use the rate for the above inelastic scatterings as the thermalization rate Γ$_{th}$ of the universe . This rate at the time when the infla- ton decay completes can be found by using Eqs. (49,50): For typical values of α ∼ 10 − $^{2}$− 10 − 1 and m ≤ 10 − $^{5}$m$_{pl}$ , we find Γ$_{th}$ ≥ Γ. Therefore the universe reaches full thermal equilibrium immediately after the completion of perturbative inflaton decay. (where m is the inflaton mass). Then, from conservation of energy, the number density of decay products is found to be Hence, perturbative inflaton decay results in a dilute plasma that contains a small number of very energetic particles. Reaching full equilibrium requires re-distribution of the energy among different particles, kinetic equilibrium , as well as increasing the total number of particles, chem- ical equilibrium . Therefore both number-conserving and number-violating reactions must be involved. The most important processes for kinetic equilibration are 2 → 2 scatterings with gauge boson exchange in the t -channel. (Scalar exchange in t -channel diagrams are usually suppressed, also vertices that arise from a Yukawa coupling are helicity suppressed.) The cross-section for these scatterings is $^{σ$_{2}$}$→ 2 ∼ α $^{2}$| t | − $^{1}$, where α is a gauge fine structure constant and the variable t is related to the exchanged energy Δ E and momentum, Δ p through t = Δ E 2 − | Δ p | $^{2}$. Due to an infrared singularity, these scatterings are very efficient even in a dilute plasma [47, 48]. Chemical equilibrium is achieved by changing the num- ber of particles in the reheat plasma. From (50) it follows that in order to reach full equilibrium the total number of particles must increase by a factor of n$_{eq}$/n , where n ≈ ρ/m and the equilibrium value is: n$_{eq}$ ∼ ρ 3 / $^{4}$. This can be a very large number, e.g. n$_{eq}$/n ∼ O (10 $^{3}$). It was recognized in [47, 49] (see also [50, 51]) that the most rel- evant processes are 2 → 3 scatterings with gauge-boson exchange in the t − channel. The cross-section for emit- ting a gauge boson whose energy is E ∼ (Γ m$_{pl}$ ) 1 / 2 (where as in earlier sections of this review Γ is the inflaton decay rate), from the scattering of two fermions (up to a loga- rithmic “bremsstrahlung” factor) is $^{σ$_{2}$}$→ 3 ∼ α $^{3}$(Γ m$_{pl}$ $_{)}$− $^{1}$. When these scattering become efficient, the number of The purely perturbative considerations of the previous subsection are subject to the same criticisms as the orig- inal perturbative analysis of the initial stages of reheat- ing. Hence, we must turn to non-perturbative analyses. Some analytical approaches were pioneered in [6] and [7]. Numerical studies, however, have proved more powerful. Since the occupation numbers of the excited modes are typically very high after the initial stages of preheating, a classical field theory analysis should be justified. Ini- tial numerical studies were pioneered in [31–34, 53]. Two numerical package to perform such simulations are pub- licly available [54, 55]. Detailed numerical simulations of tachyonic preheating are given in [40]. Here, we will focus on numerical studies of preheating in models with narrow or broad parametric resonance [11]. The resonant phase of the reheating process described in Section 3 produces either (in the case of narrow res- onance) field fluctuations in a narrow interval about the resonant frequency (which is set by the mass of the in- flaton field), or else (in the case of broad or tachyonic resonance) field fluctuations at all wavenumbers smaller than the critical value given in (23), whose magnitude is set by the inflaton mass and amplitude. Once the occupation numbers of the resonant modes become sufficiently large, re-scattering of the fluctuations begins. As first studied in [31–34] this terminates the phase of exponential growth of the occupation numbers. In the case of narrow resonance, new peaks in the spec- trum of the number density n ( k ) develop at harmonic frequencies. Soon, the spectrum shows excitations in a continuum band which reaches to k = 0 and has an ultra- violet (UV) cutoff whose value increases with increasing time. As studied in detail in [11], the evolution of the field fluctuations evolves to a regime of turbulent scaling 3 Decay processes which were considered in [51], are helpful, but in general they cannot increase the number of particles to the required level. B. Non-Perturbative Considerations Γ$_{th}$ ∼ α 3 ( m$_{pl}$ m ) Γ . (51) ρ ≈ 3 (Γ m$_{pl}$ ) 2 , and 〈 E 〉 ≈ m ≫ ρ 1 / $^{4}$, (49) n ≈ ( ρ m ) ≪ ρ 3 / $^{4}$. (50)