Lectures on Reheating after Inflation

This lecture note provides a generic introduction to reheating after inflation, covering theoretical, phenomenological and observational aspects of the process.

Acknowledgments I would like to thank my PhD supervisors, Mustafa Amin, who has taught me everything I know about reheating, and Anthony Challinor, who introduced me to the field of Cosmology. Without them, the writing of these notes would have not been possible. I am also immensely grateful to Eiichiro Komatsu for giving me the opportunity to continue my research in the area of reheating and to give short lecture series on it.

Avant propos
Over the past decades, the understanding of our cosmic history has improved spectacularly. Precise measurements of the temperature anisotropies in the cosmic microwave background (CMB) [1,2] have revealed a homogeneous and isotropic universe on large scales with tiny (if any) spatial curvature, Ω K , and nearly scale-invariant primordial curvature perturbations. The observations provide compelling evidence for an inflationary phase in the early universe [3][4][5][6][7][8], implying a connection between gravity and quantum mechanics. Measurements of the relative abundances of light-elements are also in excellent agreement with our expectations [9]. The predictions of big-bang nucleosynthesis (BBN), based on the well-understood physics of nuclear reactions, point towards a hot and dense universe, in local thermal equilibrium at late times. They tell us that the first light-elements started forming at cosmic time t BBN 1 s at equilibrium temperature T BBN ∼ E BBN ∼ 1 MeV. Connecting these two remarkable epochs, however, could be challenging. Since the energy scale at the end of inflation can be as high as E inf ∼ 10 16 GeV, with the duration of inflation corresponding to ∆t inf 10 −36 s, there is a huge range of energy (and time) scales which is poorly understood and poorly observationally constrained. Current cosmological experiments cannot probe the period between inflation and BBN. This is because typically the effects from the interesting post-inflationary dynamics are on subhorizon scales due to causality and are washed out by the later non-linear evolution of structure. The thermal state of the universe, required for BBN, also hides information about earlier times. Collider experiments can shed light on some intermediate-energy phenomena, e.g., the electro-weak symmetry breaking and QCD phase transitions, E EW ∼ 10 2 GeV and E QCD ∼ 10 2 MeV, respectively, but they will not be able to cover the entire energy range in the foreseeable future.
Nevertheless, it is crucial that we try to understand the period between inflation and BBN better, for both theoretical and observational reasons. In the standard lore, the universe at the end of inflation is cold and dark, virtually empty of particles and dominated by the approximately homogeneous inflaton field. The energy of the inflaton that drove inflation, must be somehow transferred to other species of matter, eventually populating all relevant degrees of freedom of the Standard Model, leaving a hot, thermal, radiationdominated universe, setting the scene for BBN. Importantly, this reheating process explains not only the cosmic origin of the matter that we are made of, but it also accounts for the production of cosmic relics such as photons and neutrinos, and perhaps dark matter and gravitational waves, as well as the generation of the observed matter-antimatter asymmetry in our universe (baryogenesis). Any unified theory of high-energy physics must include a complete understanding of inflation, reheating and the later evolution of the universe.
Given the current and planned advances in observational cosmology and the improvement of constraints on inflation, reheating will be an integral part of research in the coming years. Arguably, one of the most important observational implications of the postinflationary dynamics is its effect on the expansion history of the universe between inflation and BBN. It determines how we map perturbation modes from their exiting the Hubble horizon during inflation to horizon re-entry at late times. Thereby, the poorly constrained and understood post-inflationary expansion history influences directly the predictions for cosmological observables of specific models of inflation. For instance, it leads to significant uncertainties in the predictions for the spectral index, n s , and the tensor-to-scalar ratio, r, in different models. It is critical that we understand these uncertainties better, if we wish to narrow the range of the observationally-allowed models of inflation. Works on reheating have also shown the possibility of formation of relics such as solitons and cosmic defects, helping us further constrain the variety of scenarios.
The initial stage of reheating, also known as preheating, can involve highly nonperturbative processes, during which the universe gets populated via parametric resonances. They cannot be described with the usual perturbative expansions in coupling constants, even in cases with weak couplings [10][11][12][13][14][15][16][17]. Such resonances arise as the inflaton condensate (or generally any light scalar, that has attained a non-zero vacuum expectation value during inflation) begins to oscillate about the minimum of its potential, soon after inflation. The oscillations induce an effective time-dependence in the couplings of the inflaton to the other species of matter. While the background inflaton field dominates the energy budget of the system, the evolution of the remaining fields it is coupled to can be linearised. As the effective frequencies of the individual Fourier modes of the daughter fields change nonadiabatically every time the inflaton crosses the origin, we observe 'explosive' (or resonant) particle production. This can be quite efficient, since it involves the collective decay of many inflatons from the condensate. When the energy of the newly-populated degrees of freedom becomes comparable to the background, back-reaction effects become important. Typically, the condensate fragments and the subsequent evolution is non-linear. It can be studied in the classical approximation, using classical lattice simulations, since all relevant modes have large occupancies and hence quantum effects are negligible. 1 The non-linear dynamics can lead to many interesting phenomena, e.g., the production of solitons that can delay the thermalisation required as an initial condition for BBN, field configurations evolving self-similarly in a turbulent manner [18,19], non-thermal phase transitions and the production of cosmic defects [20][21][22].
These notes are meant to serve as a generic introduction to the field of reheating after inflation, starting with a brief summary of the inflationary paradigm (Section 4), followed by a review of the different preheating mechanisms (Section 5) and the ensuing non-linear evolution (Section 6), and finally considering some high-energy physics (Section 7) and observational (Section 8) aspects of reheating. 1 However, since the universe at the beginning of BBN is in local thermal equilibrium, late-stage reheating analysis should eventually include a full quantum mechanical computation of the approach of all relevant degrees of freedom to states with maximal local entropy such as Bose-Einstein or Fermi-Dirac distributions [18,19].

Notation and conventions
We use natural units in which = c = k B = 0 = 1. In these units, the reduced Planck mass is given by m Pl = 1/ √ 8πG.
Greek indices µ, ν and so on go over the four space-time coordinates x µ = [x 0 , x 1 , x 2 , x 3 ] T with x 0 for the time coordinate.
Latin labels i, j, k and so on go over the three spatial coordinates.
Spatial vectors are written in boldface.
Summation over repeated indeces is assumed unless otherwise stated.
The Ricci tensor, defined in terms of the Christoffel symbols, is 1) and the Ricci scalar is R = g µν R µν .
The spatial Fourier transform of a field f (x) is f k = f (x)e −ik·x d 3 x/(2π) 3 and the inverse transform is f (x) = f k e ik·x d 3 k.
4 Inflation and initial conditions for reheating 'With the new cosmology the universe must have been started off in some very simple way. What, then, becomes of the initial conditions required by dynamical theory? Plainly there cannot be any, or they must be trivial. We are left in a situation which would be untenable with the old mechanics. If the universe were simply the motion which follows from a given scheme of equations of motion with trivial initial conditions, it could not contain the complexity we observe. Quantum mechanics provides an escape from the difficulty. It enables us to ascribe the complexity to the quantum jumps, lying outside the scheme of equations of motion. The quantum jumps now form the uncalculable part of natural phenomena, to replace the initial conditions of the old mechanistic view.' P. A. M. Dirac (1939)

Standard cosmology, its puzzles and why we need inflation
Standard cosmology is based on the empirical observation that the universe is homogeneous and isotropic on large scales [1]. In the context of General Relativity, it means that the space-time metric takes the Friedmann-Robertson-Walker (FRW) form where t is the cosmic time and a(t) is the Robertson-Walker scale factor. The term in brackets represents the line element of the three-dimensional homogeneous and isotropic space. For positive, zero and negative K this hypersurface can be considered as a 3-dimensional sphere embedded in a 4-dimensional Euclidean space, a 3-dimensional Euclidean space and a 3-dimensional hypersphere embedded in a 4-dimensional pseudo-Euclidean space, respectively [52]. The positive, zero and negative cases are better known as the closed, flat and open universes, respectively. The evolution of a(t) is determined by the Einstein equations Here Λ is the cosmological constant, introduced by Einstein to make the universe static. Henceforth, we set Λ = 0 2 . Near the origin of locally Cartesian co-moving coordinates, the components of the energy-momentum tensor, T µν , in a homogeneous universe are functions of t only. Isotropy also imposes the additional constraints T i0 = 0 and T ij ∝ δ ij . 2 Observations favor a small, but non-zero value of Λ, unjustifiable by Quantum Field Theory (QFT) if interpreted as the energy of the vacuum. This drawback of standard cosmology and QFT is known as the Cosmological constant problem. The contribution of the Λ term to the energy budget of the universe (known as dark energy) becomes significant only at very late times, at a redshift of about 1, so it is safe to ignore Λ at earlier epochs.
Conventionally, the energy density, ρ, and the pressure, p, of the perfect fluid filling the homogeneous and isotropic universe are defined locally as T 00 = ρ(t) , T ij = −a(t) −2 δ ij p(t) .

(4.3)
Then only the 00 and the ii Einstein equations, eq. (4.2), do not vanish 4) and are known as the Friedmann and Raychaudhuri equations, respectively. The Einstein equations, eq. (4.2), also imply the conservation of the energy-momentum tensor T νµ ;µ = 0. Due to isotropy the momentum-conservation law T iµ ;µ = 0 is automatically satisfied. The energy-conservation T 0µ ;µ = 0 yieldsρ + 3H(ρ + p) = 0 . In standard cosmology the typical sources of gravity are non-relativistic matter (dust) and relativistic matter (radiation). If one of these components is dominant, then for dust : a ∝ t 2/3 , ρ ∝ a −3 , w = 0 , for radiation : a ∝ t 1/2 , ρ ∝ a −4 , w = 1/3 . The second column also applies for individual species, even if subdominant, provided they are self-interacting only. Note that the universe always decelerates,ä < 0. The Friedmann equation, eq. (4.4), can be rewritten as where the energy density makes up a fraction Ω = ρ/ρ c of the critical energy density, ρ c = 3m 2 Pl H 2 . Similarly, one often writes the spatial curvature term as Ω K = −K/(a 2 H 2 ). Cosmological observations give very tight constraints on this quantity, consistent with zero. The 95% limit, from the most recent measurement of the anisotropies in the CMB [1], on the spatial curvature today is |Ω K,0 | < 0.005. This small value leads to one of the fine-tuning problems in standard cosmology. Sinceä < 0, aH =ȧ increases when going backwards in time. Hence, Ω K becomes even smaller at earlier times, or in other words, the energy density tends to the critical one, Ω → 1, with unnaturally high precision. To get some rough idea about the degree of fine-tuning in the initial condition for Ω at some early time, t early , required by current measurements, consider the ratio a 0 H 0 a eq H eq 2 a eq H eq a early H early 2 ∼ 1 + z eq (1 + z early ) 2 . (4.9) We assume the universe to be radiation dominated at early times, t early < t < t eq , and matter dominated after the epoch of radiation-matter equality t eq < t < t 0 . For t early ∼ t BBN , |Ω BBN − 1| 10 −17 , assuming T BBN ∼ 1 MeV and z eq ∼ 10 3 . At the GUT epoch, Unless the initial conditions are set very precisely, the universe either collapses too quickly or expands too fast before large-scale structure can form. This is known as the flatness problem.
The deceleration of the scale factor in standard cosmology also leads to contradictions with measurements of the causal structure of the observable universe. The physical length, l phys = al, of a given co-moving length-scale l, increases as the universe expands. On the other hand, its ratio with the Hubble radius l phys /H −1 =ȧl decreases with time, sincë a < 0. Hence, any co-moving length-scale becomes much greater than the Hubble scale at sufficiently early times. This implies that the causally-connected region in the universe today should lie deep inside the Hubble volume. We can easily calculate the expected size of this region at the epoch of recombination. Ignoring spatial curvature, the physical size of the region at recombination is equal to the particle horizon d rec = a rec trec t early dt/a(t). For a ∼ t n , 0 < n < 1 and t rec t early the upper bound of the integral determines its value. For dust or radiation or any mixture of the two, d rec ∼ H −1 rec . We also have d rec = a rec r rec ∆θ, where the co-moving distance to the CMB is r rec = t 0 trec dt/a(t) ∼ 1/(a 0 H 0 ) and ∆θ is the angular size of the causal region. Assuming matter domination between recombination and today, ∆θ ∼ (1 + z rec ) −1/2 , which corresponds to about 1 o for z rec ∼ 10 3 . This is in conflict with observations of the microwave sky, showing the same temperature to high accuracy in all directions [1]. Within standard cosmology this isotropy of the CMB cannot be accounted for, since there is no way for points separated by more than a degree to be in thermal equilibrium (and in causal contact) before the epoch of last-scattering. This constitutes the horizon problem.
The tiny anisotropies measured in the CMB have a nearly scale-invariant power-spectrum, even on large, causally disconnected at the time of last-scattering, scales [1], see Fig. 1. They reflect small variations in the matter density at the epoch of recombination. Later on, these density fluctuations act as seeds for the formation of large-scale structure. Hence, standard cosmology also fails to explain the deviations from the FRW metric. 3 3 The appearance of well-defined acoustic peaks on small, causally connected at the time of last-scattering scales, is also a compelling evidence for deviations from standard cosmology. Only if all temperature perturbation modes of a given k, irrespective of the direction of the wavevector k, oscillate in phase can we observe these prominent acoustic peaks. This phase coherence suggests that fluctuations were seeded in the early universe, before horizon entry (the moment they started oscillating). Alternative models in which fluctuations are sourced continuously during the radiation and matter domination eras, e.g., by a network of cosmic defects, do not give rise to phase coherence. The role of such sources for structure formation can be only subdominant. Figure 1. A full-sky map of the CMB temperature (top) and its power spectrum (bottom) [1] with the mean temperatureT CMB = 2.7255 K subtracted. The angular size of the particle horizon at recombination, ∆θ ∼ 1 o , corresponds roughly to l ≈ 200.
High-energy physics theories take the view that the physical laws look simpler at higher energies. This implies that gauge symmetries inevitably get broken during the evolution of the early universe, leaving behind them topological defects such as monopoles. The density of monopoles is bounded from below, due to the existence of a maximal correlation length determined by the causal length, i.e., the particle horizon, during the phase transition. This presents a serious problem for standard cosmology. If these relics do not annihilate efficiently, their abundance on sub-horizon scales can become large after the phase transition. Furthermore, these massive relics behave as dust. Their energy density can become the dominant component at dangerously early times, e.g., before or around BBN, since it decays more slowly with a(t) than that of radiation. This is known as the monopole problem 4 .
All of the above problems can be shown to have the same origin -the expansion with time of the co-moving Hubble sphere, (aH) −1 =ȧ −1 , following fromä < 0. One can resolve all of these puzzles by postulating the existence of an earlier stage of inflation, during which the universe undergoes accelerated expansion and the co-moving Hubble sphere shrinks [3][4][5]. The most common expansion history of inflation is the quasi-exponential one, i.e., a ∼ exp(Ht), with H varying very slowly with time. Another possibility is power-law inflation, a ∼ t n , n > 1. For n 1, H again varies very slowly with time. The measured tiny value of |Ω K,0 | is in fact a prediction of inflation. If we assume that |Ω K,init | at the beginning of inflation is of order unity then from eq. (4.8) follows For a constant expansion rate during inflation, H inf , the number of e-folds of expansion until the end of inflation, N = ln(a end /a init ), is bounded to be Making the tentative assumption of t end ∼ t early , see eq. (4.9), i.e., of the universe becoming radiation dominated immediately after the end of inflation, we find that for t early ∼ t BBN we need at least about 20 e-folds to resolve the flatness problem, whereas for t early ∼ t GUT , N has to be over 60.
To account for the isotropy of the CMB, we need to make sure that the co-moving particle horizon at recombination is greater than the co-moving distance photons travel after recombination until today, i.e., Note that dt/a(t) = d ln a/ȧ. Hence, each integral is dominated by the smallestȧ, i.e., the moment when the co-moving Hubble sphere is the largest. It implies for the right-hand side of eq. (4.12) a value ∼ 1/(a 0 H 0 ) as already shown above and ∼ 1/(a init H inf ) for the left-hand side. Thus, the condition in eq. (4.11) also applies to the horizon problem. The scale-invariance of the power-spectrum of the small density fluctuations, imprinted on the CMB as tiny anisotropies, is also a consequence of the accelerated expansion during inflation. At the beginning of inflation, small-scale perturbations, lying deep inside the Hubble radius, are generated by Minkowski space-time quantum fluctuations (since the space-time curvature can be neglected). As the universe undergoes accelerated expansion, perturbations of fixed co-moving wavelength cross outside the shrinking co-moving Hubble sphere. As they become superhorizon, the Hubble friction term starts to dominate and they become over-damped. Since H is approximately constant during inflation and is the only scale determining the evolution of perturbations after Hubble exit, the perturbations which leave the Hubble sphere during inflation have an almost scale-invariant power-spectrum. As the co-moving Hubble radius begins to increase after inflation, perturbations of a given co-moving wavelength start to re-enter the horizon, accounting for the observed approximate scale-invariance of density perturbations over a range of different scales. By letting the Hubble sphere shrink during inflation and then begin to expand as the universe becomes radiation and then matter dominated we provide a causal mechanism for producing seemingly-acausal correlations in the density perturbations. 5 Basically, inflation opens up the past light cones for fundamental observers (those who are stationary with respect to the cosmic grid).
Inflation also provides a straightforward solution to the monopole problem. If the phase transition occurs during or before inflation, we need to make sure that there is sufficient number of e-folds of accelerated expansion after the transition to dilute the concentration of the relics. If the gauge symmetry is broken after inflation, the correlation length at the phase transition, equal to the particle horizon, is substantially increased in comparison to standard cosmology and we can again put a lower bound on the duration of inflation. For instance, let us consider the generation of monopoles after spontaneous symmetry breaking at the GUT scale, E GUT = 10 16 GeV, immediately after the end of inflation. The number density of monopoles is determined by the correlation length of the scalar fields, which is set by the particle horizon and is ∼ a end /(a init H inf ) = e N /H inf . On the other hand, Assuming the universe becomes radiation dominated and reaches local thermal equilibrium soon after the formation of monopoles, the number density of photons at that time is ∼ E 3 GUT . Ignoring subsequent processes that can change the number of photons and annihilation of monopoles, the ratio of the number densities of monopoles and photons remains constant, since each ∝ a −3 . Thus, n mon 0 /n γ 0 ∼ e 3N E 3 GUT /m 3 Pl = 10 −9 e 3N . For less than 10 −39 monopoles per photon [52], as suggested by terrestrial experiments, N > 23. An identical bound is obtained if GUT scale monopoles are generated during or before inflation. A possibility that cannot be resolved by inflation is the production of massive particles during the reheating process after inflation. If the reheating temperature is low enough, thermal particle production of dangerous massive relics, that can ruin the successful BBN in standard cosmology, can be evaded. On the other hand, non-thermal production has to be dealt with on a case-by-case basis.

Dynamics of inflation and setting the scene for reheating
The semi-classical theory of inflation provides not only a solution to the horizon, flatness and monopole problems, but more importantly predicts the generation of density perturbations exhibiting a nearly scale-invariant power-spectrum. These act as seeds for the large-scale structure in the late universe.
The inflationary paradigm can be interpreted in at least two different ways. We could think of it as an approximate description (some sort of a parametrisation), which does not capture the actual physical laws, due to its semi-classical nature. We could also argue that the universe genuinely underwent a stage of accelerated expansion, driven by a scalar condensate whose origin can be traced back to any of the high-energy models going beyond the Standard Model of Particle Physics. This way of thinking makes inflation the link between quantum gravity or/and extensions of the Standard Model, and the well-understood physics of BBN in standard cosmology.
In this note, we consider the most common (and consistent with observations [2]) models of inflation, in which a single scalar field φ, called the inflaton, sources the accelerated expansion of the universe, with action We limit ourselves to models, minimally coupled to gravity, with canonical kinetic terms 6 . The matter action term, S matter , contains the entire information regarding the other constituents of the matter sector, including the Standard Model Lagrangian as well as the terms describing the couplings of the inflaton to other fields.

Homogeneous dynamics of inflation
Isotropy and homogeneity require that the dominant component of the scalar field depends on t only,φ(t). This scalar condensate provides the classical background configuration during inflation (and the initial stages of reheating). The energy density and the pressure of the isotropic and homogeneous scalar fluid are simply (4.14) To have acceleration,ä > 0, the Raychaudhuri equation, eq. (4.4), demands ρφ + 3pφ < 0.

Slow-roll inflation
Having derived the equations of motion and shown that accelerated expansion is possible, we need to find what conditions V has to satisfy to have enough number of e-folds of inflation to solve the horizon and flatness problems. Note thatä = a(Ḣ + H 2 ), hence has to hold to haveä > 0. Inflation ends when H = 1, corresponding toä = 0. Since H −1 is the characteristic time-scale for one e-fold of expansion (recall dN = d ln a = Hdt), known as the Hubble time, H = −d ln H/dN < 1 implies that the time-scale over which 6 There are examples of non-minimal and non-canonical models in which a conformal transformation or a field redefinition can reduce the action to the form given in eq. (4.13).
the fractional decrease in H is significant is greater than a Hubble time. Or in other words the rate of decrease of H must be slower than the rate of expansion of the universe in order to have inflation. To achieve sufficiently many e-folds (at least 40 up to 60, see Section 4.1) of inflation we need H to be much less than 1 for a long enough period, implying that for most of the time H 1. The parameter quantifying the rate of change of H is For H to increase slowly, over many Hubble times, |η H | has to be less than unity. But because of the large number of e-folds required by observations, |η H | 1 has to hold 7 . All of this basically means that H ≈ const for most of inflation, and that the scale factor increases quasi-exponentially. That is why this period is also called quasi-de Sitter expansion. Current observational constraints are roughly H < 0.01 and η H ≈ 0.03 ± 0.01 [2] and are in support of this picture.
Using the Friedmann equation, eq. (4.4), the energy conservation equation, eq. (4.5), and the expressions for the energy density and pressure of the inflaton condensate, eq.
(4.14), we find that H = 3φ 2 /(2ρφ). Thus, H 1 implies a negligible contribution of the kinetic energy density to the total energy density of the condensate, which also means thaṫ during inflation, consistent with our comments under eq. (4.14). Hence, V has to be very flat, forφ to roll sufficiently slowly. This is called slow-roll inflation. For slow-roll inflation to last long enough, we need the kinetic energy density to remain small. This means that the fractional change inφ during one expansion time H −1 , |φ/φ|H −1 , has to be much less than 1. It implies, given H 1, that |η H | 1, since one can show that η H = 2 H + 2φ/(φH).
To see that within the slow-roll approximation, the expansion during inflation can be 7 If ηH ∼ O(1), inflation still takes place for roughly ∆N ∼ ln(1/ H ). If H is several orders of magnitude less than unity, the number of e-folds of inflation is much less than the one required by observations. exponentially large consider 8 where we assume 0 < V (φ end ) < V (φ init ) and ∆φ ≡φ init −φ end > 0, implying a positive argument in the last exponential. The slow-roll condition V 1 leads to N ∆φ/m Pl . If the value of the inflaton changes by ∼ m Pl , we definitely get a huge number of e-foldings. Note that such large field values do not mean that the quantum nature of gravity becomes important. For this to happen the energy density of the condensate must be ∼ m 4 Pl . This can be easily avoided, even forφ m Pl , if V (φ) is proportional to a sufficiently small coupling constant. None of the potential slow-roll parameters, eq. (4.18), depends on it.
We can find the approximate trajectory of the inflaton during inflation. Since, during slow-roll H ≈ V 1,φ ≈ −m Pl ∂φV / √ 3V . As a test case, we consider the simplest form for the inflaton potential, i.e., V (φ) = m 2 φ 2 /2. In general, all monomial potentials satisfy the slow-roll conditions, eq. (4.18), for some ∆φ m Pl . These models belong to the class of Chaotic inflation [53,54]. It encompasses all models having V that supports slow-roll inflation for ∆φ ∼ O(m Pl ) or smaller. In these lecture notes, we will concentrate on Chaotic inflation models 9 . For the quadratic potential the slow-roll trajectory and the expansion law take the simple approximate analytic form (t init = 0) (4.20) It breaks down towards the end of inflation,φ end ≈ m Pl for which the slow-roll conditions in eq. (4.18) are violated. This solution is also known as the attractor solution, since one can show that for a broad range ofφ init [55,56], even for largeφ init such that H (t init ) > 1, the field velocity decays very rapidly andφ(t) and a(t) quickly approach eq. (4.20). This goes to show how broad the set of initial conditions is that can lead to an inflationary stage in chaotic models. Speaking of initial conditions, the term 'chaotic' derives from the possibility of having initially a scalar field varying randomly with position, i.e., having almost arbitrary initial conditions for the inflaton, and still getting slow-roll inflation after that. Even if the value of the inflaton varies from one spatial region to another, there should be a patch of space in which the inflaton looks uniform enough and has a value for which the slow-roll conditions in eq. (4.18) are satisfied, e.g., a super-Planckian value for monomial potentials [55]. One 8 The total number of e-folds of inflation is defined as the first integral with H (tinit) = H (t end ) = 1 and H (t) < 1 for tinit < t < t end . 9 Examples of single-field models that do not belong to Chaotic inflation include Small-field models in which necessarily ∆φ m Pl , models in which higher-order kinetic terms or higher-order curvature terms, instead of V , drive inflation and models in which phase transitions stop or trigger inflation, e.g., Old and New inflation, respectively.
can easily show that the initial physical size of the homogeneous patch, L init = a init l, has to obey H init L init φ init /m Pl for the gradients to be negligible. For monomial potentials, this implies that a sufficiently large uniform patch has to be super-Hubble initially 10 (and super-Planckian after imposing the condition of sub-Planckian V init ).

End of slow-roll inflation
For inflation to be successful it must feature a graceful exit into the deceleration stage of standard cosmology; otherwise the homogeneity and isotropy of the universe are destroyed.
A famous example of non-graceful exit is Alan Guth's Old inflation [3] in which the inflaton is initially trapped in a false vacuum. As the inflaton leaks through the potential barrier and forms bubbles of true vacuum, the energy released in the transition ends up concentrated within the bubble walls. If the bubbles are able to merge, a homogeneous and isotropic universe emerges. However, the bubbles never collide, since the background false-vacuum space in which they formed, never stops inflating. Hence, for an observer located inside a bubble the universe would appear highly anisotropic and inhomogeneous, since structure has to grow out of the energy concentrated in the bubble walls. The graceful exit problem is naturally avoided in Chaotic inflation. For power-law potentials, the homogeneous inflaton background simply begins to oscillate about the potential minimum. One can easily determine the oscillatory attractor solution for a quadratic minimum. We put φ = √ 6H(m Pl /m) cos θ andφ = √ 6Hm Pl sin θ to satisfy the first expression in eq. (4.15). After differentiating it with respect to time and using the second expression in eq. (4.15) we findḢ = −3H 2 sin 2 θ. This implies that H decays during the oscillatory stage as t −1 . Taking the time derivative of the new definition ofφ in terms of θ and using the definition ofφ in terms of θ we findθ = −m − (3/2)H sin(2θ). The second term on the right decays with time, so θ ≈ mt up to a constant for mt 1. We can use this result in the expression forḢ. After integration we find This decaying scalar field condensate provides the classical background during the reheating phase. Note that a ∝ t 2/3 up to subdominant decaying oscillating terms, which implies that the universe is in a dust-dominated state of expansion. Ultimately, the universe has to reheat itself to reach eventually a radiation-like state of expansion, with the inflaton energy transferred into radiation, baryons and leptons. Also note that even if the oscillating terms are very small, they can play an important role for the space-time curvature (neglecting the spatial curvature for simplicity) Before moving forward, we should point out that the state of expansion of a universe dominated by a homogeneous oscillating scalar, about the minimum of its potential, depends on the form of V . For simple power-laws, V ∝ |φ| 2n , where n need not be an integer, one can easily determine the temporal mean equation of state during the oscillatory phase. Ignoring expansion, since H decays with time after inflation, and multiplying byφ the second expression in eq. (4.15), we find φ ∂φV = − φφ = φ 2 . The angle brackets represent time averaging over many oscillations. The last equality follows from virialization, i.e., 0 = d(φφ)/dt = φ 2 + φφ . Thus, assuming pφ/ρφ ≈ pφ / ρφ , we find [59] w ≈ For quadratic potentials, n = 1, we have the expected w = 0 matter-like equation of state, whereas for quartic potentials, n = 2, we have a radiation-like equation of state, w = 1/3. When n ≤ 1/2, w ≤ −1/3, and eq. (4.6) tells us thatä > 0 -the universe inflates. However, consistency requires that the inflaton oscillates around a potential that has a non-singular first derivative at its minimum, for the equation of motion to be well-defined for all field values, implying the condition n > 1/2, i.e., the potential has to be steeper than linear at the minimum. Oscillations about such minima always lead to a decelerating stage of expansion with w > −1/3.

Cosmological perturbations from inflation
Having described what the homogeneous and isotropic universe looks like at the end of inflation, we now consider the small deviations from the FRW approximation. After all, these small departures enable us to distinguish between different models. As mentioned above, the tiny anisotropies measured in the CMB, as well as the seeds for Large Scale Structure (LSS) can be explained within the inflationary paradigm, as being microscopic quantum fluctuations, stretched to cosmic sizes during inflation. As we will see in later sections, they also laid down the seeds for particle production during reheating. To understand the initial conditions for this process, we need to consider their evolution during inflation. We expand the metric and the energy momentum tensor about their background values [60] g µν (x α ) =ḡ µν (t) + δg µν (x α ) , whereḡ µν (t) is given in eq. (4.1) and we set the spatial curvature, K, to zero, since during inflation it quickly becomes negligible. While we know the non-zero components ofT µν (t), see eq. (4.3), it is easier to work with a co-variant form for the tensor. Since the FRW universe is filled with a perfect homogeneous fluid, i.e., a fluid that looks the same in all directions for all co-moving observers at equal cosmic times, the rank-2 tensorT µν (t) has to be a linear combination ofρ(t)ḡ µν (t),p(t)ḡ µν (t),ρ(t)ū µ (t)ū ν (t) andp(t)ū µ (t)ū ν (t), wherē u µ (t) is the 4-velocity of a co-moving observer. The only linear combination that respects homogeneity and isotropy isT µν (t) = (ρ +p)ū µūν −ḡ µνp , where we used [1, 0, 0, 0] T for u µ (t). More generally, in an arbitrary gravitational field, a perfect fluid is a medium, with energy momentum tensor We now write the actual forms of the perturbations δg µν (x α ) and δT µν (x α ). We use the conformal time, dτ ≡ dt/a(t), which simplifies the background metric,ḡ µν (τ ) = a 2 (τ )η µν . The most general metric perturbations are Consistency requires the energy density, pressure and 4-velocity fields to be also perturbed Since the 4-velocity of an observer is normalized, i.e.,ḡ µνū µūν = 1 and g µν u µ u ν = 1 one can show to linear order that δu 0 = aϕ. In deriving this expression, we have used δg µν = −ḡ µα δg αβḡ βν , which also holds to linear order forḡ µνḡ να = δ α µ and g µν g να = δ α µ . This also implies to first order in perturbations. The perturbations in the energy momentum tensor then take the form (4.28) Decomposing perturbations into scalars, divergence-free vectors and traceless transverse tensors, also known as the scalar-vector-tensor decomposition, is very useful since the Einstein equations decouple the three kinds of modes to linear order. This is a consequence of the symmetries of the FRW background. The Einstein equations are also invariant under diffeomorphisms, i.e., space-time coordinate transformations (4.29) We assume ξ µ = [ξ 0 , ∂ i ξ + ξ ⊥ i ] T to be small, of the order of the metric and energy momentum tensor perturbations. Since the metric transforms as the metric perturbations at x α transform, to linear order, as 11 (4.31) 11 It is understood thatḡµν (x α ) =ḡ µν (x α ), i.e., eq. (4.29) yields δgµν ( Similarly, From this follows that perturbations depend on our choice of space-time coordinates, e.g., This issue can be resolved by working in diffeomorphism invariant quantities that take the same values for all choices of coordinates, i.e., in all gauges. For instance, the quantity is gauge invariant. It is known as the co-moving curvature perturbation. We can construct other gauge-invariant quantities, e.g., To make further progress, we need to specify the energy momentum tensor. In singlefield inflation, see eq. (4.13), the contribution from the matter action term to the energy momentum tensor vanishes at linear order. Hence, the scalar metric perturbations are coupled to linear order only to the perturbation in the inflaton φ(x α ) =φ(t) + δφ(x α ) . (4.36) The unperturbed energy momentum tensor is from which follows that Given the expressions for the background pressure and energy density in eq. (4.14), the scalar velocity perturbation can be found directly by using eq. (4.28) to be implying R = ψ +Hδφ/∂ τφ . We can calculate this quantity during inflation, working in the slow-roll approximation. The most suitable gauge for analytical analysis is the spatially-flat gauge defined as ψ = E = 0 (fixing the two scalar perturbations removes the gauge freedom due to ξ 0 and ξ ). This implies that the scalar modes in δg ij vanish. The second order terms in the action in eq. (4.13) which couple the non-zero scalar metric perturbations, δg 00 and δg 0i , to the inflaton perturbations, coming from the g µν ∂ µ φ∂ ν φ/2 term are proportional tȱ φ, whereas those coming from √ −gV are ∼ ∂φV . This means that all couplings between metric perturbations and inflaton perturbations are slow-roll suppressed and can be ignored during inflation. The effective mass term due to the inflaton potential also vanishes in the slow-roll limit, Thus, the second order action for the Fourier components of the inflaton perturbations, without the slow-roll suppressed terms, reduces to (4.40) Note that due to the b sr pre-factor the kinetic term in the action is not canonicallynormalized. If one wishes, b sr can be absorbed into a field re-definition which makes the kinetic term canonical.
Having derived the second order action for the inflaton fluctuations, we can now quantize them. First we need the conjugate momentum density where we have taken the functional derivative of the Lagrangian and made use of the reality of the inflaton fluctuations, i.e., δφ * k = δφ −k . In the Heisenberg picture, the field operators δφ k and their conjugate momenta operatorsπ sr k satisfy the equal time commutators: (4.42) The only non-vanishing commutator and the expression for the conjugate momentum in eq. (4.41) yield (4.43) The quantized perturbations δφ k can be written in terms of operatorsâ k and mode functions The two mode functions, u k (τ ) and its complex conjugate, span the space of solutions to the classical equation of motion for δφ k obtained by varying the action in eq. (4.40), i.e., Given a set of initial conditions for the mode functions, we can evolve them forwards in time. To calculate δφ and ultimatelyR k we need to know not only the initial conditions for the mode functions, but also the commutators forâ k andâ † k . Since, during inflation the co-moving Hubble sphere shrinks, observationally-relevant co-moving modes lie inside the sphere at early times and cross out of it at some point before the end of inflation. Hence, at very early times, the physical wavelength of these modes is much shorter than the Hubble radius, and they are not affected by the space-time curvature. At these early times, the mode functions should resemble those of free fields in flat space-time, while theâ k and a † k operators should be interpreted as creation and annihilation operators. Note that the latter are time-independent, meaning their commutation relations remain the same even after Hubble exit. We can show all of this rigorously.
Consider eq. (4.45). Since we work in the slow-roll approximation, we can set H = ∂ τ a(τ )/a 2 = const, implying a(τ ) = −1/(Hτ ). We can then find the general solution exactly The two constant coefficients, c k1 and c k2 , can be found by imposing that at early times, when modes are still sub-Hubble, k/aH = −kτ 1, δφ behaves as a free, massless 12 scalar with creation and annihilation operators satisfying Using eq. (4.43), we find c k1 = −(2π) −3/2 H/ √ 2k and c k2 = 0. This means that mode functions of co-moving modes that are still sub-Hubble are of the form e −ikτ /a, i.e., plane waves with a decaying amplitude, scaling inversely with a(τ ). The scale factor is simply a consequence of the non-canonical kinetic term in the action in eq. (4.40). One can easily check that the rescaled field aδφ has a canonically-normalized action and the equation of motion for its mode functions is of the form given in eq. (4.45), with b sr = 1 and ω 2 sr = k 2 −∂ 2 τ a/a. The plane wave factor can also be checked to enforce the vacuum state that is annihilated byâ k , |0 , as the ground state of the second order Hamiltonian. This state is known as the Bunch-Davies vacuum. Conventionally, it is assumed that observationallyrelevant modes started in the Bunch-Davies vacuum while lying deep inside the Hubble sphere during inflation. As the universe expands these modes become super-Hubble and according to eq. (4.46), for k/aH = −kτ → 0, u k ≈ i(2πk) −3/2 H/ √ 2, i.e., the inflaton perturbations freeze. The canonically-normalized field, aδφ k , then simply grows linearly with scale factor, whereas its conjugate momentum is equal to aHδφ −k , and therefore the two effectively commute. Hence, the superhorizon inflaton perturbations behave classically, and can be treated as a classical stochastic field. The quantum expectation value with respect to the Bunch-Davies vacuum translates into the classical ensemble average over field realizations drawn from a Gaussian probability distribution. Since we work in the spatially-flat gauge, the curvature perturbation,R = Hδφ/φ, on super-Hubble scales is Gaussian, too. Its power-spectrum, ∆ 2 R , is defined as follows Note that the expressions on the right do not depend on the arbitrary space coordinate x. It can be understood in terms of the isotropy and homogeneity of the universe. On super-Hubble scales during inflation the power-spectrum is scale-invariant, i.e., independent of k, This result was derived in the de Sitter approximation in which H is approximated to be constant. The conservation of the co-moving curvature on super-Hubble scales can be proven to hold more generally, independently of the assumption of de Sitter expansion. We will discuss it further below. However, it is important to note that since H and H vary slowly with time during inflation, the conserved value for ∆ 2 R is slightly different for different k. Essentially, the value of the conserved power-spectrum is determined by the value of H and H at horizon exit, k = aH. This introduces a weak scale-dependence. The power-spectrum of the co-moving curvature perturbation generated during slow-roll inflation is This result is obtained in the spatially-flat gauge, ignoring interactions of the inflaton perturbations with those in the metric due to slow-roll suppression. The approximation breaks down towards the end of inflation, when the slow-roll coefficients become of order unity. However, this does not affect eq. (4.49) for observationally-relevant modes, since the expression is evaluated at the time of Hubble horizon exit, when the slow-roll approximation still holds. One can derive the above results by working in gauge-invariant variables only. Under the diffeomorphism given in eq. (4.29), the inflaton perturbation transforms as ∆δφ = −∂ τφ ξ 0 . We then define the gauge-invariant inflaton perturbation (4.50) The gauge-invariant co-moving curvature perturbation can be then expressed only in terms of gauge-invariant quantities, R = Ψ + Hδφ/φ; see eq. (4.35). The linearised equation of motion for the gauge-invariant field perturbation is

51) and the linearised Einstein equations yield
(4.52) These equations can be most easily derived in the Newtonian gauge, B = E = 0, in which the only non-zero scalar metric perturbations ϕ and ψ are equal to Φ and Ψ, respectively, whereas δφ = δφ and therefore, the gauge-invariant quantities should obey the same equations of motion as the perturbations in the Newtonian gauge. In the equation of motion for the gauge-invariant inflaton perturbation, eq. (4.51), the couplings to the Bardeen variables are slow-roll suppressed. Similarly, the source terms involving the inflaton perturbation in the Einstein equations, eq. (4.52), are also slow-roll suppressed (they also vanish in the limit k/aH 1). Thus, the evolution of the gauge-invariant δφ is identical to the one of δφ in the spatially-flat gauge. Furthermore, since the source terms for Ψ vanish during slow-roll inflation and the contribution of δφ to R dominates due to division by √ H , we find the same value for ∆ 2 R (k) as in eq. (4.49), but this time using gauge-invariant variables. We should also point out that Ψ (as well as Φ) plays the role of an auxiliary field. One can see this most easily by substituting for the scalar metric perturbations in eq. (4.51), using eq. (4.52) (4.53) This is a second-order ordinary differential equation for δφ. Its quantized solution is of the form given in eq. (4.44). This means that the scalar metric perturbations do not have their own creation and annihilation operators. They can be expressed in terms of δφ according to eq. (4.52) and do not represent gravitational radiation. It is also obvious that during slow-roll, only the first three terms in eq. (4.53) are important, as expected, so the same considerations as before apply to the initial conditions for the mode functions and ultimately the expression in eq. (4.49) can be shown to hold.
So far we have shown that R is conserved on super-Hubble scales during single-field slow-roll inflation. Using eqs. (4.51) and (4.52), we can obtain the equation of motion for the co-moving curvature perturbation This equation, often rearranged in a different form, is referred to as the Mukhanov-Sasaki equation. In the limit of k/aH 1 it has a constant solution and a decaying solution going like dτ /(a 2 H ). The constant solution is the relevant one for observations. It remains constant even after the end of slow-roll of inflation.
In fact, a theorem due to Weinberg [52] states that no matter what the constituents of the universe are, for scalar and tensor perturbations about an FRW background, in the limit k/aH 1 there always exist two adiabatic solutions, one constant and one decaying. Adiabatic solutions have the same ratio δs/ṡ for any 4-scalar, s. 13 In single-field inflation, there is only one degree of freedom, δφ (the scalar metric perturbations are auxiliary fields), which implies that there are two solutions to its second-order differential equation. Since there are only two solutions, they must approach the adiabatic limit for k/aH 1 according to Weinberg. One can check this by considering the gauge-invariant quantity known as the non-adiabatic pressure, δp nad = δp−δρṗ/ρ. In single-field models, it can be shown to vanish on super-Hubble scales. However, in multi-field models, the non-adiabatic pressure does not vanish necessarily. When it does not, the equation of motion for R has an additional source term, due to the non-adiabatic (entropy) perturbations. There are more than two solutions for R, implying that R is not generally conserved in these cases.
The tensor metric perturbations given in eq. (4.26),h ij , are gauge invariant and evolve independently of the matter instabilities. Theh ij represent gravitational waves. There are no constraint equations on them and they represent the gravitational degrees of freedom. Since the 3-tensorh ij is traceless and transverse, it has two degrees of freedom only. They are frequently denoted as h + =h 11 2 and referred to as the + and × polarizations, respectively. In this notation h 3i = h i3 = 0 (the transverse plane waves are propagating along the z-direction). Perturbing the Einstein-Hilbert term in the action in eq. (4.13), one can show that the second order action governing each polarization state is of the form given in eq. (4.40) with b = m 2 Pl a 2 and ω 2 = k 2 . Note that in deriving the gravitational waves action we do not make the slow-roll assumption. We can then follow the standard quantization procedure, eqs. (4.41-4.45), separately for each polarization state. The equation of motion governing the mode functions reduces to manifesting the free nature of the tensor perturbations. The calculation of the mode function evolution during inflation is then identical to the one in the spatially-flat gauge for the scalar perturbations after ignoring slow-roll suppressed terms. At early times, for modes lying deep inside the Hubble sphere, one can show that u . Later on, as k/aH = −kτ 1, the mode function freezes to a constant u . Like in the scalar perturbations case, one can again show that on superhorizon scales the canonically-normalized tensor perturbation operatorsĥ (+,×) c = am Plĥ (+,×) effectively commute with their conjugate momenta. Hence, the gravitational waves become classical and Gaussian. Their total power is conventionally given by where in the last line we define the tensor power-spectrum. Similarly to the power-spectrum of the co-moving curvature perturbation, see eq. (4.49), the tensor power-spectrum generated during slow-roll inflation is The tensor perturbations are generally conserved in the limit k/(aH) 1, just like the co-moving curvature perturbation. One can see that most easily from eq. (4.55), which shows that the mode functions become overdamped in the super-Hubble limit. Hence, again there is a constant and a decaying solution.
The weak scale-dependences in ∆ 2 R (k) and ∆ 2 t (k) are characterised by their logarithmic derivatives In other words, one can approximate the scale-dependences by simple power-laws The quantities n s and n t are known as the scalar and tensor spectral indices, respectively, k is the pivot scale, and A s and A t are the amplitudes of the scalar and tensor power-spectra, respectively. Normally, the tensor amplitude is normalized by the scalar amplitude The quantity is known as the tensor-to-scalar ratio. Slow-roll inflation predicts small values for the logarithmic derivatives in eq. (4.58), n s − 1 = −2 H − η H and n t = −2 H . All slow-roll parameters are evaluated at Hubble exit of the pivot scale, k = aH, during inflation. In deriving these expressions, we have made use of the identity d/d ln k| k=aH ≈ H −1 d/dt| k=aH , which follows from the assumption that during slow-roll inflation d ln a/d ln k| k=aH ≈ 1 as H ≈ const. The scalar and tensor amplitudes can be also written in terms of the Hubble slow-roll parameters during inflation, Pl H ), A t = 2H 2 /(π 2 m 2 Pl ) and r = 16 H with again all quantities evaluated at k = aH. Note that r = −8n t and is known as the consistency relation for slow-roll inflation. To connect with the shape of the inflaton potential in models of single-field slow-roll inflation, we recall that with all potential and potential derivative terms evaluated atφ =φ , corresponding to the inflaton value at the Hubble exit of the pivot scale. CMB observations [2] yield A s = 2.2 × 10 −9 , n s = 0.968 ± 0.006 and the constraint r < 0.11 at 95 % confidence level. They are consistent with adiabatic primordial fluctuations, as predicted by single-field inflation.
In the above analysis of cosmological perturbations, we made several approximations. We ignored the contribution to δT µ ν due to anisotropic stresses, Π µ ν . They are a complimentary source of perturbations to the isotropic pressure term, i.e., The anisotropic stress in single-field inflation is zero at linear order. Multi-field models involving scalar fields only, also have Π i j = 0 at the linear level. More complicated models with vector fields in some homogeneous and isotropic background configuration for instance, can feature a non-negligible Π i j . Even in the presence of anisotropic stresses, according to the Weinberg theorem, there always exist a constant and a decaying solution for the scalar and tensor perturbations on super-Hubble scales. We have also not talked about the vector metric perturbations. The reason is that according to Einstein equations, the vector metric perturbations are always redshifted away in the absence of sources.
The aim of this section was to show that inflation can make the universe homogeneous and isotropic at the level required by observations. However, this comes at a price. At the end of inflation, the universe is in a cold and non-thermal state. On the other hand the successful theory of big-bang nucleosynthesis calls for a universe very close to thermal equilibrium at temperatures at least around 1 MeV. That is why reheating is an integral part of inflationary cosmology. Any successful theory of inflation must give an account of the production of Standard Model matter out of the energy stored overwhelmingly in the oscillating inflaton condensate at the end of the period of accelerated expansion. Reheating should also include baryogenesis and perhaps the production of dark matter. In the remaining sections of the lecture notes we review our current understanding of reheating. The early transfer of energy, from the inflaton condensate to the fields it is coupled to, is the subject of the next section. The main focus is on preheating -the exponential particle production due to non-perturbative resonances and tachyonic instabilities. Section 6 discusses the non-linear dynamics ensuing after the fragmentation of the inflaton condensate, and the approach to thermalization. Sections 7 and 8 connect phenomenological models of reheating with High-Energy Physics models and cosmological observations. We should point out that the details of the reheating process depend on the underlying particle physics theory beyond the Standard Model. Since there are so many possible extensions of the Standard Model, it makes more sense to begin by studying simple toy models inspired by High-Energy Physics to clarify the relative importance of different reheating mechanisms. Many toy models of reheating allow for a thermal universe at the epoch of big-bang nucleosynthesis. To some this is disappointing, since it shows that the current precision of observations does not let us distinguish between different models of inflation and reheating, but to others it is encouraging, since it advocates the inflationary scenario.
-25 -5 Preheating: the decay of the inflaton condensate 'The career of a young theoretical physicist consists of treating the harmonic oscillator in ever-increasing levels of abstraction.'

Sidney Coleman
Around the end of inflation, H = 1, the homogeneous inflaton begins to oscillate about the minimum of its potential. The inflaton condensate must decay into other forms of matter and radiation, eventually giving the particle content of the Standard Model and perhaps dark matter. These more familiar forms of matter and radiation must eventually reach thermal equilibrium at temperatures greater than 1 MeV in order to recover the successful big-bang nucleosynthesis scenario. The transition of the universe from the supercooled state at the end of inflation to the hot, thermal, radiation dominated state required for big-bang nucleosynthesis is called reheating. The subject of this section is the early transfer of energy from the inflaton condensate to the fields it is coupled to. We begin with the perturbative theory of reheating -historically, the process was first treated this way. We then show the importance of non-perturbative effects arising from the coherent nature of the inflaton condensate. They include parametric resonances and tachyonic instabilities, all of which lead to exponential growth in the occupation numbers of the fields the inflaton decays to (i.e., the decay products). These kinds of rapid decay are called preheating, with the decay products in a highly non-thermal state. Finally, we discuss the implications from coupling these decay products to additional matter fields for the energy transfer from the inflaton condensate.

Perturbative treatment of reheating
Originally, reheating was studied as a perturbative process [61] in which individual inflaton particles were assumed to decay independently of each other. Interaction rates and decay rates were calculated in the usual manner, using perturbative coupling expansions. For illustrative purposes we consider inflaton interactions of the form S matter ⊃ d 4 x √ −g(−σφχ 2 − hφψψ), where χ and ψ are some scalar and fermion decay products. These sort of couplings arise in gauge theories with spontaneously broken symmetries. We avoid tachyonic instabilities in χ by assuming that its mass, m χ , is greater than σ|φ|. The inflaton potential is assumed to be V = m 2 φ 2 /2. To tree-level order, for decay products much lighter than the inflaton quanta, the decay rates are [62] The total width, Γ tot ≡ Γ φ→χχ + Γ φ→ψψ , is supposed to determine the decay rate of the number of inflaton quanta in a fixed co-moving volume where nφ = ρφ/m is the number density of inflaton particles in the condensate. Hence, a 3 (t)nφ(t) ∼ exp(−Γ tot t). Since after inflation m H ∼ t −1 , the homogeneous inflaton undergoes many oscillations during one Hubble time. If Γ −1 tot m −1 , then we can approximateφ ≈Φ(t) cos(mt), whereΦ(t) varies much more slowly than the phase. Using eq. (4.14), we then find that nφ ≈ mΦ 2 /2. Thus,Φ(t) ∼ a −3/2 (t) exp(−Γ tot t/2), which agrees with eq. (4.21) to leading order up to an extra exponential factor. We can check that this additional exponential decrease due to particle production can be roughly taken into account by including a friction term into the background equation of motion Having m H ∼ Γ tot , one can write the WKB ansatzφ ≈Φ(t) cos(mt) assuming the phase varies much faster than the amplitude. NeglectingΦ and HΦ terms, we then find that 2Φ + (3H + Γ tot )Φ = 0 as required. Even if m H Γ tot one can still show that the second order WKB solution is a −3/2 (t) exp(−Γ tot t/2) cos(mt).
For small coupling constants, as required for radiative corrections to not spoil the flatness of the potential during inflation, typically Γ tot H. At the beginning of the oscillatory phase, the inflaton condensate mainly loses energy due to the expansion of space. Once the Hubble rate has decreased to H Γ tot , the particle production becomes effective. Thus, the energy density transferred into decay products is ∼ 3m 2 Pl Γ 2 tot . Note that H Γ tot is one of the conditions for establishing thermal equilibrium between the inflaton particles and (at least one of) the decay products. Setting the decay rates into individual species to be comparable to each other, i.e., Γ φ→χχ ∼ Γ φ→ψψ , all decay products can be in thermal equilibrium provided they have sufficiently high concentrations. Thereby, we can find an upper bound on the reheating temperature. It is safe to assume that most of the energy has been transferred into the light (with respect to m) decay products. Assuming they are relativistic as well, the energy density of the universe is π 2 g * T 4 /30, where g * is the number of relativistic degrees of freedom, of order 10 2 for the Standard Model. The maximal reheating temperature is Recalling eq. (4.49) and the CMB bound on the tensor-to-scalar ratio we find that T reh < 10 15 GeV, implying that the GUT symmetries cannot be restored after inflation and the solution to the monopole problem is not in danger. However, this does not rule out the production of other dangerous massive relics such as gravitinos. They could ruin the predictions of the successful big-bang nucleosynthesis by leading to an unwanted matter dominated state of expansion at the beginning of the epoch or by releasing excessive amounts of entropy close to it. One needs to make sure that in this sort of models, the reheating temperature is low enough to avoid the thermal production of such relics. We should point out that since each Γ is proportional to the square of the small coupling constants, the perturbative decay is actually quite slow and can take many e-folds of expansion after inflation before the Hubble rate becomes small enough for perturbative particle production to become efficient.

Limitations
There are many issues with the above perturbative analysis. The heuristic equation of motion in eq. (5.3), while capturing the qualitative behaviour, does not provide a consistent description of even the perturbative decay of the condensate. It violates the fluctuation dissipation theorem which states that dissipation inevitably leads to fluctuations within the system at hand. The effects of these fluctuations on the effective mass of the inflaton condensate are not included in eq. (5.3) [14].
Another problem with the above perturbative approximation is that it does not account for the Bose condensation effects. Even if the couplings of the inflaton to bosons, e.g., to χ, are small enough to allow for a perturbative coupling expansion, if the phase space of bosonic decay products, e.g., of χ particles, is densely populated Bose condensation effects can greatly enhance the decay rate. We discuss this situation in Section 5.1.2.
Most importantly, for larger couplings (but still small enough for radiative corrections to be negligible) the perturbative methods fail. Particle production has to be treated as a non-perturbative effect. The inflaton condensate is a coherent oscillating homogeneous field, implying that particle production has to be treated as a collective process in which many inflaton particles decay simultaneously, not independently of each other. Due to the large occupation number, we can treat the condensate classically. However, the decay products have to be described quantum mechanically, since they have vanishing occupation numbers at the end of inflation (due to the enormous dilution of space during the accelerated expansion). It is justified to use their vacuum state as an initial condition for the ensuing quantum mechanical particle production in the classical inflaton background. The periodic time-dependence of the effective masses of the decay products in the classical oscillating background can have a powerful effect on their production rates in the form of a parametric resonance, which will be the subject of Section 5.2.
Despite all of these problems, the perturbative analysis in this section can be applied to the late stages of reheating, e.g., to the decay of remnant inflaton particles after most of the energy has been transferred into relativistic species. Note that such decay channels are crucial to include, so that the energy transfer can be completed. Otherwise, we can face another relic problem.

Bose condensation of decay products in the perturbative limit
We finish this section with a short discussion of the Bose condensation effects in the perturbative limit, σ m 2 /Φ. By a perturbative limit, we mean that the tree-level order Feynman diagram gives the dominant contribution to the decay of the condensate into χ particles. Higher-order Feynman diagrams are subdominant. They can describe the simultaneous decay of more than one inflaton particles from the condensate and are negligible in the perturbative limit to leading order. To avoid significant radiative corrections to the Lagrangian we also put σ m. Taking into account that the condensate is comprised of particles at rest with large occupation number n φ 0 , the decay rate to a pair of χ particles at tree-level order is proportional to whereas the rate of the inverse process is proportional to Note that the occupation number is the number of occupied states per (2π) 3 phase space volume. The only exception is the inflaton condensate for which n φ 0 /V com = nφ = mΦ 2 /2, whereas for the χ particles the number density, n χ , is related to the occupation number in the standard way n χ = d 3 k n χ k /(2π) 3 . Note that n χ k = n χ −k and are independent of the direction of k. From now on we put them to be equal to n χ k . Roughly speaking, due to energy and momentum conservation, a stationary inflaton particle decays into a pair of χ particles, each of which has energy m/2 and momentum [m 2 /4−m 2 χ −2σφ(t)] 1/2 . Since m m χ > σΦ, all particles are produced within a thin spherical momentum shell in phase space, centred near m/2 and of width 4σΦ/m. Hence, n χ k=m Then the rate of change of χ particles within a given co-moving volume is where n φ 0 {n χ k , 1} and |k| = m/2. For n χ k > 1, i.e, n χ > σnφ/(π 2Φ ), the second term inside the brackets in the last line in eq. (5.7) becomes important, which is a manifestation of Bose condensation effects becoming relevant. Since ρ χ /n χ ∼ m = ρφ/nφ, Bose effects should be considered for fractions of energy stored in the decay product satisfying ρ χ /ρφ > σ/Φ. For small coupling constants and large amplitudes, the right-hand side can be much less than unity and the equality can be satisfied shortly after inflation. Bose effects become important and the perturbative treatment presented in the beginning of this section breaks down. For high occupancies, n χ k 1, after ignoring the expansion of space, and using eq. (5.1) we can integrate eq. (5.7) to get Bose effects lead to an exponential increase of the decay efficiency. We have shown it for small enough couplings which allow for a perturbative treatment. When couplings are increased, non-perturbative effects become important, but the exponential increase in the decay efficiency remains. This is shown in the next two sections. A discussion of the effects on the particle production rate due to the expansion of the universe is also included.

Parametric resonance
As shown at the end of the previous section, Bose condensation effects can exponentially enhance the rate at which energy is transferred from the oscillating inflaton condensate to the bosonic fields it is coupled to. We worked in the perturbative limit in which the coupling is small, e.g., σΦ m 2 in the trilinear interaction model V (φ, χ) = m 2 φ 2 /2+m 2 χ χ 2 /2+σφχ 2 , with m m χ > √ σΦ. In this limit a perturbative coupling expansion makes sense. If the amplitude of inflaton oscillations and/or the coupling constant become large, e.g., σΦ > m 2 in the trilinear case, high-order Feynman diagrams give comparable predictions to the lowest-order ones and the problem has to be approached non-perturbatively. Note that simultaneous decays of more than one inflaton particles from the condensate are described by high-order diagrams. Such decays are a consequence of the coherent nature of the oscillating condensate and the non-perturbative calculation presented below captures them, unlike the perturbative one in the previous section. It turns out that Bose effects still exponentially enhance the rate of energy transfer. It is more efficient than in the perturbative limit, due to contributions from the simultaneous decays of more than one inflaton particle. The phenomenon can be understood most easily in the language of parametric resonance. Of course, the method can be applied to the perturbative case as well.
At the end of inflation, matter fields can be treated as fluctuations on top of the oscillating homogeneous inflaton background. Typically, they start in the vacuum state, since inflation has diluted the corresponding particle densities to vanishing values. Ignoring the expansion of space for now, the linearised equations of motion take the form where the angular frequency is periodic, i.e., ω 2 (k, t) = ω(k, t + T ) 2 ; T is the period of oscillations of the condensate. In the trilinear model, ω 2 (k, t) = k 2 + m 2 χ + 2σΦ cos(mt) and T = 2π/m. Unlike the previous section, here we do not assume anything about the relative values of m, m χ and √ σΦ. The equation of the form given in eq. (5.9) with ω a periodic function of time is known as the Hill's equation [63,64]. In the triliniear case ω evolves harmonically and the equation can be reduced to the Mathieu equation form with A k , q and z dimensionless and determined by the form of ω. In the trilinear model, A k = 4(k 2 + m 2 χ )/m 2 , q = 4σΦ/m 2 and z = mt/2. 14
The Floquet theorem [63] states that the most general solution of the Hill's equation is given by where µ k is called the Floquet exponent and P k± (t) = P k± (t + T ). If (µ k ) = 0 one of the two terms increases exponentially with time. This is called parametric resonance. Let's prove eq. (5.13) [17,65] and show how to find numerically µ k [17,66], knowing the form of ω(k, t).
If u k (t) is a solution to eq. (5.12), then so must be u k (t + T ). Hence, if u k1 (t) and u k2 (t) are two linearly independent solutions, their time-shifted counterparts must be linear combinations of them, i.e., u ki (t + T ) i are the two eigenvalues of B ij and v ki (t) are independent linear combinations of u ki (t). From this follows that v ki (t+T ) = λ B i v ki (t), i.e., a time sift t → t+T leads to a rescaling by an eigenvalue. The most general solutions having this property are This completes the proof of eq. (5.13). The Floquet exponent is simply µ k = ln(λ B )/T , whereas each of the periodic functions P k± (t) is some linear combination of P k1,2 (t).
To find the Floquet exponent, we just need to calculate the eigenvalues of B ij , which, as we just showed, has a unit determinant. To do that we choose two orthogonal initial Hence, after evolving the Hill's equation forward for one period T for the two sets of initial conditions we can find the eigenvalues 15 (5.14) The initial conditions are relevant for the efficiency of the parametric resonance. Essentially, if both the initial field and field velocities are zero, parametric resonance does not lead to any growth. We can see this most easily by re-writing eq. (5.13) as a linear combination of the linearly independent v k1,2 (t) = exp(±µ k t)P k1,2 (t), i.e., u k (t) = c 1 v k1 (t) + c 2 v k2 (t). If both u k (t 0 ) andu k (t 0 ) are zero, then the only possibility for the constant pre-factors is c 1 = c 2 = 0. Hence, unlike ordinary resonance where the forcing term leads to a rapid growth even if initially the field displacement and velocity are zero, parametric resonance does not allow for any resonant excitations if no energy is stored in the fluctuations initially. That is why vacuum fluctuations, albeit small, play a crucial role for particle production after inflation as seeds for parametric resonance. 15 Using the fact that for our choice of initial conditions W [u k1 , u k2 ](t0) = 1 and that W [u k1 , u k2 ](t0 +  The dark areas correspond to vanishing real part of the Floquet exponent and are regions of stability. Narrow parametric resonance occurs for A k = n 2 and |q| → 0 (n is an integer) and broad resonance takes place when A k < 2|q| as well as A k − 2|q| |q| 1/2 . Resonant decay of the inflaton condensate into χ particles in the trilinear model, V (φ, χ) = m 2 φ 2 /2 + m 2 χ χ 2 /2 + σφχ 2 , can be understood in terms of the Mathieu equation -the equation of motion for χ k (after ignoring the expansion of space) can be mapped onto the Mathieu equation with A k ≥ 0 and q ≥ 0, i.e., the square region with northeast white lines (see also Fig. 3). Resonant production of χ particles in another common toy model, V (φ, χ) = m 2 φ 2 /2 + m 2 χ χ 2 /2 + g 2 φ 2 χ 2 /2, can be also mapped onto this chart -for g 2 > 0, the wedge-shaped region with vertical white lines, whereas for g 2 < 0 the region with northwest white lines, cover the relevant ranges of A k and q (see also Figs. 4 and 5).

Narrow resonance
As an exercise, we can now calculate the dimensionless Floquet exponent,μ k , of the Mathieu equation eq. (5.10). The magnitude of the real part ofμ k is plotted in Fig. 2 as a function of the parameters A k and q. We call plots of this type instability charts. There is a series of regions of stability in which (μ k ) = 0. They are surrounded by 'unstable' regions in which (μ k ) > 0. For |q| 1 and A k > 0, the regions of instability become narrow and approach A (n) k = n 2 as q → 0 (n is an integer). In the first narrow band the peak value of the Floquet exponent is (μ k ) (1) max ≈ |q|/2, while A k ≈ 1 ± |q| [63]. For the triliniear model, this corresponds to resonant production of χ particles with momentum in the range m/2 ± σΦ/m (assuming m χ = 0) and mode functions growing as exp(σΦt/m) (see also Fig.  3). Since the χ particles are described with the action of a time-dependent simple harmonic oscillator, eq. (5.11), the energy stored in a given k mode is simply where n χ k can be interpreted as the mean occupation number (mean, because it is evaluated by taking the expectation value of the Hamiltonian with respect to the Bunch-Davies vacuum). Hence, for q = 4σΦ/m 2 1, modes lying near the peak in the first narrow instability band have occupation numbers growing as n χ |k|≈m/2 ∝ exp(2σΦt/m). This is in Figure 3. The instability chart featuring the real part of the Floquet exponent normalized by the inflaton mass (left) and the Hubble rate (right), characterizing the χ particle production rate in the trilinear model, V (φ, χ) = m 2 φ 2 /2 + m 2 χ χ 2 /2 + σφχ 2 . The equation of motion for χ k can be reduced to the Mathieu equation, eq. (5.10), with A k = 4(k 2 + m 2 χ ) 1/2 /m, q = 4σΦ/m 2 , wherē Φ is the amplitude of inflaton oscillations (see also Fig. 2). In FRW space-timeΦ ∝ a −3/2 and k ∝ a −1 , implying that a given co-moving mode flows towards the bottom left corner of the chart as the universe expands as indicated with the white lines in the second chart (drawn for m χ = 0 for simplicity). Note that resonance is efficient if (µ k )/H ∼ σm Pl /m 2 1.
good agreement with the perturbative treatment of Bose condensation from the previous section, see, e.g., eq. (5.8). Thus, in the perturbative limit, σΦ/m 2 1, the Bose effects, due to the population of χ modes, in the leading order φ → χχ Feynman diagram can be described as a parametric resonance due to the first, n = 1, narrow, q 1, instability band. Higher order, n > 1, narrow bands lead to production of particles with momentum in the range k (n) = nm/2 ≥ m. They correspond to higher order Feynman diagrams describing the simultaneous decay of n φ particles from the condensate into a pair of χs, taking into account the Bose effects due to the dense populations of the χ modes, n χ |k|=nm/2 > 1. Since this happens in the perturbative limit, one should be able to describe it using the methods from the previous section, leading to eq. (5.8). In summary, resonance from the narrow bands, |q| 1, describes perturbative decays of particles from the inflaton condensate in the trilinear model, taking into account the occupation of χ modes. This type of parametric resonance is known as narrow resonance. It corresponds to |q| 1 for the Mathieu equation, but for the general Hill's equation it corresponds to the parametric resonance in some region of parameter space which features a narrow instability band.

Broad resonance
Similarly, the term broad resonance is used to describe parametric resonance in broad instability bands in parameter space. For instance, it occurs if |q| 1 for the Mathieu equation, see Fig. 2. This corresponds to the non-perturbative limit in the trilinear model. In this limit, the only means for calculating the particle production is by solving the mode equation, eq. (5.12), and a very intuitive way of describing its solutions is the Floquet Figure 4. Same as Fig. 3, but for V (φ, χ) = m 2 φ 2 /2 + m 2 χ χ 2 /2 + g 2 φ 2 χ 2 /2, for which A k = (k 2 + m 2 χ )/m 2 + 2q, 2q = g 2Φ2 /(2m 2 ). The charts are for g 2 > 0, see Fig. 2. Broad resonance occurs only for A k − 2q q 1/2 , i.e., for low-momentum modes (k 2 + m 2 χ )/m 2 gΦ/m, for specific ranges of gΦ/m. analysis we have developed. Broad resonance is much more efficient than narrow resonance since a broad, continuous range of k modes is excited. The typical rate of excitation is comparable to the background oscillation rate, | (µ k )| ∼ T −1 , and is much greater than in narrow resonance. The reason why the period of inflaton oscillations is the characteristic time-scale for particle production can be understood from the fact that in broad resonance, particles are produced in bursts, rather than smoothly as in the narrow resonance. Those bursts are separated in time by ∼ T . They occur every time the adiabadicity conditioṅ is violated. Since broad resonance occurs in the non-perturbative regime, where interactions with the inflaton background determine ω(k, t), and since their magnitude varies with period T , the adiabadicity condition is violated each time the background value of the inflaton is such that the interaction terms vanish -thenω(k, t) ω 2 (k, t). For an oscillating field, this happens twice a period, implying a rate of particle production comparable to T . Note that in the narrow resonance the adiabadicity condition is always satisfied, since interactions are weak (they can be treated perturbatively) and ω 2 (k, t) ≈ k 2 = const always. The only reason for resonance is the dense occupation of χ modes, which leads to a smooth exponential increase in the occupation numbers of particular modes. The reason why the case of broad resonance is different can be understood qualitatively by considering the mode functions in the adiabatic and non-adiabatic regimes. In the adiabatic limit, the WKB solutions to eq. (5.12) are The vacuum state mode functions which minimize the Hamiltonian corresponding to the action in eq. (5.11) and which are such that the commutators in eqs. (4.42) and (4.47) are Figure 5. Same as Fig. 4, but for g 2 < 0, (see also Fig. 2). Note the prominent broad resonance bands, corresponding to A k 2q < 0, not present in the instability chart given in Fig. 4 for g 2 > 0.
satisfied forχ k andπ χ k , and a χ k and a χ † k , respectively, have |α k | = 1 and β k = 0 [67]. More generally, mode functions satisfying the field-momentum commutator in eq. (4.42), given eq. (4.47), correspond to |α k | 2 − |β k | 2 = 1. One can show that this expression is consistent with the constancy of the Wronskian, since 3 and is equal to i/(2π) 3 if we start in the vacuum. The mean occupation number, see eq. (5.15), is simply n χ k = |β k | 2 , i.e., an adiabatic invariant and equal to zero in the vacuum state. We should note that for |β k | > 0, the Bunch-Davies vacuum is no longer an eigenstate of the Hamiltonian.
For instance, the adiabaticity condition, eq. (5.16), for the Mathieu equation reduces to 2q sin(2z) (A k + 2q cos(2z)) 3/2 1 , implying that if A k 2|q|, the inequality is not satisfied near z j = π/4, 3π/4, ... and the WKB solution, see eq. (5.17), does not hold. Away from these z j , the WKB solution is a good approximation. If A k 2q > 0, adiabadicity can be also violated for A k − 2q q 1/2 near z j = π/2, 3π/2, ... (similar expressions hold for A k −2q > 0; this provides a qualitative explanation of the broad bands in Figs. 2, 3, 4, 5). In general, since |α j k | 2 − |β j k | 2 = |α j+1 k | 2 − |β j+1 k | 2 = 1, where the superscript j labels the coefficients between the jth and (j + 1)th violation of adiabadicity, etc., the connection between these Bogolyubov type coefficients is where θ j k = t j t 0 ω(k, t)dt is the accumulated phase until the jth violation of adiabadicity, and the reflection and transmission coefficients must obey |R j k | 2 + |D j k | 2 = 1, to preserve the Bogolyubov nature of the αs and the βs. The calculation of the reflection and transmission coefficients is tedious -one has to derive connection formulae for the WKB solution on both sides of the non-adiabatic region [68] for a given form of ω(k, t). Nevertheless, one can use the general expression in eq. (5.19) to study particle production. If we assume that we start in the vacuum state, i.e., β j=0 k = 0, there will be particle production after the first violation of adiabadicityn χ,j=1 20) where ∆θ j k = arg(R j k α j k β j * k ). In the limit n χ,j k 1, we can write n χ,j+1 The argument of the logarithm can take values greater or smaller than unity. µ j k > 0 corresponds to non-adiabatic particle production at event j. We note that violation of adiabaticity is a necessary, but not sufficient condition for particle production. The actual form of ω(k, t) must be such that R j k and θ j k + ∆θ j k allow for µ j k > 0, at least on average. That is why regions in parameter space in which the adiabatic condition is not satisfied can still contain stability bands, e.g., see the narrow regions of stability for A k 2q in Fig. 2 for the Mathieu equation.

Classical limit
The last point we wish to make for parametric resonance after inflation and preheating in general, is that the exponentially amplified modes can be treated classically. Intuitively, this can be understood from the large occupation numbers of these modes. Quantitatively, we can see it by considering the field-momentum commutator The expectation values of the operator products on the left and right hand sides of the expression grow as e 2| (µ k )|t and can become much greater than unity. Their difference, however, remains small and constant. It is equal to the delta-function term. We can check this by evaluating the commutator explicitly Since the Wronskian W [u k (t), u * k (t)] = const for any equation of the form given in eq. (5.12), and since we start with vacuum fluctuations (β k = 0) W [u k (t), u * k (t)] = i/(2π) 3 always, even if each of the terms has grown exponentially. This proves that even though the commutation relation is respected, the quantum correction delta-function term affects the expectation value ofχ k (t)π q (t) negligibly for the resonantly amplified k. Hence, the quantum expectation value with respect to the Bunch-Davies vacuum of any function of the densely populatedχ k mode can be treated as a classical ensemble average over field realizations drawn from a Gaussian probability distribution. The variance of a (zero-mean) field in the WKB regime at some time t is where (2π) 3 |u k (t)| 2 = [1+2n χ k +2 n χ k (n χ k + 1) cos(γ k )]/(2ω(k, t)) ≈ 2n χ k cos 2 (γ k /2)/ω(k, t) for n χ k 1, where γ k = arg(α k β * k ). Similar considerations apply to more complicated functions which depend on time derivatives ofχ k (t) as well.
In the following section we discuss how gravity and additional oscillating background fields can affect the resonant particle production described here.

Stochastic resonance
In the previous section we showed that parametric resonance can play an important role in the preheating phase. We considered the growth of matter fields, represented as fluctuations in an oscillating background, by applying Floquet analysis to the linear equations of motion with periodic coefficients. In doing so, we made several simplifying assumptions. In this section we re-introduce some of the ignored effects and show that they lead to a phenomenon known as stochastic resonance.
Neglecting gravity and assuming that the inflaton is the only field that has a background value allows for the possibility of having strictly periodic linear equations of motion, with exponentially growing solutions. One expects that any extension beyond this set-up can spoil the exact periodicity and, in general, counteract the growth of perturbations.

Metric fluctuations
Actually, incorporating gravity is not difficult. The metric perturbations remain negligible while particle production takes place. One can see that from the generalized Poisson equation (which follows from a combination of the Einstein equations) where, Ψ, is the Bardeen potential, see eq. (4.35), and δρ m ≡ δρ + [ρ (τ )/a(τ )]δu is the comoving, gauge-invariant, density perturbation. 16 After defining the co-moving overdensity field δ m = δρ m /ρ, we can say that the linearized equations of motion governing the perturbations hold for small δ m 1. The Fourier transform of eq. (5.25) is Ψ k = (3/2)(aH/k) 2 δ mk , implying Ψ k → 0 as aH/k 1 for small δ mk . Hence, metric perturbations remain vanishingly small on sub-Hubble scales during the preheating phase. During this phase, the super-horizon metric perturbations also do not grow in single-field models of inflation according to Weinberg's adiabatic theorem [52].

Expansion of space
Unlike the metric perturbations, the background space-time curvature cannot be easily neglected during preheating. The FRW expansion of space causes the amplitude of inflaton oscillations to decay, while co-moving wave-numbers are red-shifted to smaller physical values. Going back to our parametric resonance approach, we can see that the equation of motion for the scalar matter fields, eq. (5.9), can still be reduced to the form of a simple harmonic oscillator with a time varying frequency. Using the canonically-normalized field χ c (t) = a(t) 3/2χ (t), where t is cosmic time, we obtain In the trilinear model, see Sections 5.1 and 5.2, ω 2 (k, t) = (k/a) 2 + m 2 χ + 2σ 2Φ (t) cos(mt) − (3H/2) 2 − (3/2)Ḣ, implying that this is not the Hill's equation any more. Nevertheless, one can depict qualitatively the effects from FRW expansion on particle production by adding flow lines to the Floquet chart, tracing the evolution of particular co-moving modes. Since, in m 2 φ 2 /2,Φ(t) ∼ a(t) −3/2 and 3H 2 ≈ −2Ḣ a given co-moving mode k flows exactly alonḡ Φ ∼ k 3/2 phys ≡ (k/a) 3/2 curve in the k phys −Φ plane, see right panel in Fig. 3 (see also Figs. 4, 5 for other models). Empirically, a condition for parametric resonance (both narrow and broad) to result in significant particle production is for sufficiently long times. This is another way of saying that particle production occurs only in those bands in which the resonant growth is rapid on the Hubble time-scale. Using the heuristic picture of Floquet theory, we can conclude that depending on the model, broad resonance can be enhanced or shut off by the expansion of space. When more and more co-moving modes are redshifted towards a broad instability band, we observe a temporary increase in the net particle production, see bottom left corner of right panel in Fig. 3, but as they eventually leave the instability regions the resonance gets completely shut-off. As we showed in the previous section, broad resonance can be described as a series of particle creation events in which the adiabaticity condition, eq. (5.16), is violated. Taking into account the effects of the expansion of space, implies that the quantities appearing in eq. (5.21) will be time-dependent. The reflection coefficient, R j k , should have some model dependent and usually monotonic time-dependence, whereas the phase, θ j k + ∆θ j k , can be assumed to vary randomly in the interval [0, 2π) . The fact that the Floquet index µ j k in eq. (5.21) can change stochastically between successive particle creation events is the reason why broad resonance in an expanding space is called stochastic resonance. On average implying an increasing number of particles, in agreement with entropic arguments. Note that due to the randomness of the phase, µ j k on average can be smaller than in the Minkowski space-time limit. This is a curious feature of stochastic resonance, where particle production occurs on time-scales much shorter than the Hubble time, but still the expansion of space affects the final result.
On the other hand, the efficiency of narrow resonance is severely degraded by the FRW expansion. As one can see in the Floquet charts in Figs. 3, 4, 5, co-moving modes cross the narrow instability bands much faster than in the broad resonance regime. Thus, expansion takes particles out of the thin resonance layers and the occupation numbers boosting the Bose condensation effect become smaller than in the Minkowski limit. If the rate of escape of particles is greater than the rate of their production, i.e., eq. (5.27) does not hold, then Bose effects play no role. The efficiency of narrow resonance is sensitive to other suppressing effects such as the re-scattering of the newly created particles out of the resonance layer, as well as the shift of the resonance region from its original location due to the change of the inflaton effective mass as a consequence of particle production.
We also note that after including the expansion of space we are still allowed to treat the heavily populated modes classically. In particular, the analysis after eq. (5.22) still holds for the canonically-normalized fieldχ c (t).

Multi-field preheating
The periodicity of the time-dependent background can be violated also if there are several oscillating homogeneous fields. Even without expansion of space, unless the motion at the background level in the multi-field space occurs along special trajectories such as Lissajous curves or effectively one-dimensional oscillatory trajectories, the time-dependent coefficients in the linear equations of motion governing the fluctuations are not exactly periodic. This can again lead to stochastic resonance if the adiabaticity condition, eq. (5.16), is violated [69]. Note that this time both the reflection coefficient, R j k , and the phase, θ j k + ∆θ j k , can be assumed to vary randomly between successive non-adiabatic events. Even the length of the time intervals separating such events can vary randomly. Nevertheless, just like in the case of an expanding space, we could approximate the motion in field space at the background level as being periodic to check if substantial instability (both broad and narrow) bands exist.
We should point out that there is an alternative description of resonant particle production when the number of oscillating homogeneous fields is much greater than one. In this case the effective masses of the daughter fields evolve with a random component to a very good approximation. This reduces the efficiency of the particle production, but resonance still takes place. It occurs at all wavenumbers, not only within particular resonance bands. The alternative way to see why this happens is to note that there is a duality between the equation of motion of daughter fields, see eq. (5.26), and the time-independent one-dimensional Schrodinger equation. The duality interchanges time and space, the modefunction with the wavefunction, the time-dependent effective mass squared with the spacedependent one-dimensional potential energy and k 2 with the eigenenergy. Then recalling the celebrated condensed matter phenomenon of Anderson localization, in which small random impurities make eigenfunctions exponentially localized in space, we expect that in the case of preheating, time-dependent masses with random components give rise to exponentially growing modes at all wavelengths; for more details on the condensed matter analogue and the random resonance see [69].
We have shown that realising a strictly periodic motion at the end of inflation is difficult. The FRW expansion and the possibility of having more than one oscillating homogeneous fields can lead to a quasi-periodic motion at the background level. This can lead to stochastic resonance if the adiabaticity condition, eq. (5.16), is not respected. Even if it is, there could be still some particle production due to perturbative decays. However, as opposed to the strictly periodic case, the Bose enhancement of decays into scalar fields is normally not significant. Despite all that, Floquet analysis remains an important first step towards understanding the instabilities in the evolution of matter fields during preheating.

Tachyonic decay
So far we have assumed that the effective frequency, ω 2 (k, t), of the matter fields, χ c , changes (quasi) periodically with time due to the inflaton oscillations. This need not be the case always. For instance, towards the end of Hybrid inflation [70], , as the inflaton becomes smaller than a critical value, φ 2 < λ χ v 2 /g 2 , but is not oscillating, the sign of ω 2 (k, t) changes from positive to negative for long-wavelength modes and can remain such for an extended period of time. Since one of the two imaginary frequency solutions to eq. (5.26) is exponentially growing with time, χ c ∝ e |ω|t , we again have exponential particle production. A negative squared frequency, ω 2 (k, t) = (k/a) 2 + m 2 χeff < 0, implies an imaginary effective mass, m 2 χeff < −(k/a) 2 < 0. That is why this mechanism for particle production is dubbed tachyonic preheating. Importantly, all modes whose momenta are less than the magnitude of the imaginary effective mass are unstable, and in the limit k → 0 the exponential index approaches the maximal value of |m χ,eff |. Tachyonic instabilities in fluctuations always occur in symmetry breaking models for small background field values, e.g., in Hybrid inflation for small enough φ. Tachyonic instabilities can be also observed in the fluctuations of the inflaton field itself, e.g., when it has a symmetry breaking self-interaction potential or in field ranges where the self-interaction potential is shallower than quadratic.
Just like in the case of resonant particle production, to have efficient tachyonic decay of the inflaton condensate, |m χ,eff | H 1 (5.28) must hold for a sufficiently long time. Otherwise, the expansion of space drives m 2 χ,eff to its equilibrium, positive value (implying positive ω 2 (k, t)) before substantial energy can be transferred from the condensate to fluctuations.
In general, tachyonic instabilities can be achieved in models with negative couplings. For instance, in the trilinear model in Sections 5.1 and 5.2, the interaction term σφχ 2 implies that even if the inflaton is oscillating, half of the period small k modes will be tachyonic. This corresponds to the 0 < A k < 2q region in the Mathieu instability chart in Fig. 2, and explains why there the stability bands are so narrow (see also Fig. 3). They correspond to the small parameter region in which effectively only the exponentially decaying imaginary frequency solution is excited. Note that the expansion of space blurs the boundaries between different regions in the Floquet chart and the narrow stability bands in the tachyonic region go away. Another example of negative coupling resonance is the models with a g 2 φ 2 χ 2 /2 interaction, where g 2 < 0, see Fig. 5. This implies q < 0 in the notation of the Mathieu equation, see Fig. 2. Note that in models like this, where interaction terms are always negative to ensure stability we should add higher order positive potential terms. In this case, we can add quartic potential terms, that dominate at large field values, but are unimportant during preheating. In terms of the Mathieu instability chart the tachyonic region corresponds to 2|q| ≥ A k ≥ 2q, where the latter bound comes from the q-dependence of A k in this model, see Figs. 2, 5. Compared with the standard resonant preheating scenario (g 2 > 0) where µ max k m, see Fig. 4, tachyonic preheating can be much more efficient, with maximal exponential index ∼ |g|Φ. Even if the couplings are small, |g| 1, to ensure negligible radiative corrections, we can still have |g|Φ m at the end of inflation. Typically,Φ ∼ m Pl , and even with small couplings it can take less than one oscillation of the condensate for the tachyonic growth of the long-wavelength modes to lead to interesting non-linear dynamics.

Instant preheating
The time-dependent nature of the effective mass of fluctuations can give rise to another preheating mechanism. Normally, for a coupling of χ to some fermion ψ of the Yukawa form, h χ χψψ, the decay χ →ψψ is kinematically forbidden if the corresponding bare masses are such that m χ < 2m ψ . However, if the scalar is coupled to the inflaton via g 2 φ 2 χ 2 /2 (assume g 2 > 0) then the effective mass, m 2 χeff = m 2 χ + g 2 φ 2 , can become significantly bigger. And even for a scalar of vanishing bare mass, the decay can be kinematically allowed. For an oscillating inflaton with large enough amplitude,Φ > 2m ψ /g, the decay rate, see also eq. (5.1), vanishes whenφ ≈ 0, and is maximal as the oscillating inflaton reaches its maximal value |φ| =Φ. In the large coupling limit, √ q = gΦ/m 1, we have broad resonance, or in other words non-adiabatic particle production every time the non-adiabaticity condition given in eq. (5.16) is violated. This happens whenφ ≈ 0, implying that Γ χ→ψψ is maximal half-way between two consecutive particle production events. Hence, even if a significant amount of χ particles are produced at each creation event, they can all decay into fermions before the next one. This mechanism is called instant preheating. In it, the back-reaction of χ particles on the φ condensate is slowed down and the efficiency of the resonance maintained for very long times. Furthermore, for g ∼ 10 −2 andΦ ∼ m Pl the light inflaton, m ∼ 10 −6 m Pl , can decay to heavier scalars and fermions, as heavy as the GUT scale ∼ 10 16 GeV. The return of the GUT scale into play obviously presents a threat to inflationary models. Far-fromequilibrium production of topological defects can take place, thus allowing cosmological observations to place bounds on different preheating scenarios.
6 Non-linear reheating 'Using a term like non-linear science is like referring to the bulk of zoology as the study of non-elephant animals.' Stanislaw Ulam As inflation ends, non-perturbative phenomena such as stochastic resonances and tachyonic preheating can amplify quantum fluctuations of the matter fields, creating particles in a far-from-equilibrium state. The instabilities grow exponentially fast on cosmological timescales. Such exponential growth cannot proceed forever. Eventually, the produced particles back-react on the preheating process. Mode-mode couplings and non-linear interactions become important. Soon the inflaton condensate fragments and non-linear dynamics takes over. The subsequent evolution of the bosonic fields can be rather non-trivial and a lot of interesting things can happen. Towards the end of this out-of-equilibrium evolution, the fields must thermalize, marking the end of reheating and setting the scene for big-bang nucleosynthesis.
This section begins with a discussion of the end of preheating. We talk about the various places back-reaction can arise in and terminate preheating. We then focus on the non-linear dynamics following the initial burst of particle production and the fragmentation of the inflaton condensate. We survey the different numerical techniques available for tackling the non-linear evolution, and also review the various non-trivial structures that have been studied. We finish with a discussion of the approach to thermalization which can include the turbulent evolution of scalar fields.

Back-reaction: the end of preheating
Resonant particle production and tachyonic instabilities can be terminated in various ways. If the expansion of space does not intercept the non-perturbative particle production, then the back-reaction of the produced particles eventually shuts it off. Back-reaction effects are associated with higher order in field fluctuations correction terms to the equations of motion in the approximate picture of preheating in which the inflaton condensate is treated as a time-dependent classical background with quantum field fluctuations on top of it.

Back-reaction at the background level
The equation of motion describing the evolution of the classical inflaton background, eq. (4.15), can have corrections due to non-vanishing spatial averages of interaction terms. For instance, in the V (φ, χ) = m 2 φ 2 /2 + g 2 φ 2 χ 2 /2 model, the presence of χ particles alters the effective squared mass of the inflaton condensate oscillations by ∆m 2 φ = g 2 χ 2 . The angle brackets represent a volume average of the classical χ (classical in the sense described at the end of Section 5.2). From now on, when discussing back-reaction and non-linear dynamics, we shall treat all bosonic fields classically and drop their hats. If there are exponentially unstable modes, then where µ is some effective growth index, close to the maximal one µ max k . 17 The coefficient of proportionality varies slowly with time (apart from an oscillating modulation, it decays monotonically due to the expansion of space). Hence, back-reaction effects become important, ∆m 2 φ ∼ m 2 , within, up to logarithmic factors, ∆t br ∼ µ −1 , (6.2) from the beginning of particle production. For broad resonance, µ ∼ m H, and so ∆t br is very short in comparison to the Hubble expansion time-scale. The effect on the condensate from the increase in its effective mass is that its amplitude of oscillations,Φ, decreases whereas its frequency increases. In the Mathieu equation notation, the resonance parameter q = g 2Φ2 /m 2 φ rapidly decreases and soon it is difficult for the resonant production of χ particles to continue further.

Re-scattering and non-linearity
The equations of motion describing the field fluctuations are also affected by the particle production. Working in the mean-field/Hartree approximation in which different modes and fields evolve independently (are uncorrelated in time), i.e., χ * k−q χ k time ≈ 0 if q = 0, δφ * k−q χ k time ≈ 0 for all q, etc., 18 there are only correction mass terms, ∆m 2 χ = g 2 δφ 2 and ∆m 2 δφ = g 2 χ 2 , for χ and the inflaton fluctuations, respectively. They may change the evolution of the field fluctuations slightly, e.g., shift χ particles out of resonance bands. However, as the number of particles increases, the mean-field/Hartree approximation stops being a good description. The coupling between different Fourier modes becomes important, heralding the true beginning of the non-linear stage. The mode-mode coupling between different momentum modes is called re-scattering and is what actually leads to the fragmentation of the inflaton condensate. For instance, there is an additional non-vanishing source term in the equation of motion for the inflaton fluctuations Having an inhomogeneous equation with exponentially growing source term, implies that its particular solution also grows exponentially, i.e., δφ q ∝ e 2µt . Thus, due to the interactions of pairs of χ particles with particles in the condensate, inflaton fluctuations grow twice as fast. The growth is a manifestation of inflaton particles being scattered out of the inflaton condensate. They are low-momentum excitations, predominantly. When δφ 2 Φ 2 , we say that the condensate is substantially fragmented and if δφ 2 Φ 2 we say that it is completely fragmented or destroyed. Rescattering also re-distributes the energy stored in the χ particles. Parametric resonance leads to the excitation of χ momentum modes lying in instability bands. The re-scattering transfers energy from the amplified modes to modes with momenta lying in the stability regions. This may slow down the resonant particle production. Even if it completely shuts off the resonance, re-scattering becoming important is a sign of the ensuing non-linear evolution and fragmentation of the inflaton condensate. 17 Note that according to the Ergodic theorem [52], the spatial average of χ 2 is also equal to the ensemble average over realizations of the stochastic field. This is what the vacuum expectation value of the quantum field tends to, since 0|χ(t) 2 |0 = k 3 |u k (t)| 2 d ln k/(2π 2 ) ∝ e 2µt , where we integrate the mode function. 18 The time average is taken over several oscillations of the more slowly oscillating Fourier transform.
We should point out that having a fragmented inflaton condensate, δφ 2 Φ 2 , does not necessarily imply that the energy stored in it is negligible. However, we can say with certainty that re-scattering and fragmentation kick in when the energy stored in interaction terms and/or fluctuations is comparable to the energy of the classical background. Hence, non-perturbative particle production ends and non-linear evolution begins with either most or at least a non-negligible fraction of the total energy being stored in field fluctuations.
Before moving forward to different approaches for studying the non-linear stage, we consider the possibility of having a second field that has a small, but non-vanishing background value, e.g., Φ |χ| > 0. While the back-reaction mechanisms remain largely unchanged, the preceding linear evolution during preheating can exhibit novel behaviour. Essentially, there are additional mixing terms, ∼ g 2 Φχδφ k and ∼ g 2 Φχχ k , in the equations for χ k and the inflaton fluctuations, respectively. They lead to chaotic evolution of the field fluctuations. The strong dependence on the initial value ofχ can give rise to observational signatures of preheating, as we will discuss in Section 8.3.
We now proceed with the non-linear stage of reheating, following the back-reaction of the produced particles on the inflaton condensate and the breakdown of the linear analysis.

Non-linear evolution
Preheating ends when the occupation numbers of excited bosonic field modes become large and back-reaction effects render the linearized approximation not applicable. The inflaton and the fields it is coupled to start evolving as a combined system. Non-linear interactions lead to the transfer of power between different wavenumbers. This non-linear phase is dynamically rich and can be studied numerically.

Numerical approach
The standard approach in numerical analysis is to solve the classical evolution equations, e.g., There are several publicly available codes created for this purpose. Most of them use a finite-difference method for solving the equations. The fields are discretized on a cubic co-moving spatial grid, with periodic boundary conditions. The time evolution is then a matter of evolving forward a system of coupled ordinary differential equations. For numerical integration in time LATTICEEASY [71] uses the simplest symplectic integratorthe leapfrog scheme, which is fast (no need for storage of field and field time derivatives simultaneously during a time step) and second order accurate in time. DEFROST [72] and HLATTICE [73] use higher order symplectic integrators. GABE [74] uses a second order Runge-Kutta method which stores field and field time derivatives on the same time slices, unlike symplectic integrators. Although, this requires more time to run the simulations and more physical memory, it allows for non-canonical kinetic terms. CUDAEASY [75] and PYCOOL [76] are GPU-accelerated codes based on DEFROST. A pseudo-spectral code, PSpectre [77], is also available, which evolves the Fourier transforms of the fields. In it, unlike in finite-difference codes, Laplace terms are dealt with straightforwardly. Each contributes a single term, e.g., k 2 φ k , to the Fourier transformed equations of motion with no computational cost, whereas ∆ x φ(x) in finite-difference codes is more costly, since one has to compute the differences with neighbouring points for each lattice site. However, non-linear interaction terms in the equations of motion, e.g., g 2 φ 2 (x)χ(x), are easy to deal with in finite-difference codes, whereas for pseudo-spectral codes they present a problem, since there one has to calculate multidimensional integrals.
We should also point out that most publicly available codes do not include metric perturbations, i.e., they evolve the fields in pure FRW space-time. In addition to the Klein-Gordan equations, they solve one equation for the evolution of the scale factor, a(t). Note that the Einstein equations yield two equations for the evolution of a(t), namely the Friedmann and Raychaudhuri equations given in eq. (4.4), with ρ(t) = ρ and p(t) = p averaged over the simulation box. Programs typically evolve the Raychaudhuri equation and treat the Friedmann equation as a constraint that has to be satisfied after each time step. Empirically, violations of the Friedmann equation ≥ 0.1 % indicate poor energy conservation and render the simulations unreliable. Some studies simplify matters further, by assuming a fixed time-dependence of a(t). This means that the expansion of space is not calculated self-consistently, e.g., by solving the Raychaudhuri. Common choices are a ∝ t n with n = 2/3, 1/2 for matter and radiation-dominated backgrounds, respectively. But still the Friedmann equation is treated as a constraint that has to be checked after each time step. Approximating the space-time to be FRW is justified, since the lattice size of typical simulations is sub-horizon and just like during preheating, metric perturbations are suppressed on these scales and do not affect the non-linear evolution of the fields. The reason why sub-horizon scales are of main interest are the causal mechanisms which drive the non-linear evolution of the fields. Causally disconnected Hubble patches evolve independently and almost identically, implying that it is sufficient to capture one Hubble volume in numerical simulations. Otherwise, the only publicly available code that can include metric perturbations is HLATTICE.
Of course, the FRW approximation is non-viable if large sub-horizon inhomogeneities, that can lead to the formation of primordial black holes, are present during preheating. However, such inhomogeneities rarely form due to matter field instabilities [15,66]. They either occur in models in which significant super-horizon inhomogeneities generated during inflation re-enter the horizon during preheating, or are induced by gravitational instabilities which become important long after the end of inflation.
We should also point out that all publicly available codes are written for scalar fields. The GABE code can be adapted for gauge field dynamics, but its ability to respect the gauge constraints has not been fully tested yet, especially with charged scalar fields. We should also note that DEFROST differs from the other finite-difference codes. In it, instead of directly discretizing the equations of motion, the Lagrangian is discretized and then the corresponding equations of motion are evolved numerically.

Non-linear dynamics
Non-linear effects can become important even in the simplest models of reheating in which the interactions of the inflaton with other fields are negligible, see Fig. 6. If self-interaction terms, e.g., ∝ φ n , n = 2, become important, the inflaton condensate can fragment after selfresonance. It can form non-trivial field configurations such as oscillons which can lead to long periods of matter-dominated state of expansion [36], or form Q-balls if the inflaton is a complex scalar [78]. Oscillons (as well as Q-balls) can also affect predictions in baryogenesis models with a complex inflaton [79]. If the inflaton is very light, but self-interacting, it inevitably fragments and attains a radiation-like equation of state [34]. Gravitational waves can also be generated due to fragmentation induced by self-interactions [44,80]. Even if the inflaton is not self-interacting, the condensate inevitably fragments due to gravitational instabilities [81].
Coupling the inflaton to other fields can lead to very rich phenomenology, see Fig.  6. Interaction with scalar fields can lead to the formation of metastable bubble-wall-like configurations [84], whose size is not much smaller than H −1 . After they collide the density is transferred to much smaller scales. This could be interpreted as upscattering of modes to higher momenta due to non-linear interactions. Scalar field theories can also feature the formation of solitons and defects, such as domain walls and metastable global cosmic strings [20]. Models that include gauge fields lead to new phenomena as well [22]. Non-equilibrium phase transitions can lead to the formation of stable topological defects, whereas models with conformal couplings feature large-scale magnetic fields and a faster approach to a radiation dominated state of expansion [31,37].
The reason why topological defects can be produced during preheating, in theories which allow them, is non-equilibrium restoration of broken symmetries [21,85,86]. Just like in phase transitions, a negative bare mass squared can receive significant positive contributions from large field variances leading to a positive effective mass squared and temporary restoration of a symmetry. This can happen either during back-reaction and re-scattering, or after the fields enter the full non-linear regime. Once the expansion of space dilutes the energy enough, the symmetry is re-broken and topological defects can be produced. For instance, in models allowing the formation of cosmic strings, one has to wait for the mean energy density to become less than the potential energy in the central unstable maximum, filling the ring at the bottom of the potential. This mechanism leads to the production of strings on both sub and super horizon scales [20]. The sub-horizon strings are transient -they evaporate due to the emission of classical radiation. In general, the density of topological defects arising in such non-equilibrium phase transitions is determined by the correlation length of the fields shortly before the symmetry is restored. For inhomogeneous configurations, the correlation length scale can be considerably smaller than the Hubble scale unlike in the Kibble mechanism. Thus, preheating provides a mechanism for copious production of dangerous topological defects, even at GUT scales, which can have implications for inflationary models.  [82] and global strings (right) [66]. Middle row: formation of bubbles (left) [72] and their decay products (right). Bottom row: domain walls (left) [83] and cosmic strings (right) [22].

Turbulent scaling
The early stage of the non-linear evolution, following back-reaction and re-scattering, in models with scalar fields is dynamically rich and chaotic. Various transient non-trivial field configurations can form, wiping out details on initial conditions from inflation and preheating. Unless long-lived objects form, e.g., oscillons, stable defects, black holes, etc., the state of the fields eventually enters a highly inhomogeneous phase which can persist for very long times, much longer than the preheating and transient phases. It is characterised by a slow, but steady transfer of energy to higher momenta. Essentially, straight after the transient stage the field occupation numbers in the infrared quickly saturate to a power-law, n(k) ∝ k −3/2 , with a UV cut-off not much greater than the typical wavenumber of excited particles during preheating [18,19]. The power-law is non-thermal -for a thermal distribution of relativistic weakly interacting bosons we expect n(k) ∝ k −1 in the infrared. It then slowly propagates towards higher momenta. Typically, the cascading of the distribution towards the UV can be characterized as turbulent scaling in which the occupation numbers evolve self-similarly, n(k, τ ) = (τ /τ 0 ) −q 1 n 0 (kτ −q 2 ), where q 1 and q 2 are some positive powers, determined by the form of the interactions and τ 0 is the conformal time when the scaling regime begins. This slow fragmentation proceeds until the occupation numbers of the highest k-modes belonging to the power-law distribution become of order unity. Then the classical description breaks down and quantum effects become important. Note that the energy density in a given mode is ρ(k) ∝ k 4 n(k), implying that the high-k modes belonging to the power-law dominate the energy budget. That is why, if present, their quantum behaviour cannot be neglected.

Two stages
None of the preheating mechanisms described in Section 5, nor the subsequent non-linear evolution yield a thermal spectrum of decay products. However, measurements of the anisotropies in the CMB and the relative abundances of light elements tell us that the Standard Model degrees of freedom were in thermal equilibrium at the beginning of the big-bang nucleosynthesis and that the universe at that time was in a radiation-dominated state [1]. The moment when the universe achieves thermal equilibrium for the first time after the end of inflation, at some reheating temperature, T reh ≥ T BBN ∼ 1 MeV, in a radiation-dominated state of expansion, w ≈ 1/3, marks the end of thermalization and the reheating epoch. The value of T reh can have an impact on the production of dangerous relics, such as gravitinos, or on the formation of topological defects from thermal phase transitions and the gravitational waves they generate. The expansion history of the universe during reheating, and in particular the moment when the equation of state approaches 1/3 can have important implications for the uncertainties in predictions of inflationary models [34,87].
Thermalization can be a long process, much longer than the preceding preheating, transient and turbulent phases. In principle, the universe can attain a radiation-like equation of state during or shortly after the turbulent stage, i.e., it can satisfy one of the two criteria for thermalization quite early. However, reaching a state of Local Thermal Equilibrium (LTE) can take much longer and involve particle fusion and off-shell processes. We say that the universe is in a prethermalized state if w ≈ 1/3, but LTE is not established yet. Prethermalization can be delayed by the formation of long-lived objects like oscillons and Q-balls. They behave as pressureless dust and therefore must decay into relativistic matter to achieve a radiation-like equation of state before BBN. Similar considerations apply to massive scalar field condensates. For instance, if there is some remnant inflaton condensate, even if subdominant in energy during thermalization, it can make the universe re-enter a matter-dominated state of expansion before BBN. To avoid this, one must ensure the complete decay of the condensate. Introducing perturbative decays through three-leg interactions like Yukawa couplings, hφψψ, proves to be a reliable way for the absolute removal ofφ. That is why, albeit unimportant during the early non-perturbative stages of reheating, perturbative decays ofφ are vital for the late stage of thermalization.
In a state of LTE the local value of the entropy, i.e., the entropy per unit volume, s, is maximized. LTE is achieved by particle species which are both in kinetic and chemical equilibrium. This requires both the re-distribution of momentum and energy between different particles, as well as an increase in their total number. Hence, both numberconserving and number-violating (off-shell process, particle fusion) reactions are involved. Negligible interactions between different species lead to Bose-Einstein and Fermi-Dirac distributions for bosons and fermions, respectively, in kinetic equilibrium. Kinetic equilibrium entails efficient exchange of energy and momentum between particles, i.e., it is sufficient to have number-conserving interactions only. On the other hand, chemical equilibrium can be achieved only by changing the number of particles. If number-violating interactions are suppressed (this could occur if the particles mediating the number-violating interactions acquire a large mass at early times) the state of kinetic equilibrium is also known as a quasi-thermal state. However, as number-violating processes become efficient and particles flow to lower chemical potentials until the sum of chemical potentials of reacting particles becomes equal to the sum of the chemical potentials of the products in every reaction, s can be truly maximized and full LTE reached.

Perturbative limit
If after the end of inflation Bose effects and non-adiabatic particle production are unimportant, i.e., the inflaton condensate undergoes perturbative decays as described at the beginning of Section 5.1, then the decay is completed when ρ ∼ Γ 2 m 2 Pl . In the trilinear model, V (φ, χ) = m 2 φ 2 /2 + σφχ 2 , the perturbative decay rate is Γ = Γ φ→χχ , see eq. (5.1), and to ensure no parametric resonance q = σΦ/m 2 σm Pl /m 2 1. Hence, ρ ∼ (σ 4 /m 2 )m 2 Pl m 4 (m/m Pl ) 2 m 4 . The momentum of the massless χ particles will be m/2, implying a particle energy E = m/2 ρ 1/4 . On the other hand, the χ particles number density is n = ρ/ E ρ 3/4 . Note that in LTE E LTE ∼ T and ρ LTE ∼ T 4 , implying LTE . Thus, perturbative preheating leads to a non-equilibrium dilute universe containing very energetic particles. To ensure the completion of thermalization, we should now introduce number-conserving and number-violating interactions, which are efficient even in a dilute plasma. When the universe reaches a radiation-dominated state of expansion, and the rate of these interactions is >H, reheating is completed.

Non-perturbative effects
Resonant and/or tachyonic decays of the inflaton condensate are also highly non-thermal processes. They yield non-equilibrium spectra, with peaks lying in instability bands n(k) ∝ e 2µ k t at the end of preheating. After the phase of exponential growth is terminated by back-reaction and re-scattering, and after a brief period of chaotic evolution of inhomogeneous field configurations, the spectrum of a scalar field typically relaxes into a continuous band, n(k) ∝ k −1 , going all the way to k → 0 and having an increasing UV cut-off, k c (τ ) = k c0 (τ 0 )(τ /τ 0 ) q 2 , as described in Section 6.3. After the front of the distribution, which dominates the energy budget of the universe, starts behaving quantum mechanically, n(k c (τ )) = O(1), the matter spectra should relax into Bose-Einstein and Fermi-Dirac (for fermions weakly coupled to the scalars) distributions. We can then use the results obtained from the classical field theory analysis to put a lower bound on the duration of reheating. Assuming the radiation-dominated state of expansion begins soon after the end of inflation, a(τ )/a(τ 0 ) = τ /τ 0 = ρ(τ 0 ) 1/4 /ρ(τ ) 1/4 and putting (k c (τ reh )/a(τ reh )) 4 ∼ T 4 reh (recall n(k c (τ reh )) = O(1) and ρ(k) ∼ (k/a) 4 n(k)) and k c0 /a(τ 0 ) ∼ m τ reh τ 0 we find that Putting m = 10 −6 m Pl and ρ(τ 0 ) 1/4 ∼ 10 15 GeV and q 2 = 1/7 [18,19] yields T reh ∼ 10 3 eV. This estimate gives an unacceptably low reheating temperature, implying that additional interactions, e.g., decays into fermions, become important before the non-linear evolution of the scalar fields drives them into thermal equilibrium. The calculation of the reheating temperature in such models with highly inhomogeneous scalar field configurations remains an open challenge.

Reheating and High-Energy Physics models
'Is the universe 'elegant', as Brian Greene tells us? Not as far as I can tell, not the usual laws of particle physics, anyway. I think I might find the universal principles of String Theory most elegant -if I only knew what they were.'

Leonard Susskind
Accelerator experiments such as the LHC have given us information about the governing particle theory up to O(10 TeV). This is many orders of magnitude below the highest reheating scale allowed by observations. Measurements of the CMB anisotropies [2] constrain r < 0.11, implying that the energy scale at the end of inflation, see eq. (4.61), must be V 1/4 end < 10 15 GeV. This is also the upper bound on T reh . At such high energy scales we could ignore GUT-mass particles (approximately) and stringy states, but have to include all other degrees of freedom. Hence, there is a huge theoretical uncertainty regarding the actual model of reheating.
Since there is a great ambiguity regarding the degrees of freedom and their interactions at the high energy scales relevant to reheating we can just focus on the particle content of the Standard Model for simplicity. Ignoring the effects from extensions that account for baryogenesis and the generation of dark matter, one can study the evolution of the Standard Model degrees of freedom during preheating, assuming they were all spectator fields during inflation. Furthermore, coupling non-minimally the Standard Model Higgs field to gravity allows it to play the role of the inflaton field. This model is known as Higgsinflation [88]. It is quite interesting since in it all couplings of the inflaton to the Standard Model degrees of freedom are known, allowing for a complete calculation of the thermal history of the visible universe. Studies of the non-perturbative preheating dynamics (in the linear approximation) [23,24] have shown that non-linear effects become important soon after the end of Higgs-inflation and a further detailed numerical investigation is required for the calculation of T reh .
The exploration of extensions of the Standard Model, motivated by e.g., supersymmetry and/or supergravity, can introduce many new degrees of freedom. Unfortunately, they come with new interactions and parameters, many of which are poorly constrained, if at all. The vast landscape of string theory is a good example of the level of theoretical uncertainty one has to deal with when building models of reheating. To make further progress in constraining the particle physics of reheating we should turn to observations. A determination of the exact model of inflation through observations of the CMB could give us some insight into the physical laws governing the dynamics of reheating. Another possibility is the detection of a reheating signal that cannot be mimicked by any inflationary model. The observational consequence of reheating are the subject of the next section.
Despite the fact that we do not know the exact particle physics model describing reheating, it is safe to say that in any realistic scenario there will be a large number of scalar fields, fermions, vector fields, and perhaps non-minimal couplings to gravity and operators that are suppressed below some energy cut-off scale. In the rest of this section we briefly discuss their effects on various aspects of the non-perturbative linear dynamics Figure 7. The instability chart featuring the real part of the Floquet exponent normalized by the effective inflaton mass (left) and the Hubble rate (right), characterizing the χ particle production rate in the Vanilla model of preheating, V (φ, χ) = λφ 4 /4 + g 2 φ 2 χ 2 , g 2 > 0 [89]. In FRW space-timē Φ ∝ a −1 and k ∝ a −1 , implying that co-moving modes do not flow across the chart as the universe expands unlike in Figs. 3, 4, 5. The resonance is virtually unaffected by the expansion of space.
of preheating by considering simple toy models. We also talk about some generic models of (p)reheating.

Scalar fields
Historically, non-perturbative effects during preheating were first studied in the context of scalar field dynamics [10][11][12]14]. In Sections 5.2 and 5.3 we also used scalar fields to introduce the concept of resonant particle production. We showed that oscillations of the scalar condensate induce a time-dependence in the effective mass of the daughter scalar fields. Depending on the strength of the interactions, the occupation numbers of the decay products can grow either gradually in narrow momentum ranges for small couplings or in bursts in broad momentum ranges for large couplings, but in both regimes exponentially fast. We showed that this kind of particle production can be understood qualitatively using Floquet theory and, in particular, if the inflaton is massive -in terms of the instability chart of Mathieu equation. In this section we just wish to point out that although it is often enough to approximate the inflaton potential during the oscillatory phase of preheating by Taylor expanding around its minimum to quadratic order, there are models with massless inflatons that do not fall into this category. A famous example is the Vanilla model, The interesting thing about this model is that the Floquet analysis provides a virtually exact description of the resonant particle production [89]. Essentially, a certain choice of field and time redefinitions can lead to the absorption of all terms containing the scale factor, a(t), and higher order derivatives of it (provided we ignore the small oscillations in a(t)). The expansion of space can be transformed away, converting the equations of motion into Minkowski form. In this case the background equation is strictly periodic, whereas the equation for the field fluctuations is known as the Lame equation. The instability chart of the Lame equation, see Fig. 7, describes quite accurately the instabilities in the daughter fields and there is no need for the introduction of any flow lines, unlike the case of the Mathieu equation.
As we discussed in Section 5.3, stochastic resonance arises when we include the background expansion of space (apart from the Vanilla model). The effective masses of the daughter fields vary quasi-periodically with time and are modulated by powers of a(t). The momentum range of unstable modes is increased, at the expense of decreasing the rate of particle production. The same phenomenon is observed if several scalar fields are oscillating at the background level. If any of the ratios of their frequencies is different from one the effective masses of the daughter fields rarely go to zero (non-adiabatic events are rare) and if the ratio is an irrational number, the motion is not periodic at all. When the number of the oscillating background fields is 1, we still get resonant particle production, see Section 5.3.3. We note that it is quite natural for a large number of scalar fields to acquire VEVs during inflation in supersymmetric models.

Fermions
Interactions of the inflaton with fermion fields is a natural thing to consider. As already discussed at the end of Section 6.4, they can have important implications for the last stage of reheating. Interactions of, e.g., Yukawa form, are needed to ensure that any massive remnants of the inflaton decay into pairs of fermions and anti-fermions at late times, making it possible for the universe to become radiation-dominated.
We should point out that while the inflaton condensate oscillates, the fermions acquire a periodically varying mass and this can lead to fermionic preheating [90]. The resonance is not as efficient as in the case with daughter scalar fields, since the Pauli exclusion principle enforces the occupation number of a given mode to be ≤ 1. Nevertheless, the resonance can excite a broad range of modes, enhancing the decay rate in comparison with the standard perturbative estimate. In supergravity models, the gravitino can be non-perturbatively produced during reheating. Thermal production can take place after that as well. The danger of overproducing this massive relic can put constraints on its interactions and T reh [91].

Gauge fields
Similar to scalars, gauge bosons can be resonantly amplified quite efficiently during preheating. If the inflaton is a gauge singlet, it can be coupled to gauge fields through conformal factors without violating the gauge invariance of the action. The first term can lead to very efficient resonant transfer of energy to the massless gauge fields [37] during preheating. The second term violates parity and can generate chiral gravitational waves [92]. The second term naturally arises in models where the inflaton is an axion, e.g., in Natural inflation, with W 2 (φ) ∝ φ. Axions are the Goldstone bosons, appearing whenever an axial symmetry is spontaneously broken. An axion possesses an almost exact shift symmetry, so it naturally couples to total derivative terms such as µνησ F µν F ησ with µνησ the totally anti-symmetric tensor. Note that the mass dimension of both terms in eq. (7.1) is > 4 and they must be suppressed by some energy cut-off. If the inflaton is charged under a gauge symmetry, the covariant derivative can give rise to novel types of interaction. For instance, if the inflaton is a complex scalar, charged under an Abelian U (1) symmetry, the kinetic term in the action 2) yields a term that couples the complex phase of the inflaton with the gauge fields, in addition to a g 2 φ 2 χ 2 type of term. It turns out that the two transverse components of the spatial part of the gauge field appear only in the final term and their evolution during preheating is identical to that of χ and has been studied extensively [28,[93][94][95]. However, the longitudinal spatial component of the gauge field, the complex phase of the inflaton and A 0 are all coupled through the 2g A (φ∂ µ φ * )A µ term. Due to the complexity of the interaction, their evolution during preheating used to be approximated or ignored [28,[93][94][95] until very recently. Our paper [32] provided the first accurate treatment of the resonant particle production of these degrees of freedom, taking into account the redundancy introduced by the gauge freedom. It showed that all approximate treatments were insufficient for capturing the dynamics.
Preheating of a U (1) gauge field can be applied to the generation of the observed large scale magnetic fields [96]. Conformal couplings like the ones in eq. (7.1) can generate strong magnetic fields soon after the end of inflation [97,98]. They can act as primordial seeds for the galactic dynamo mechanism which can amplify them to the observed values today [96]. On the other hand, parametric resonance in models with an electrically charged inflaton fail to produce strong enough seed fields [32,[93][94][95]. It is worth pointing out that at the high energy scales relevant to preheating, the U (1) symmetry of electromagnetism is unified with the weak force, making it necessary to consider the full electroweak gauge theory SU (2) × U (1) [32,94].

Non-minimal couplings to gravity
The inflaton and the rest of the matter fields can have non-minimal couplings to gravity, which can become important at the high energies relevant to inflation and preheating. The simplest interaction one can consider is of the form ξχ 2 R, where χ could be the inflaton or a daughter scalar field, R is the Ricci scalar and ξ is a dimensionless coupling constant. If it was the inflaton, than for field values m Pl / √ ξ the expansion of space will be affected by the interaction term, e.g., Higgs-inflation [88]. Otherwise, in general, since R oscillates during preheating, see eq. (4.22), this type of interaction provides a new way for amplifying scalar field fluctuations. If |ξR| ∼ |ξ|H 2 the oscillating effective mass squared of field fluctuations (the mass induced by non-gravitational interactions), the parametric resonance could be affected. And if the gravitational interaction provides the dominant contribution to the effective mass, than it can induce resonant particle production on its own. We note that models with inflaton interaction terms that are linear in R can be studied in the Einstein frame, which is related to the original frame, also known as the Jordan frame, via a conformal transformation. The conformal transformation automatically makes coupling constants time-dependent during preheating. The resonant particle production in the Einstein frame has been studied in [23,24] for the Standard Model degrees of freedom, in [99][100][101] for scalar field fluctuations after multi-field inflation and in [102] for gauge fields.

Non-conventional interactions
We refer to all operators that have some cut-off scale, Λ UV , below which they are suppressed, as non-conventional. They can be important during preheating. For instance the nonminimal coupling to gravity discussed in Section 7.4 falls in this category and becomes negligible for field values m Pl / √ ξ ≡ Λ UV . All terms of higher than 4 mass dimension also have a cut-off scale, e.g., the terms in eq. (7.1). If the inflaton is an axion, we can write W 2 (φ) = φ/Λ UV ; Λ UV should be associated with the Peccei-Quinn scale and the gauge fields with the gluons to resolve the Strong CP problem. Another example is the case of noncanonical kinetic terms, f (φ/Λ UV )∂ µ φ∂ µ φ, with lim x→0 f (x) = 1. They automatically arise in models with non-minimal coupling to gravity, see Section 7.4, after transformation to the Einstein frame. Similar patterns are observed in models with non-local interactions where the Fourier transformed kinetic terms are non-canonical f (k/k UV )|∂ τ φ k | 2 . The longitudinal component(s) of gauge field(s) in models with a charged inflaton have kinetic terms of this form with the role of the cut off played by the Compton wavenumber of the gauge field [32]. Another situation where suppression can occur is when high-derivative interactions are present (∂ µ φ∂ µ φ) 1+n /Λ n UV , e.g., DBI inflation [103]. A common feature of all non-conventional interactions is that they modify the effective mass of the field fluctuations (after canonical normalization) above the cut off scale. This typically changes the resonance structure -it alters the shapes of the instability bands, but never degrades the efficiency of the resonant particle production. In fact non-conventional interactions can provide an alternative channel for resonant preheating.

Miscellaneous
Inflation and the subsequent stage of reheating allow for the testing of low energy models of particle physics (e.g., supersymmetric models), constrained by colliders. Even if all fields are negligibly coupled to the inflaton during inflation and reheating, they can still exhibit non-perturbative preheating dynamics or even lead to new phenomenology. Fields that have negligible interactions with the inflaton sector are known as spectator fields.
Light spectator scalar fields during inflation (having masses H) develop an effective non-zero vacuum expectation value. Basically, the equation of motion for long-wavelength modes is overdamped during inflation. Once vacuum fluctuations cross outside the Hubble radius, they freeze and their amplitude is determined by the Hubble rate which is approximately constant, yielding a nearly scale-invariant power-spectrum The mechanism is similar to the one which generates the curvature perturbations. In fact, the fields evolve identically to the inflaton fluctuations in the spatially-flat gauge in the slow-roll approximation, see eq. (4.40).
As inflation ends the preheating dynamics of the spectator field χ is reminiscent of that of the inflaton. On small scales, comparable to or shorter than the Hubble radius at that time, the field can be approximately separated into background and inhomogeneous parts, i.e., χ(x) =χ + δχ(x), withχ drawn from a Gaussian distribution with variance 19 ∆ 2 χ , δχ(x) being vacuum sub-horizon fluctuations and |x| < H −1 . This is known as the separate universe approach. The homogeneous value of the spectator,χ, does not evolve until the Hubble rate becomes smaller than its effective mass. Then it starts to oscillate about the bottom of its potential.
If χ is coupled to other fields, e.g., if it is the Standard Model Higgs, which is coupled to the charged leptons, W and Z bosons, its oscillations can lead to resonant particle production [28], followed by a non-linear period [26,29], generating a stochastic gravitational wave background [27,104]. More generally, a complexχ, embedded within, e.g., a supersymmetric model, can lead to baryogenesis, according to the Affleck-Dine mechanism [105], and the non-linear dynamics following the resonant stage can involve the formation of Q-balls [106,107].
Ifχ eventually comes to dominate the energy budget of the universe (e.g., if it oscillates about a quadratic minimum with its energy being redshifted as ∼ a −3 ), it has to decay into radiation to be in agreement with the big-bang nucleosynthesis scenario. If we assume that χ = 0 on cosmological length scales, than the radiation will have the inhomogeneities of the field imprinted on it. This is the essence of the curvaton scenario [108][109][110]. In it, the final primordial density fluctuations are generated after inflation and depend on the physics during reheating. Observations of the CMB give a constraint on the combination of the initial fluctuations from inflation and the post-inflationary ones coming from the decay of the curvaton field, χ.
Another possibility for a light spectator, χ, to lead to the generation of primordial curvature perturbations after inflation is if it affects the decay rate, Γ, of the inflaton into other fields. The fluctuations in χ will lead to a spatial variation in Γ. Hence, the final primordial curvature perturbation imprinted on the decay products is a consequence of the spatially varying couplings and any initial fluctuations generated during inflation. This mechanism is known as modulated reheating [111] and was first discussed in the context of Superstring theory models.

Observational implications and signatures of reheating
'The recent developments in cosmology strongly suggest that the universe may be the ultimate free lunch.' Alan Guth Despite being a very important and phenomenologically rich period, reheating and the high-energy physics laws governing it are hard to constrain observationally. Just like inflation and all other epochs preceding recombination, reheating cannot be observed directly, since it is hidden by the opaque thermal baryonic plasma. Similarly to the case of inflation, one should look for observational signatures of reheating that survive thermalization and could be inferred from various cosmological measurements. Inflation predicts the stretch of microscopic quantum fluctuations to super Hubble scales, generating a superhorizon curvature perturbation, which is conserved and eventually imprinted on the CMB. Unfortunately, during the decelerating phase of reheating, co-moving modes re-enter the horizon and only the sub-horizon scales are affected by the non-linear dynamics. The length-scales on which the curvature perturbation is affected are so short, that the change is completely concealed by the later non-linear evolution of cosmic structure, making it impossible (for now) to be inferred from the CMB. Another reason why reheating is difficult to connect with observations is that by the time of BBN at the latest, all Standard Model species must be thermalized, hiding away the details of the earlier stages when they were produced.
Still, reheating can yield signatures, potentially observable in the future. These include the generation of relics and metric perturbations, which could be observed directly. In effect, the early universe takes the role of an accelerator for poor people, allowing us to probe roughly, yet freely, the fundamental physics at otherwise virtually inaccessible energy scales.
Indirect signatures are also possible. For instance, the mapping of co-moving modes between horizon exit during inflation and re-entry at later times depends on the entire expansion history between the two events. Thus, the confirmation of a particular model of inflation can give us information about reheating, e.g., constrain its expansion history. Or vice versa, a better understanding of reheating can reduce the uncertainties in predictions of simple inflationary models.
In the rest of this section we discuss different observational implications of reheating, starting with the indirect expansion history effect.

Expansion history of reheating and the CMB
The expansion history of reheating is largely uncertain. We only know that between the end of inflation and the time the universe thermalizated completely (i.e., achieved chemical and local thermal equilibrium), the mean equation of state is t th t end dt w(t)/(t th − t end ) ≡w int > −1/3. This implies various possibilities for N -the number of e-folds of expansion before the end of inflation, when the pivot scale crossed outside the Hubble radius, k = a H . N can take a range of different values, depending on the expansion history of reheating. Hence, the predictions of any model of inflation have an inherent uncertainty, due to the poorly constraint period of reheating. Before discussing inflationary observables in detail we consider the uncertainties in N . Given a co-moving pivot scale that has re-entered the horizon at late times, having some fixed physical wavenumber today, say k ,phys0 = k /a 0 , we discuss how N depends on the details of the inflaton potential,w int , ρ th and g ththe last two being the energy density and the number of relativistic degrees of freedom at thermalization.
We start with the free parameters where {q i } are the parameters entering the inflaton potential, i.e., where the last expression follows from eq. (4.18) and the discussion above it. The value of the inflaton at the end of slow-roll inflation, φ end , is to a very good approximation insensitive to the initial conditions and the inflationary dynamics, implying φ end = φ end ({q i }). 20 Hence . This and eq. (4.61) imply that the measured magnitude of curvature perturbation where we can determine one of the potential parameters, e.g., Note that for a single-parameter model, e.g., V = m 2 φ 2 /2, this implies a one-to-one correspondence between V and N (and likewise for derivatives of V ), since A s is known from observations. For a two-parameter model of inflation, e.g., V = Λ 4 tanh 2 (φ/M ), there is a one-parameter set of solutions, etc. To find N and V (given {q i − 1}) we need to match the pivot scale today to the time it left the horizon, i.e., k /(a 0 H 0 ) = (a H )/(a 0 H 0 ). Taking the log of both sides of this equality, one can show that making the only assumption that entropy, s ∼ gT 3 , is conserved, sa 3 = const, after thermalization, a > a th . The derivation is given in, e.g., [87]. Since ρ end = ρ end ({q i }), after substituting V from eq. (8.3) into eq. (8.4), we find that Reheating affects N only through the specific combination of quantities appearing in the last argument. Note that the g th dependence contains the information about the For simplicity we fix ρ th and a th , and vary only the mean equation of state of reheating, w =w int (reheating takes place while a th > a > a end ). We consider single-parameter models of inflation, e.g., V = m 2 φ 2 /2, for which N = N (w int ) and V = V (N (w int )) = V (w int ). This implies that choosing an energy scale of inflation, V , uniquely determines N andw int , unlike in multi-parameter models of inflation where we have additional degrees of freedom, see Fig. 9. For plotting purposes, we have approximated the inflationary and the dark energy stages as de Sitter expansions, and the two periods preceding and following radiation-matter equality as radiation and matter dominated, respectively. None of the conclusions depend on these simplifications. Note that H ≡ d(ln H −1 )/d(ln a) = 3(1 + w)/2.
time of thermaliztion, a th . Essentially, the conservation of entropy implies that g th = ((π 2 /30)g 4/3 0 T 4 0 /ρ th ) 3 (a th /a 0 ) 12 , where T 0 = 2.725 K and the effective number of relativistic degrees of freedom in entropy is g 0 = 43/11. In other words, N depends on a combination of the energy scale, ρ th , and the time, a th , of thermalization, as well as the mean equation of state,w int , of reheating (holding for a th > a > a end ). We depict this effect in Fig. 8 for a single-parameter model of inflation, e.g., V (m, φ) = m 2 φ 2 /2. For simplicity, we fix ρ th and a th (fixing these two quantities fixes g th ) and vary onlyw int . We plot the Hubble radius, H −1 , in orange and the physical wavenumber corresponding to the pivot scale, aλ , in blue. Both have some fixed values today, H −1 0 and a 0 λ , respectively. For plotting purposes, we approximate the dark energy dominated universe today and inflation as stages of de Sitter expansion, i.e, for a > a de and a < a end , H = H 0 = const and Figure 9. Same as Fig. 8, but for two-parameter models of inflation, e.g., V = Λ tanh 2 (φ/M ). The extra parameter in the inflaton potential introduces an additional degree of freedom, N = N (q 1 ,w int ) and V = V (q 1 , N (q 1 ,w int )) = V (q 1 ,w int ). This means that unlike in the singleparameter models of inflation, choosing an energy scale of inflation, V , allows for a range of N andw int . H = H = V /3/m Pl = const, respectively. We also assume w = 0 between radiationmatter equality and dark energy domination, a de > a > a eq , and w = 1/3 between thermalization and radiation-matter equality a eq > a > a th . Note that for the single-parameter m 2 φ 2 /2 inflation, {q i −1} ∈ ∅. Given that A s = 2.2×10 −9 , eq. (8.5) implies N = N (w int ). This is shown in the figure with the three orange solid, dashed and dotted lines corresponding to three different choices ofw int . Note that we also have V = V (N (w int )) = V (w int ), see eq. (8.3). Hence, in single-parameter models of inflation,w int uniquely defines V and N (provided ρ th and a th are fixed). This is a peculiar feature of single-parameter models. In multi-parameter models of inflaton there is a degeneracy, as we discuss below. Note that even in the single-parameter models of inflation, the uncertainty in the equation of state of reheating translates into an uncertainty in the energy scale of inflation. This could be turned the other way round -a possible confirmation of a single-parameter model of inflation with a given V uniquely determines N and hencew int if ρ th and a th are known (if they are not, it at least uniquely determines the combination of ρ th , a th andw int on which N depends). So pinning down the model of inflation could give us information about reheating.
We now repeat the analysis for a two-parameter model of inflation, e.g., V (M, Λ, φ) = Here the inflaton mass is not a free parameter; instead m = m(N , A s ). Light and dark blue: the 68 % and 95 % confidence level regions for n s and r from Planck and other data sets [2]. Λ 4 tanh 2 (φ/M ). This time, we have one additional degree of freedom, i.e., {q i − 1} = q 1 . For fixed ρ th and a th this implies N = N (q 1 ,w int ) and V = V (q 1 , N (q 1 ,w int )) = V (q 1 ,w int ). The dependence is depicted in Fig. 9. As mentioned above, unlike the singleparameter models, multi-parameter models possess additional degeneracy due to the extra parameters. It explains the additional lines in Fig. 9. Essentially, there is not a one-to-one correspondence between V and N , i.e., a particular V gives a range of N . Note that in two-parameter models, having fixed V and N , uniquely determinesw int . Conversely, one has to measure or calculate V andw int separately, to uniquely determine N in twoparameter models of inflation.
Hence, the uncertainty in the expansion history of reheating leads to an uncertainty in N . This has consequences for spectral observables such as n s and r, see eq. (4.61). The slow-roll potential parameters at the time of horizon exit are The reason for having one fewer q i is because we take ratios of the inflaton potential and corresponding derivatives in eq. (4.18). Substitution for However, the measured magnitude of curvature perturbation reduces the number of free parameters by one, e.g., q 1 = q 1 ({q i − 1}, N , A s ). The first two expressions in eq. an additional degree of freedom. This transforms the bars in the n s -r plane into bands as shown in Fig. 11 for Λ 4 tanh 2 (φ/M ) inflation.

Relics
The out-of-equilibrium dynamics of reheating may lead to the generation of various relics. These include stable ones, e.g., topological defects, that have not been observed yet, tantalizing ones, e.g., dark matter and primordial magnetic fields, for which there is some experimental evidence, and observed ones, such as the baryon asymmetry. In the following we briefly talk about each of these applications of reheating.

Baryon asymmetry
Models of baryogenesis try to explain the observed baryon-to-photon ratio where n b = n B − nB and n γ are the (net) number densities of baryons and photons, respectively. There are many high-energy physics models that explain the value of η with some dynamical mechanism, none of which is singled out by observational tests. It is indeed very interesting to try to connect baryogenesis to reheating. We should point out that since the net number of baryons in the late universe, certainly after BBN, is conserved, then n b ∝ a −3 . Furthermore, since after the epoch of electron-positron annihilation, which also happened around BBN, n γ ∝ a −3 , this means that η = const. At earlier times η was still conserved, apart from the moments when a relativistic particle species went out of thermal equilibrium and became non-relativistic -then it changed in a step-like manner by a factor that is determined by the ratio of the old and new number of relativistic degrees of freedom. However, what remained always constant in thermal equilibrium was the ratio n b /s which is ∼ η 0 . The Standard Model of particle physics ensures the net conservation of baryons and cannot explain why there are so many photons (or so much entropy) per baryon today, assuming that the universe started in a natural state with no baryons n B = nB = 0, or with no net baron number, n b = 0. To explain the puzzling observed asymmetry in the amount of matter and antimatter, one has to invoke physics beyond the Standard Model. It must allow for physical processes that meet the following three criteria, known as Sakharov's conditions [112]: (i) non-conservation of the baryon number, b; (ii) violation of C and CP invariance; (iii) departure from thermal equilibrium. While condition (i) is obvious -to generate a net baryon number starting from n b = 0 we need reactions that violate baryon number conservation -it is not enough. Condition (ii) is necessary to ensure different decay rates into baryons and anti-baryons. Consider a baryon number violating reaction X → Y + Z. The violation of C invariance ensures that the rates of the reaction and its charge conjugated counterpart are different i.e., the b violating process that creates more baryons than anti-baryons is not counterbalanced by its conjugate that creates more anti-baryons than baryons. Taking into account the helicity of the baryons, e.g., a reaction of the form X → q L + q L , CP violation ensures that Γ(X → q L + q L ) = Γ(X →q R +q R ) and Γ(X → q R + q R ) = Γ(X →q L +q L ).
Otherwise, the amount of produced left-handed baryons equals the amount of produced right-handed anti-baryons and vice versa, implying that the net baryon number does not change Γ(X → q L + q L ) + Γ(X → q R + q R ) = Γ(X →q L +q L ) + Γ(X →q R +q R ) , (8.9) despite C invariance being violated and (i). The reason for condition (iii) is slightly less obvious. It comes from the fact that the equilibrium number densities of particles and anti-particles depend on their chemical potentials. In thermal equilibrium µ q L = −µq L , etc (recall that baryons and anti-baryons can annihilate into photons). However, since the baryon number is not conserved by the interactions, µ q L = µq L = 0. Hence, particles and anti-particles in thermal equilibrium have equal number densities despite conditions (i) and (ii). Many high-energy physics models that satisfy the three criteria have been put forward to explain the observed baryon asymmetry. GUT-scale baryogenesis models rely on superheavy particles (with GUT-scale masses) decaying into baryons through C and CP violating reactions. As discussed in Section 5.5, such superheavy particles can be produced non-perturbatively after inflation, out of thermal equilibrium, making (p)reheating the ideal setting for these models. Leptogenesis, the generation of a number asymmetry, n l = n L −nL, between leptons and antileptons, can also account for the observed baryon-to-entropy ratio. The lepton number, l, can be converted into the baryon number via Standard Model sphalerons -transitions between degenerate topologically different SU (2) electroweak gauge field configurations -they become suppressed at temperatures < 300 GeV. Essentially, for sphaleron transitions (you can think of them as reactions) b − l = const, but b = const and l = const. Thus, if one starts with n b = 0, by the end of the transition it is converted into n b ∼ n b−l . This is an example of electroweak baryogenesis and preheating can provide a way for generating the initial lepton asymmetry. Another possibility is the Affleck-Dine baryogenesis mechanism. It involves a complex scalar field, X, that carries a baryon number and whose non-equilibrium dynamics does not conserve b. The field normally starts in a spatially homogeneous configuration in which invariance under C and CP may or may not be spontaneously broken. Then it evolves into a non-thermal configuration, in which the C and CP symmetries are spontaneously broken, with a final non-zero b, which is eventually converted into baryons. A version of this model where the inflaton plays the role of the scalar field is the subject of [79]. There we show that the non-linear dynamics of reheating can play an important role for the prediction of η.

Magnetic fields
Magnetic fields are abundant in our universe [96]. They have been observed in galaxies B 10−10 2 kpc ∼ 10 −5 G and galaxy clusters B 0.1−1 Mpc ∼ 10 −6 G. There is a (conservative) lower bound on the strength of magnetic fields with cosmic scale correlation lengths B >1 Mpc > 10 −17 G. While galactic fields can be accounted for by the amplification of seed fields via the dynamo mechanism [113], the origin of those seeds, as well as the largecorrelation-length fields that are unaffected by magnetohydrodynamic processes remains an open problem. It can be explained by a primordial magnetic field component. CMB observations have put upper bounds on it B prim 1 Mpc < 10 −9 G [114], whereas the seed amplitude needed for the dynamo mechanism is model and scale-dependent.
It is difficult to connect the causal non-linear dynamics of reheating with the large scale magnetic fields. However, the linear stage of preheating can provide the perfect setting for magnetogenesis. Low momentum magnetic field modes can be resonantly amplified [93] or undergo tachyonic instabilities [97,98]. The biggest challenge is to avoid backreaction of small-scale modes before low-momentum modes have been sufficiently amplified [32,115]. Tachyonic instability can be achieved quite easily. A conformal coupling of the form L Maxwell = −f (τ )F µν F µν /4 yields where A T k (τ ) = a(τ )f (τ )A T k (τ ) are the canonically-normalized transverse (Fourier) modes. f tends to 1 at late times, but if it is ∝ τ α earlier on, certain choices of α and the magnitude of f could lead to a successful magnetogenesis via tachyonic preheating [116].

Miscellaneous
In many dark matter models, the relic abundance is determined by the self-interactions of a thermalized dark matter sector. After the inflaton resonantly excites the Standard Model and dark matter degrees of freedom during preheating, the two sectors can attain different equilibrium temperatures -a phenomenon known as asymmetric reheating, which can be sensitive to the non-linear dynamics of reheating [117]. A detection of a temperature difference can put constraints on the inflaton mass and couplings [118].
The non-linear dynamics of reheating can lead to the formation of stable topological defects, see Section 6.2, for which there are no observational evidence [1]. Overproduction of such defects could overclose the universe or affect CMB anisotropies. In fact, CMB measurements provide the tightest constraints [1].

Metric fluctuations
Departures from the FRW universe described by matter and metric perturbations are at the heart of modern cosmology. Within current observational limits, an adiabatic curvature perturbation (a scalar mode) with Gaussian statistics can explain the measured CMB temperature anisotropies [1]. Furthermore, the detection of polarization B-modes generated by primordial tensor fluctuations in the metric (gravitational waves) is one of the main goals of the upcoming Stage-4 CMB experiments [119]. The linear and non-linear stages of reheating can give rise to gravitational waves, as well as entropic and non-Gaussian contributions to the curvature perturbation. Their non-detection constrains different reheating scenarios, as we discuss in the remainder of this section.

Gravitational waves
Shortly after it was appreciated that non-perturbative particle production during preheating can lead to the fragmentation of the inflaton condensate, it was shown that the non-linear dynamics can give rise to a stochastic gravitational wave background [120] in addition to the one generated during slow-roll inflation. Unlike the gravitational waves from inflation [121], whose origin is quantum mechanical and power-spectrum scale-invariant, the gravitational waves from reheating are sourced by the classical evolution of inhomogeneities on sub-horizon scales and their power-spectrum is strongly peaked around a single frequency. Typically, the frequency of the peak is determined by the fragmentation lengthscale, which can be estimated from the linear analysis of preheating. Taking into account the expansion of the universe between reheating and today, one can show, see [17], f 0 ∼ β −1 H br m Pl × 4 × 10 10 Hz , Ω GW,0 ∼ 10 −6 β 2 , where f 0 and Ω GW,0 are the peak frequency and gravitational energy density per logarithmic frequency interval normalized by the critical energy density, respectively. Both quantities are evaluated today. H br is the Hubble rate at back-reaction -the time when most of the signal is generated, and βH −1 br gives the physical wavelength of the excited mode causing the back-reaction on the condensate. Typically, β = O(10 −2 − 10 −3 ), thus for efficient preheating after GUT-scale inflation f 0 ∼ 10 10 −10 11 Hz and Ω GW,0 ∼ 10 −10 −10 −12 . These frequencies lie above the highest frequency ranges 10 3 − 10 4 Hz of planned gravitational wave detectors [122], see Fig. 12. Decreasing the back-reaction energy scale drives the peak frequency towards the observable range, but the small amplitude of the signal is outside the reach of any of the upcoming gravitational wave observatories. Figure 12. The gravitational wave energy density per logarithmic interval in frequency in the current universe [17]. The peaked thick curve on the right is what is expected typically from preheating soon after 10 15 GeV-scale inflation.
When the inflaton condensate fragments as a result of the resonant particle production of light scalar fields, the gravitational wave background can get tiny modulations on large scales [123,124]. In addition to the prominent peak corresponding to the fragmentation (sub-horizon) lengthscale, the gravitational wave power-spectrum features a small component on low frequencies, too. The latter is a consequence of the superhorizon scale-invariant power-spectrum of light degrees of freedom developed during inflation, see eq. (7.3), which can be interpreted as the light fields having non-zero vevs that vary between different Hubble patches as discussed in Section 7.6 (recall that individual, causally disconnected, patches evolve independently of each other). The preheating and subsequent non-linear dynamics, including the amplitude of the generated gravitational waves, can be sensitive to these vevs, leading to a super-horizon modulation of the stochastic gravitational wave background. When a light scalar remains a spectator during inflation and reheating, i.e., remains decoupled from the inflaton, similar effects are observed if it is allowed to decay nonperturbatively. This was shown for the gravitational waves produced out of the resonant decay of the Standard Model Higgs into W and Z bosons and their subsequent non-linear evolution [27], assuming no coupling between the inflaton and the Standard Model sector.
Stochastic gravitational wave backgrounds from reheating with additional features can be generated as a consequences of the formation of defects [125] and non-topological solitons [44,78,80]. Even if non-linear effects never become important during reheating, different expansion histories affect the spectrum of the gravitational wave background generated during inflation [126].

Non-Gaussianities
Scalar metric perturbations can also be generated during reheating. While the adiabatic curvature perturbations are unaffected by reheating, see eq. (4.54) and the subsequent discussion, the generation of an entropy (or isocurvature) perturbation, S, during reheating could modify the total curvature perturbation, R. A significant growth of super-Hubble modes of R occurs if on these scales ∆ 2 S ∆ 2 R [16]. However, in models with interacting fields (prone to peheating), e.g., V = λφ 4 /4 + g 2 φ 2 χ 2 /2, the super-Hubble power-spectra at the end of slow-roll single-field inflation are ∆ 2 S ∆ 2 R ∼ 10 −9 , where S = (H/φ)χ [15,66]. Even if the entropy perturbation is resonantly amplified during preheating, back-reaction takes place while on super-Hubble scales ∆ 2 S ∆ 2 R [43] and the observationally-relevant part of the power-spectrum of R is affected weakly (at most).
On the other hand, important statistical properties of the curvature perturbation, such as the bispectrum of its non-Gaussianities [127,128], can be affected significantly by the resonant entropy production during preheating and the subsequent non-linear dynamics. The same mechanism responsible for large-scale modulations in the gravitational wave spectrum from preheating (see the above discussion of gravitational waves) also leads to strong non-Gaussianities in the curvature perturbation [43]. Essentially, extreme sensitivity is shown to the vevs of light fields within individual Hubble patches by the expansion of these patches. The latter is equivalent to the curvature perturbation, implying, e.g., where R G is the standard nearly Gaussian adiabatic mode from single-field slow-roll inflation and the last term comes from the back-reaction and non-linear dynamics following the resonant amplification of the nearly Gaussian and scale-invariant (at the end of inflation) χ G , see eq. (7.3). As shown in [43], the transfer function F NL is highly non-linear and describes non-Gaussianities very different from the standard (weak) local ones The F NL term can in principle lead to non-Gaussian components that could be observable in the CMB [43] and their non-detection [129] constrains preheating scenarios. Other reheating scenarios that lead to potentially observable levels of non-Gaussianity include curvaton reheating (where curvature perturbations are generated by the decay of a slightly inhomogeneous curvaton field after inflation) and modulated reheating (where curvature perturbations are generated due to the dependence of the decay rate of the inflaton on the local value of a spatially varying field). Evolution of non-Gaussianity during reheating after multi-field inflation was studied in [130,131].
We should point out that the modification of large-scale curvature perturbations during preheating is consistent with causality, since it involves no transfer of energy across super-Hubble scales. Entropy perturbations are simply resonantly amplified and then converted into curvature perturbations. 21 Even when the reheating dynamics is perturbative and no strong resonances take place, the local non-Gaussianity prediction in single-field inflationary scenarios [132] f NL ∼ n s − 1 4 , (8.14) depends on the expansion history of reheating, see eqs. (8.5) and (8.6).

Afterword
By examining the dynamics in realistic models of reheating, we can tie together the wellunderstood and well-tested high-energy physics from laboratory experiments and the more speculative physics of inflation. We can determine not only how the Universe was populated with ordinary matter, but also the origin (and perhaps the nature and the fundamental properties) of cosmic relics such as dark matter, the baryon asymmetry, stochastic gravitational wave backgrounds, etc. The study of ever-more realistic models of reheating has been successful in recent years, and this looks set to continue. With the upcoming Stage-4 CMB experiments set to provide superb data to constrain further inflationary observables [119], understanding reheating-related uncertainties will become increasingly important for narrowing the range of viable models of inflation. We hope and expect that reheating will continue to be at the forefront of research in theoretical cosmology in the coming years.