The Secret Higgstory of the Highest Temperature during Reheating

We study the role of the Standard Model Higgs condensate, formed during cosmological inflation, in the epoch of reheating that follows. We focus on the scenario where the inflaton decays slowly and perturbatively, so that there is a long period between the end of inflation and the beginning of radiation domination. The Higgs condensate decays non-perturbatively during this period, and we show that it heats the primordial plasma to much higher temperatures than would result from the slowly-decaying inflaton alone. We discuss the effect of this hot plasma on the thermalization of the inflaton's decay products, and study its phenomenological implications for the formation of cosmological relics like dark matter, with associated isocurvature fluctuations, and the restoration of the electroweak and Peccei-Quinn symmetries.


I. INTRODUCTION
The discovery of the Higgs boson in 2012 not only provided the final piece of the Standard Model of particle physics, but also stimulated the realization that the Higgs is of fundamental significance to cosmology. Its deep connections to cosmic inflation, and in particular its dynamics during inflation, have attracted special attention [1][2][3][4][5][6][7].
If the Higgs is a light spectator field during inflation [8], then its quantum fluctuations accumulate on superhorizon scales and locally displace the field away from the minimum of its potential [9][10][11][12]. A Higgs condensate is formed, which does not survive in the universe today but rather is destroyed during the epoch of reheating, when the Higgs and inflaton transferred their condensate energy to a thermal bath of Standard Model particles.
In this paper, we present a new physical prediction for how reheating of the Standard Model proceeds when the inflaton field φ decays slowly, meaning that its decay rate Γ φ is small enough that Here m φ is the inflaton's mass, and M Pl = (8πG) −1/2 is the reduced Planck mass. If the inflaton decays, for example, through a dimension-5 Planck-suppressed operator L = (φ/M Pl )O 4 , then the decay rate is parametrically Γ φ ∼ m 3 φ /M 2 Pl andΓ ∼ 1. Models withΓ 1 then occur when the coupling that mediates the decay of the inflaton to matter is weaker than gravity [13].
A tiny coupling is technically natural, and couplings of gravitational strength or weaker for the inflaton are in * samuel.passaglia@ipmu.jp general desirable on theoretical grounds to avoid spoiling the flatness of the inflaton's potential [13,14], and on phenomenological grounds to avoid overproducing gravitinos during reheating [15].
In this regime, the decay of the inflaton is a perturbative process which can be understood analytically. As the inflaton decays to relativistic particles during reheating, it sources the primordial plasma with an energy density ρ r evolving according to until reheating is completed. We use the e-folds N ≡ ln a counted from the end of inflation as our time coordinate throughout this work, and H = H(N ) is the Hubble parameter. Under the assumption that the inflaton's potential is quadratic near its minimum, the inflaton condensate's coherent oscillations drive an effectively matter-dominated cosmological expansion with H ∝ e −3N/2 . The radiated energy density then has the asymptotic solution As the Hubble rate drops, this radiation loses energy in absolute terms but gains energy relative to the inflaton condensate itself. Reheating is generally taken to complete when Γ φ ∼ H, at a time where H end denotes the Hubble rate at the end of inflation. The plasma temperature at this time, i.e. the reheating temperature, is taken by convention as [16,17] T RH = 90 where g * is the effective number of relativistic degrees of freedom of the plasma. Since the radiation energy density decreases throughout the reheating epoch, the maximum temperature of the thermal component of the plasma can be much larger than T RH [16,17]. If the Higgs condensate is neglected, the maximum temperature during reheating is achieved when a majority of the radiation produced by the inflaton has thermalized, giving the relation [18] T φ-max 30ρ r π 2 g * 1/4 where the thermalization time N φ-therm can be computed as a function ofΓ [19]. Since the radiation energy density evolves as ρ r ∝ H ∝ e −3N/2 during matter domination, the maximum plasma temperature scales as T φ-max ∝ e −3N φ-therm /8 . We mark N φ-therm and T φ-max with φ's to denote they are associated with the decay of the inflaton.
In this work, we show that T φ-max is not necessarily the maximum temperature of the primordial plasma during reheating. The decay of the Higgs condensate formed during inflation can lead to a higher temperature.
After inflation the Higgs condensate oscillates around the minimum of its potential. We assume that the amplitude of the Higgs condensate is not so large as to probe the classically unstable part of its potential. Avoiding this instability requires either that the energy scale of inflation is below H unstable ∼ 10 10 GeV or that the Higgs potential is stabilized at high energies [3,20,21]. As the Higgs condensate oscillates in its approximately quartic potential, its energy density redshifts like radiation. Therefore, well after the end of inflation, the Higgs' energy density can be simply parameterized as whereρ h is a constant, dimensionless parameter which encodes all the details of how the Higgs condensate was formed and how it began to oscillate. We will show that it is generally of order unity or larger. The Hubble rate at the end of inflation, H end , will mainly scale all of our results together rather than affect the relative importance of the various processes. After a few oscillations, the Higgs condensate decays by parametric resonance in a completely Standard Model process [22,23]. This produces an effectively thermal plasmas at some time N h-therm not long after the end of inflation, with temperature Comparing the maximum temperature of the Higgs condensate's contribution to the plasma (8) to the maximum temperature of the inflaton condensate's contribution (6), we see that even though the radiation from the inflaton will eventually dominate the energy density of our universe, the radiation from the Higgs controls the maximum temperature of our universe if exceeds unity. We denote parametric relations with ∼ here and throughout. Withρ h generically of order unity or larger, and m φ ∼ H end as is typical for relatively large field inflation *1 we see that the Higgs contribution dominates over the inflaton condensate contribution so long as it sources a thermal population sufficiently before the inflaton does andΓ is small. DecreasingΓ both suppresses the energy density of the inflaton decay products and delays their thermalization time N φ-therm . Technical naturalness allows it to be very small, and empirically it is constrained only by ensuring that reheating completes before big bang nucleosynthesis around 1 MeV [24][25][26], yielding a lower bound, We will therefore be able to show that for much of the allowed parameter space, the maximum temperature of our universe is provided by the decay of the Higgs condensate after reheating, with our results joining a growing body of work highlighting the importance of the Standard Model Higgs in a wide variety of reheating scenarios: when the inflaton decays quickly [27,28]; when the inflaton decays to a hidden sector [29] or has a stiff postinflationary equation of state [30]; and when the condensate relaxation contributes to leptogenesis [31,32].
This article is organized as follows. Table I follows this introduction and summarizes the main parameters of this work and gives the equations in which they are defined.
In §II, we track the Higgs condensate from its formation during inflation, through its oscillations during reheating, to its resonant decay to effectively thermal Standard Model radiation. We estimate the parametersρ h and N h-therm which control the maximum temperature of the Higgs contribution to the primordial plasma.
In §III, we review how the inflaton condensate decays to a population of underoccupied hard particles, slowing down thermalization, to relate the inflaton-decay product thermalization time N φ-therm to the inflaton decay rateΓ. We compute for the first time the effect of the Higgs decay products on the thermalization of the inflation decay products.
In §IV we compare the Higgs and inflaton contributions to the plasma and show that for much of the parameter space the maximum temperature of our universe is controlled by the Higgs. We show that this temperature will generally have inhomogeneities on large scales *1 One would generally expect that if inflation ended on a quadratic potential H end /m φ ∼ φ end /M Pl .  uncorrelated with the adiabatic fluctuations sourced by the inflaton.
In §V, we discuss the theoretical and observational implications of having a maximum temperature controlled by the Higgs, which include the enhancement of the production of relic particles and the restoration of spontaneously broken symmetries. We conclude in §VI.

II. RADIATION FROM THE HIGGS
In this section, we discuss the formation of the Higgs condensate during inflation, its dynamics at the end of inflation, and its decay to effectively thermal radiation via parametric resonance after inflation.

A. Condensate formation
In their pioneering work, the authors of Ref. [12] studied the equilibrium state of a self-interacting scalar field in a de Sitter background. This is done by treating the field amplitude, coarse-grained on a fixed physical scale larger than the Hubble scale ∼ 1/H, as a random variable. Its evolution is governed by a competition between deterministic rolling and stochastic fluctuations of amplitude ∼ H/2π per e-fold as the exponentially-expanding vacuum fluctuations cross the averaging scale [33].
The probability distribution over field amplitudes can be calculated by finding the stationary solutions of a Fokker-Planck equation. For example in the case of four real scalar fields with an SO(4) symmetry in their quartic self-interaction λ ϕ 4 /4, the equilibrium distribution in de Sitter space with Hubble parameter H dS has moments ϕ = 0 and [21] where h ≡ | ϕ|. We say then that h forms a scalar condensate. This discussion carries over to the SU(2)-doublet Higgs field during inflation. Assuming that the inflationary Hubble scale is much larger than the electroweak scale v 246 GeV, the Higgs field will develop a scalar condensate. The typical condensate amplitude can be estimated in de Sitter using Eq. (11) upon identifying λ with the Higgs self coupling [3,34]. For simplicity, we neglect any running of λ and in numerical estimates take λ = 0.01. We assume throughout that the Higgs is minimally coupled to gravity and not directly coupled to the inflaton.
The cosmological inflationary epoch is only quaside Sitter, and the Hubble parameter decreases as inflation proceeds. Near the end of inflation in particular, the Hubble rate can evolve quickly enough that a condensate established stochastically early during inflation can further evolve and the de Sitter result (11) becomes an inaccurate estimate of the typical condensate amplitude [35,36]. This evolution can be tracked directly by solving the Fokker-Planck equation, but it can also be estimated by noting that it is dominated by deterministic rolling rather than stochastic fluctuations.
The local behavior of the condensate amplitude away from the origin, where we can ignore the effective angular momentum barrier associated with non-radial fluctuations, can then be understood from the Klein-Gordon equation of motion where w is the background equation of state; w −1 during inflation and w = −1/3 when inflation ends.
When the Hubble drag conditions |dw/dN | 1 and |d ln V ,h /dN | 1 are satisfied, the Klein-Gordon equation admits slow-roll solutions satisfying [37] but the drag conditions break down if the condensate amplitude is larger than   drag when h/H end ∼ λ −1/2 , and the Higgs displacement at the end of inflation is then independent of the Hubble rate earlier in inflation.

Initial Higgs Value
When this occurs the Higgs at the end of inflation will have a smaller displacement than the equilibrium distribution in the earlier phase of inflation would suggest, but a larger displacement than the equilibrium distribution with H dS = H end would imply. This is consistent with the results of Refs. [35,36] which solved the Fokker-Planck equation for a spectator field in a variety of inflationary backgrounds.
We show this effect in Fig. 1 by solving the equation of motion (12) with various initial conditions near the end of m 2 φ φ 2 inflation. In this inflationary model the Hubble rate at the end of inflation has decreased by an order of magnitude from the Hubble rate when CMB scales crossed the horizon, i.e. H end /H −60 1/12 and m φ /H end 2. Due to the dynamics discussed above, the condensate's amplitude at the end of inflation is roughly h end ∼ H end / √ λ, even if the initial displacement was much larger.
Therefore though the Higgs field amplitude at the end of inflation h end /H end is a stochastic variable, it should be within an order unity factor of the typical values where the first argument corresponds to the case where the Hubble rate does not decrease significantly near the end of inflation, and the second argument is the case where it does. H dS is now identified with the Hubble rate when the Higgs departed from the equilibrium solution. Since H dS /H end > 1 in any inflationary model, this window of typical expectations is fairly narrow in practice, at most a factor of ∼ λ −1/4 3 for λ ∼ 10 −2 .
After the end of inflation at N = 0, the Hubble rate drops and the condensate will eventually be released from Hubble drag and oscillate in its quartic potential. As seen in Fig. 1, in this regime the condensate energy density redshifts on the cycle average like radiation [38], which definesρ h as the dimensionless, time-independent asymptotic energy of the condensate. We calibrateρ h directly as a function of h end /H end by solving the equation of motion (12) numerically starting from the end of inflation with Hubble-dragged initial velocity and assuming that the inflaton potential is quadratic near the end of inflation and thereafter, V = m 2 φ φ 2 /2. From solutions in the range h end /H end ∈ [0, λ −1/2 ], we infer an empirical fitting function for the asymptotic Higgs energỹ The first factor of the first line accounts for the energy density carried by the Higgs field at the end of inflation, the second factor accounts for the Hubble drag phase until H ∼ √ λh end during which the energy density is approximately constant rather than redshifting like e −4N , and the final numerical factor 1.3 is calibrated from the numeric solutions. This fitting function agrees parametrically with the fitting functions provided for a wide range of post-inflationary backgrounds by Ref. [23].
For the typical Higgs displacements at the end of inflation given in (15), we therefore havẽ In our examples and to normalize scaling relations we takeρ h = 10 reflecting H dS ∼ H end and λ ∼ 10 −2 . Similarly the asymptotic Hubble rate after inflation ends and the inflaton oscillates on a quadratic potential can be approximated as where the coefficient is a fit to numerical results. The solution withρ h = 10 is shown in the first several e-folds after the end of inflation in Fig. 2. We compare this Higgs condensate energy density to the energy density of the inflaton and its decay products (3) for an inflaton decay rate withΓ = 10 −7 . While the Higgs condensate's energy density is much smaller than that of the inflaton, it can be large in comparison to the inflaton's decay products.
e-folds after inflation After inflation, the Higgs condensate's energy density (red,ρ h = 10) is negligible compared to the inflaton's (purple), but it can be larger than the energy of the inflaton's decay products (blue,Γ = 10 −7 ). The Higgs condensate (dashed) thermalizes (solid) some N h-therm ∼ 6 e-foldings after inflation, which can be well before the inflaton condensate's decay products thermalize at N φ-therm . The Higgs condensate thermalization time is associated with the condensate's decay at N h-decay (boxed region, shown in right panel), which takes some N h-decay ∼ 4 Higgs amplitude oscillations. We have assumed here that inflation ends in a m 2 φ φ 2 potential with m φ = 10 10 GeV, and we have not explicitly modeled the decay of the Higgs condensate here so ρ r,h continues to oscillate after N h-decay .

B. Condensate decay
After inflation, we assume that the Higgs condensate is on a purely radial trajectory, which we have checked is a good approximation due to Hubble friction. Rather than | ϕ|, h now represents a field that oscillates around zero and its equation of motion asymptotically approaches the simple form [38]ḧ where h ≡ e N h is the conformally rescaled field and overdots denote derivatives with respect to the conformal time η = dt/e N . This asymptotic equation of motion has a solution in terms of an elliptic cosine function [38] h h osc cn is the Euler gamma function. h osc is the amplitude of the conformally conserved oscillations, and η osc is an arbitrary turning poinṫ h(η osc ) = 0. Matching this solution to the asymptotic energy density (7) relates h osc toρ h as The Higgs condensate oscillates with this amplitude until it decays into Standard Model particles [22,23,[39][40][41] via parametric resonance [38,42,43]. The most efficient decay channel is to weak gauge bosons [23], since they have a large coupling to the Higgs field and they do not experience Pauli blocking [44]. Following Refs. [23,40], the Higgs decay to one such field A ∈ {Z, W ± } can be modeled using a set of 3 scalar fields χ i , one for each helicity, and a corresponding scalar-Higgs interaction L int = −g 2 A χ 2 h 2 /8 which is familiar from studies of inflationary preheating [38,42,43,45]. This modeling provides an qualitative understanding of the Higgs decay by parametric resonance analytically.
The equation of motion for the mode functions χ k in the oscillatory regime is [38] where the conformal mode function χ k ≡ e N χ k , the resonance parameter q ≡ g 2 A /(4λ), and denotes a derivative with respect to z. For a given energy scale, the Higgs self coupling λ and the various gauge couplings can be computed in the Standard Model, with typical values λ ∼ 10 −2 and g 2 Z ∼ 0.6, g 2 W ∼ 0.3 yielding q ∼ 10. At early times the Higgs condensate's energy loss from decay can be neglected, h follows the analytic solution (21), and the mode equation (23) describes a Lamé equation which exhibits resonant behavior and exponentiallygrowing solutions for certain ranges of k/( √ λh osc ) and q. When n(n + 1)/2 < q < (n + 1)(n + 2)/2 for an odd integer n, the first resonance band extends from 0 to k * ∼ ( √ λh osc )(q/(2π 2 )) 1/4 [38]. Modes in the resonance band grow exponentially in z and yield an exponentially increasing occupation number of particles with z osc ≡ z(η osc ) and with the Floquet exponent µ k a non-monotonic function of k and q bounded by 2 ln 1 + √ 2 /∆z 0.24. It can be approximated as a top hat with |µ| ∼ 0.2 [23]. The boson A is now identified with the weak gauge boson which yields the largest Floquet exponent, since it will dominate the condensate decay. Its exponentially increasing number density then yields an energy density where in the per-particle energy we have accounted for the effective mass g 2 h 2 /4 g 2 h osc /8 from the condensate displacement. A key timescale is z h-decay when the condensate has transferred an O(1) fraction of its energy to the gauge boson, ρ A = ρ h . This time is conveniently expressed in terms of a number of Higgs field oscillations, Using the fiducial values provided above, we estimate that the Higgs condensate decays after N h-decay ∼ 2 field oscillations. This number depends only logarithmically on all of the parameters in the problem except the Floquet exponent |µ|, and is independent of the energy in the Higgs condensate. *2 This simple analytic estimate that the Higgs decays in just a few oscillations is well supported by the lattice simulations performed by Refs. [22,23], which incorporate respectively abelian and non-abelian gauge structures for the Higgs. They find that the condensate decays most of its energy after N h-decay = 3 ∼ 4 oscillations.
To absorb the theoretical uncertainty in this number, we treat the number of Higgs field oscillations until decay N h-decay as free parameter in our reheating study, with the expectation that it lies in the above range.
Keeping N h-decay a free parameter also allows us to absorb the scenario where none of the Standard Model weak gauge bosons exhibit a k → 0 resonance. In this case, while there will always be a resonance band at some k and in the large q limit the resonance parameter µ k in that band will still approach the maximum possible value 0.24, in general the Floquet exponent µ k in a higher resonance band will be smaller and imply a correspondingly longer decay time N h-decay . Refs. [22,23] find that backreaction can take up to ∼ 4 times more oscillations than in the usual case. *2 One might wonder whether a resonance analysis is valid for such a small number of oscillations -just 4 zero crossings -but numerical solutions validate it within an intrinsic uncertainty of 1 zero crossing [38,46].
Finally, we need the number of cosmological e-folds that elapse between the end of inflation and the time of Higgs condensate decay. In the oscillatory regime the number of e-folds ∆N in an interval of ∆N oscillations is ∆N 2 ln (∆N ∆z/2 + 1) .
To get the absolute number of e-folds N in terms of the absolute number of oscillations N , we calibrate the mapping empirically by numerical solution of the condensate equation of motion. Identifying N as half the number of zero-crossings of the field, the number of e-folds to reach N h-decay for a condensate exiting inflation with dimen- where the first term accounts for the e-folds from the beginning of oscillations until decay (28), the second for the e-folds between the end of inflation and the beginning of the Higgs oscillations, and the final numerical term is calibrated from the numeric solutions. For N h-decay = 4, the typical decay time is where we have evaluated withρ h ∼ 10 and λ ∼ 0.01.

C. Effective Temperature and Thermalization
We have seen that in just a few oscillations of the Higgs field, most of the Higgs condensate's energy is transferred to weak gauge bosons. This radiation is very abundant in the sense that the mode occupancy per helicity is much larger than one. Specifically, combining Eqs. (24) and (27) yields an estimate of the gauge boson spectrum when the condensate decays. The occupation number is in a resonance band that extends up to a physical momentum The number density n * ∼ f * p 3 * and energy density ρ * ∼ f * p 4 * of produced particles is dominated by momenta p ∼ p * . This momentum is larger than the Hubble rate, which confirms that a particle description is valid. The weak boson radiation is already near thermal at the time of production. If the energy of the decay products were redistributed into a thermal spectrum with g * relativistic species, then the typical energy per particle would be set by the temperature T ∼ (ρ h, * /g * ) 1/4 . With E denoting the population average energy, comparing this expression to the per-particle energy of the produced particles yields and reveals that most of the energy is carried by particles with momenta p * ∼ T comparable to the eventual temperature of the thermalized system. With particles of similar number density and energy to those of a thermal distribution, the weak boson spectrum that results from parametric resonance is therefore effectively thermal. This observation motivates us to identify the effective thermalization time N h-therm with the Higgs decay time N h-decay , i.e.
At this time the effective temperature of the Higgs' decay products follows by energy conservation, where we have used g * ∼ 100 and scaled the result to typical values forρ h and N h-therm from Eqs. (18), (30), and (35). As we will see in §III, this is in sharp contrast with the situation for the inflaton condensate, which decays by producing less abundant particles but with much larger energy than that of a thermal distribution of the same energy density.
The subsequent evolution and complete thermalization of the effectively thermal plasma of Higgs decay products then depends on processes adjusting the number of particles and their momentum distribution. An important process while the Higgs is decaying is the non-linear interaction of the resonantly produced particles [47,48]. The non-abelian interactions of the Higgs' decay products are especially efficient at extending the overoccupied particles at p * to higher momenta with order unity occupancy, which is even closer to thermal [22]. At low momentum, scattering and absorption/emission processes likewise bring the distribution closer to thermal [49].

III. RADIATION FROM THE INFLATON
We now discuss the decay of the inflaton and the thermalization of its decay products. In contrast to the Higgs condensate's rapid decay via parametric resonance, which produces an abundance of effectively thermal particles, the inflaton's decay is perturbative, and its decay products first take the form of a severely underoccupied distribution of hard primaries. These slowly transfer their energy to a thermal soft population via in-medium splitting, delaying thermalization of the full energy released by the inflaton and lowering the maximum temperature of the inflaton's decay products.
Our calculations are based on the extensive literature on thermalization of non-abelian plasmas [50][51][52], applied to the cosmological context of reheating [18,19,49,53,54]. We provide the first calculation of how this thermalization process proceeds in the presence of the decay products of the Higgs condensate.

A. Inflaton condensate and decay
After inflation, the energy density of the universe is dominated by the coherent oscillations of the inflaton condensate. Radiation domination begins once the condensate fully decays. Well before that time, the condensate's energy loss is negligible, and it can be treated as a free field oscillating in a quadratic potential V (φ) = m 2 φ φ 2 /2 with a mass parameter m φ corresponding to the mass of an inflaton particle. On the cycle average the inflaton's energy density ρ φ ∝ e −3N corresponds to a matter-dominated universe H ∝ e −3N/2 as in Eq. (19). Each inflaton particle in the cold condensate carries an energy of ∼ m φ and their number density is If the inflaton decays at a rate Γ φ = m 3 φΓ /M 2 Pl into relativistic particles, then the energy density ρ r of the emitted radiation obeys The asymptotic solution is

B. Hard primaries
We assume that the inflaton decays to pairs of relativistic Standard Model particles, each with energy ∼ m φ /2. The energy of each hard primary produced by the inflaton redshifts as e −N . The number density of hard primaries obeys where the factor of 2 accounts for the pair of hard primaries produced by each inflaton particle decay. After some time their accumulated phase space distribution function is then [55,56] (also [18,19]) where g hard counts the decay products' redundant internal degrees of freedom (e.g., color and spin). The numerical coefficients ensure that the energy density of hard particles is just ρ hard = ρ r from Eq. (39), while the number density evaluates to In these calculations, the momentum integration is dominated by the UV cutoff p = m φ /2, and we drop terms suppressed by H/H end that are negligible after the end of inflation.
Though the hard particles are sourced continuously, the energy and number density of the hard particle population is always dominated by those produced in the most recent Hubble time. The per-particle energy ∼ m φ is then much larger than it would be in a thermal bath of the same energy T ∼ (ρ r /g * ) 1/4 ; explicitly, This quantity is large because the total decay product energy is suppressed byΓ 1 and decreases with time, while the energy of each particle is ∼ m φ which we assume is comparable to or larger than H end . Equivalently, the number density of particles in the hard population is much smaller than those in a thermal distribution of the same energy density.
The hard primaries are therefore very far from thermal. In order to thermalize, they must transfer their energy to an abundant population of soft particles.

C. Energy cascade toward thermal bath
The hard primaries emit lower momentum soft particles by collinear splitting in a medium comprised of the hard primaries themselves and the products of the previous splittings. This splitting leads to a cascade of energy from an underoccupied hard distribution to an abundant soft population. The total energy density of the radiation (hard + soft) obeys Eq. (39) until backreaction on the condensate occurs and reheating completes. For pedagogy and to connect to the existing literature, we neglect here the influence of the decay products of the Higgs. Their impact can change the thermal history of the plasma, as we discuss in §III E.
The abundantly populated soft particles rapidly thermalize at a temperature T soft that is much less than the hard particle energy m φ /2. Most of the energy ρ soft and particle number n soft in the soft population is carried by particles with momentum p ∼ T soft . T soft itself is determined by the energy the hard particles have lost through in-medium splitting. This transfer from hard to soft is a bottleneck that prevents immediate thermalization of the full energy released as the inflaton decays.
To write down the in-medium splitting rate, we must specify the nature of the hard particles and their interactions. For concreteness we consider a non-abelian plasma [50], and suppose that the hard primaries are gluons. Then g hard = 16 and the strong coupling is then denoted by α = g 2 s /4π ∼ 0.1. By emitting soft radiation, a gluon of momentum p in can split and form a gluon of momentum p out p in . The rate at which this splitting occurs in the medium is estimated as [54,57] Γ split (N, p in , p out ) ∼ αΓ el min 1, where Γ el is the rate for elastic gluon scattering. Γ split only depends on the incident gluon's momentum p in through the kinematic restriction p out p in , and Γ split is suppressed for high-momentum daughters p out > p LPM by the Landau-Pomeranchuk-Migdal (LPM) effect [58,59]. The LPM effect is the destructive interference of radiation produced from nearby scattering sites when the formation time for the radiation is long compared to the typical time between scatterings. In QCD this occurs for daughter particles with momentum p out above [60,61] (see also Refs. [51,62,63]) where m is the gluon's screening mass in the medium. We estimate which has contributions from the hard particles and from the thermal soft population, containing g QCD ∼ 88 quark and gluon degrees of freedom. The screening mass also acts as an infrared regulator for the elastic scattering rate At this point we remind the reader that all expressions containing the coupling α should be viewed parametrically. For example we neglect numerical factors like the quadratic Casimir, and we neglect logarithmic factors that appear when solving the exact system of Boltzmann equations describing collinear splitting in a non-abelian plasma.
Energy transfer between the hard primaries and the soft population relies on efficiently producing lowmomentum particles by splitting. This splitting is efficient if Γ split H. Particles with large momenta p out are more difficult to produce due to the LPM suppression, and their production is inefficient above a momentum scale p split where which implies Note that p split receives contributions from collisions with both the hard and soft particles, and in §III E we will include the contribution from collisions with the decay products of the Higgs. Since the splitting rate (44) does not significantly depend on the incident momentum p in , any secondaries produced with p out < p split can themselves radiate efficiently. When they do so they lose an order unity fraction of their energy, and in this way the energy cascades from the hard particles to the thermal bath. (50) The integral is dominated by the UV modes with p ∼ p split , Asymptotically p split is dominated by the contribution from the soft particles themselves, and substituting p soft split (49) and n hard (42) yields Here we have used that (n soft /g * ) 1/3 ∼ (ρ soft /g * ) 1/4 for a thermal population, and we have defined a dimensionless constant which is much less than unity for the parameters of interest. The asymptotic solution of Eq. (52) is Despite the cosmological redshifting, the thermal population's energy density increases with time. This growth continues until p split reaches the hard particle energy m φ /2 and ρ soft ∼ ρ hard ; this occurs e-folds after the end of inflation. At this point, the inflaton's decay products (hard + soft) fully thermalize. The subsequent temperature is then set by conservation of energy, and the plasma's maximum temperature is reached at the beginning of this fully thermal phase [18] T This is the maximum temperature of the inflaton decay products during reheating in the slow decay regime, neglecting the presence and influence of the hot plasma produced by the Higgs. In Fig. 3, which we discuss further in the next subsection, it corresponds to the maximum of the purple line.

E. Effect of Higgs decay products
The Higgs condensate's decay products contribute to the medium in which the hard primaries split, amplifying the hard particles' splitting rate and changing their thermalization history. Previous studies of thermalization during reheating, upon which §III C-III D are based, have neglected this Standard Model process.
As a simplifying assumption, we suppose that the decay of the Higgs condensate to an effectively thermal population of relativistic particles occurs abruptly at N h-therm . These decay products then provide scattering targets for the hard primaries, and the energy density of the soft radiation emitted by in-medium splitting ρ soft evolves subject to an extension of Eq. (51) where the splitting scale p split is now p split = p hard split + p soft split + p higgs split .
The first two terms appear in Eq. (49), and the new contribution from the effectively thermal decay products of the Higgs is obtained by including their contribution to the elastic scattering rate (47) in the definition of the splitting scale (44), yielding which is time independent. Since the soft particles will themselves interact and thermalize with the Higgs' decay products, ρ soft solved from Eq. (51) should now be viewed as the energy transferred from the inflaton to the thermal plasma with total energy ρ soft + ρ h . We show the solution for ρ soft for a typical parameter set in Fig. 3. We take an inflaton decay rateΓ = 10 −12 and an inflaton mass m φ = H end . We assume a Higgs energyρ h = 10 and a Higgs decay time N h-therm = 6, which are typical values as seen in (18) and (30). We use a coupling strength α = 0.1.
At early times, the contribution to the total plasma energy from the Higgs' decay products dominates over the contribution from all the inflation decay products. The splitting scale (58) is therefore dominated by p higgs split and so the solution of Eq. (51) is where the minimum in this expression limits splitting to occur below the energy of hard particles themselves. Due to the α 4 factor in (59), p higgs split is smaller than m φ unless H end m φ orρ h is much larger than unity. *3 Despite being enhanced by scatterings off the Higgs' decay products, ρ soft is initially a subdominant component of the thermal plasma relative to the Higgs contribution ρ h . The solution (60) is valid until ρ soft overtakes ρ h at time which is the intersection time of the red and blue lines in Fig. 3. N switch is an important time. Observables which depend on the temperature of the plasma during reheating, such as the dark matter relic abundance we will discuss in §V, inherit the fluctuations of the thermal plasma at the time they were produced. If they are produced before N switch , the Higgs contribution to the thermal plasma dominates and observables will inherit the Higgs' fluctuations. If they are produced after N switch , they inherit the inflaton's. We are therefore interested in the maximum temperature of the plasma before and after N switch . Before N switch , the Higgs contribution dominates the thermal plasma and the maximum temperature is T h-max (8). At N switch , we writeρ h in terms of p higgs split to compare the energy of the thermal plasma to the maximum energy it has when it thermalizes without the Higgs ρ r (N φ-therm ). We find and thus regardless of whether p higgs split is greater or less than m φ , by the time the inflaton contribution to the thermal plasma dominates over the Higgs contribution, the plasma is inevitably at a lower temperature than its maximum without the Higgs, T φ-max .
After N switch , the energy ρ soft transferred from the inflaton hard primaries dominates the energy density of the thermal plasma and therefore provides the dominant contribution to the splitting rate. If p higgs split was less than the hard particle energy m φ , as we expect and as shown in Fig. 3, then p soft split begins to grow as described in §III D and approaches the asymptotic solution (54). Once ρ soft reaches the hard particle energy ρ hard the plasma thermalizes completely. At this point it reaches the maximum temperature T φ-max computed without the effect of the Higgs in Eq. (57). *3 While this can in principle be the case if λ is very small (see Eq. (18)), in the very small λ regime the decay of the Higgs condensate by parametric resonance may be disrupted if the weak gauge bosons decay to fermions on the condensate oscillation timescale [46,64].
We therefore see that while the Higgs can increase the hard primary splitting rate to enhance the inflaton contribution to the thermal plasma, T φ-max still represents the maximum temperature of the plasma after N switch as long as p higgs split m φ .
If, however, p higgs split m φ , then the inflaton contribution to the thermal plasma only comes to dominate the Higgs one after the thermalization time N φ-therm (55) which defined the maximum temperature T φ-max of the plasma neglecting the Higgs. Explicitly, (63) When N switch > N φ-therm , the maximum temperature of the plasma after the inflaton contribution dominates at N switch is in fact determined by the energy at N switch itself ∼ ρ r (N switch ), lower than it was without the Higgs. We will not focus on this regime in the following.

IV. MAXIMUM TEMPERATURE DURING REHEATING
We have seen that the thermal plasma during reheating receives contributions from the decay of the Higgs and inflaton condensates, and that it can be much hotter than the ultimate reheating temperature at the onset of radiation domination. The inflaton decay products provide a maximum temperature (57) which is controlled by the inflaton parameters m φ /H end andΓ and by the Standard Model parameters of which we have retained here only α.
The decay of the Higgs condensate provides a maximum temperature (36), which is controlled by the Higgs density parameterρ h and by the Higgs effective thermalization time N h-therm . Both of these temperature scales can be much larger than the reheating temperature (5) T RH H end ∼ 10 −7 H end 10 10 GeV We show these temperatures in Fig. 4 as a function of the inflaton decay rateΓ for the parameter values above. With these parameter choices the maximum temperature of our universe is provided not by the decay of the inflaton but by the decay of the Higgs condensate provided that the inflaton decay rate is sufficiently small, or equivalently when reheating takes a sufficiently long time If the maximum temperature during reheating arises from the decay of the Higgs condensate, i.e. T h-max > T φ-max , then the plasma temperature will for some time track the Higgs condensate's spatial inhomogeneities. Since these inhomogeneities result from the Higgs' stochastic fluctuations during inflation, the maximum temperature during reheating can therefore locally inherit large-scale spatial correlations.
In de Sitter space, correlation functions must be invariant under the SO(4, 1) isometry. Consequently, the Higgs field's equal-time spatial correlation function G(R) can be computed from the unequal-time temporal correlation function G (∆N ) ≡ h( x, N ) h( x, N + ∆N ) . The leading order non-trivial behavior is [12,21] G(∆N ) ∝ e −|∆N |/Nc , with the asymptotic correlation time determined by the smallest non-zero eigenvalue Λ 1 of the Fokker-Planck equation for the Higgs' one-point probability distribution function. *4 The equal-time spatial correlation function between points separated by a physical distance R is then obtained by the mapping ∆N → 2 ln(RH). For regions probed by the CMB, which crossed the inflationary horizon some 60 e-folds before the end of inflation, the physical distance is R 60 ∼ H −1 e 60 , and we estimate the spatial correlation function to be for λ = 10 −2 . These estimates imply that the Higgs condensate during reheating was only partially correlated on large length scales that correspond to our present Hubble patch. If the Higgs condensate sets the maximum temperature of our universe during reheating, we expect order unity fluctuations in that temperature across the CMB scales.
Note that the evolution of the Hubble rate during inflation can change this expectation. As discussed in §II A, if the Hubble rate shrinks near the end of inflation then regions with stochastic Higgs fluctuations set earlier during inflation are released from Hubble drag and roll down the potential until their evolution is arrested once again. This causes regions with different values of the Higgs condensate to converge on the Hubble drag solution (see Fig. 1). In the H dS /H end 1 limit, this attractor behavior enhances the Higgs' correlation on large scales relative to Eq. (71), and can in principle reduce the inhomogeneity of the effects we will discuss in the next section.

A. Dark matter and unwanted relics
It may appear that by increasing the temperature of the primordial plasma the Higgs can source viable dark matter candidates when the inflaton alone cannot. This is not so. Any relic produced while the Higgs' decay products dominate the thermal plasma inherits the Higgs' order unity correlations on large scales (71), which can be viewed as an extreme form of matter-radiation isocurvature fluctuations. One should therefore instead be concerned that scenarios with a low reheating temperature could generically suffer from the production of catastrophic relics sourced by the Higgs contribution to the *4 We use here the eigenvalue for a real field with a Z 2 symmetry, because only radial steps decohere the condensate. The variance of the SU(2)-doublet Higgs is slightly larger relative to the radial steps than the variance of the Z 2 field, so the Higgs should be slightly harder to decohere and the correlation time (70) is a mild underestimate.
thermal plasma. Fortunately, we shall see that this can only occur if H end 10 11 GeV, a region of parameter space which is likely excluded by the Higgs instability constraint H end H unstable ∼ 10 10 GeV.
Observations of the CMB constrain the amplitude of matter-radiation isocurvature perturbations S mr to be 10 −1 R where R 5 × 10 −5 is the amplitude of the adiabatic perturbations [65]. Any stable, weakly-coupled relic χ that was produced by the thermal plasma in the epoch when it was dominated by the Higgs' decay products will carry O(1) energy inhomogeneities δρ χ ∼ ρ χ . If it is nonrelativistic at recombination it will behave as a dark matter component leading to an isocurvature amplitude S mr ∼ δρ χ /ρ tot comparable to its relic abundance Ω χ ∼ ρ χ /ρ tot . Consequently, the CMB isocurvature constraint implies an abundance constraint We now compute the relic abundance of particles produced during reheating by the thermal plasma. We separate the produced abundance into the contributions from the Higgs and inflaton decay products to see under what circumstances the bound (72) on the Higgs contribution is satisfied.
When the relic is produced during a matter-dominated reheating phase, its relic abundance can be related to its production time N * and its number density at production n χ (N * ) through [17] so long as after reheating the comoving entropy density is conserved. m χ is the mass of the particle, g * S,0 3.91 is the number of effective entropy degrees of freedom today, and H 100 = 100 km/s/Mpc. For simplicity we focus on thermal production mechanisms for n χ , though non-thermal production of relics during perturbative reheating can also generally be important due to the continual generation of hard primaries with energies ∼ m φ [66][67][68][69]. Since the Higgs provides additional scattering targets for those particles (see §III E), it may play a role in such processes as well.
The kinetic equation for n χ is [17] d(e 3N n χ ) where the thermally-averaged annihilation cross-section σv can be implicitly time dependent if it is temperature dependent, and n χ,eq is the equilibrium number density. We assume bosonic dark matter with one internal degree of freedom so n χ,eq = T 3 ζ(3)/π 2 .

Freeze out
Freeze out occurs when σv n eq ≈ 3H, and if this occurs when m χ T then n χ n χ,eq and the kinetic equation (74) yields where N FO is the freeze-out time. So long as the cross section is not enhanced at low temperatures, the particle annihilation occurs mainly near N FO and the number density at N N FO is simply where H FO is the Hubble rate when the freeze-out temperature T FO is reached and σv FO is the cross section at freeze out. The particular order unity numerical coefficient here assumes a temperature independent cross section. First let us consider the well-known case (see, e.g., Refs. [17,70]) in which the relic particle is produced by the inflaton decay products before radiation-domination is reached; this corresponds to T RH < T FO < T φ-max . Under these assumptions, the relic mass m χ , relic abundance Ω χ h 2 , freeze-out temperature T FO ≡ m χ /x FO , reheating temperature T RH , and scattering strength σv , are related through where we have parameterized the cross section as which defines the dimensionless parameter β. The value x FO ∼ 20 is familiar from studies of weak-scale freeze out, and x FO is only logarithmically sensitive to the mass scale (see, e.g., [70]). The validity of the underlying assumptions requires that the cross section be sufficiently high that the particle reaches equilibrium before freezout but not so high as to violate unitary bounds. Unitarity requires β 1, and thermal equilibrium requires σv FO n χ,eq,FO H FO , which translates to The upper and lower bounds on β then yield, respectively, lower and upper bounds on the m χ which can yield a given relic abundance at a given reheating temperature , subject to the condition T RH < T FO < T φ-max . Notice that unitarity provides a lower bound on the mass. This is because H ∝ T 4 before radiation domination; if instead freeze out occurs after radiation domination when H ∝ T 2 , then unitarity leads to the familiar upper bound on mass [71].
For Ω χ h 2 0.12 such that χ can be all of the dark matter, the lower bound on the mass is below the upper bound for all reheating temperatures below ∼ 400 TeV. The crossing point determines the maximum mass particle which can be the dark matter, m χ ∼ 10 10 GeV. On the other hand, in the limiting scenario of reheating right before BBN, T RH ∼ 1 MeV, the maximum mass for a dark matter candidate is only m χ ∼ 300 MeV.
Now if instead freeze out occurs at T φ-max < T FO < T h-max , then the Higgs contribution to the plasma determines the relic abundance and we generically expect order unity inhomogeneities in its number density on large scales.
In this case the temperature is related to the Hubble rate as H ∝ T 3/2 by Eq. (7) and Eq. (19), and so the abundance is enhanced for more massive particles. The analog of Eq. (77) *5 If the dark matter particles are relativistic at the time of freeze out, the resultant scaling would be the same as the lower limit on mχ shown here.
This limit sets the Hubble rate at the end of inflation above which the presence of the Higgs can be dangerous. If H end is above 10 11 GeV andΓ is suitably small but not too small (though the dependence is very weak), then the Higgs condensate can source relics by thermal freeze out in sufficient abundance to interfere with the CMB. Conversely if H end 10 11 GeV then no such relics can form.
Remarkably, this safe region is similar to the Higgs instability bound, which restricts the inflationary Hubble rate to be below H unstable ∼ 10 10 GeV. Thus if the Hubble rate during inflation is low enough that the Higgs was not sent to its instability, which we implicitly assume in our calculations, then by this criteria the Higgs is also not able to source dangerous isocurvature perturbations by thermal freeze out. The exact bound H unstable is on the other hand currently subject to experimental uncertainties in the top mass, as well as the stochastic history of the Higgs during inflation, and our results therefore provide an independent mechanism by which the Higgs may limit the inflationary energy scale.

Freeze in
The lower bound on the cross section (79) is the condition that the relic was once in equilibrium with the thermal bath. When relic particles are produced from the bath at a low enough rate that the process was never in equilibrium, thermal freeze in occurs instead [72][73][74][75][76][77][78][79].
With n χ n χ,eq , the kinetic equation (74) yields where we assume that m χ T (N ) while particle production occurs. Freeze-in production is inefficient in the complementary regime, m χ T (N ), where the rate is Boltzmann suppressed. The integral may be dominated by early times or late times, depending on how T , H, and σv vary in time.
With H ∝ T 3/2 the Higgs contribution comes from production at T h-max for all σv which do not decrease with increasing temperature. This UV-dominated freeze in yields a number density at production [80] where σv h-max is the production rate at T h-max and the specific order unity coefficient here assumes a constant σv .
With H ∝ T 4 the inflaton contribution, on the other hand, comes from production at the end of reheating for constant σv . This remains true for UV-enhanced cross sections σv ∝ T n for all n < 6 [81]. Only for n > 6 is the dominant contribution from T φ-max . Regardless, the relative relic abundance produced by the Higgs and inflaton depends on the relic mass and the relative values of σv at the respective production times. We focus here on the regime where σv is sufficiently larger at T h-max that the Higgs decay products provide the dominant contribution to the abundance.
That abundance is then bounded from above by the condition that it was not in thermal equilibrium, which is stronger than the relativistic version of the unitarity bound on the cross-section. As in the freeze-out case, the bound can be translated into a bound on H end given a relic density H end 10 11 GeV Ω χ h 2 10 −6 2/7 where we have scaled m χ by T h-max but note it should satisfy m χ T h-max for our approximations to be consistent. The danger zone for potentially producing isocurvature perturbations that are too large is quite similar to the freeze-out one (83): H end 10 11 GeV.

B. Symmetry Restoration
The maximum temperature of the thermal plasma during reheating has implications for the restorations of symmetries by finite temperature effects. The maximum temperature of the thermal plasma neglecting the presence of the Higgs condensate (64) can be written in terms of the reheating temperature as [18,19] T φ-max ∼ 100 GeV 10 10 GeV m φ 1/5 (88) Note that for a fixed reheating temperature it is inversely dependent on m φ . The electroweak symmetry is restored if the temperature reaches values 100 GeV, and we therefore see that, neglecting the Higgs condensate, for low reheating temperatures the electroweak symmetry is not necessarily restored by the thermal plasma if m φ 10 10 GeV. This point was discussed in Ref. [19].
The decay of the Higgs condensate, however, provides a plasma temperature (65) T h-max ∼ 10 7 H end 10 10 GeV independent of the reheating temperature. Comparing the maximum temperature from the Higgs decay to the maximum temperature from the inflaton decay, we see that if m φ ∼ H end then the decay of the Higgs is complementary to the decay of the inflaton: whenever the electroweak symmetry is not restored by the inflaton it is restored by the Higgs.
Thanks to the decay of the Higgs condensate, the electroweak symmetry is therefore in fact guaranteed to be restored in the early universe when reheating is perturbative, except when the Hubble rate at the end of inflation is low, The Higgs condensate may similarly play a role in the restoration of symmetries predicted by theories of new physics. For instance the hypothetical Peccei-Quinn axion arises as the pseudo-Nambu-Goldstone boson associated with the spontaneous breaking of a global U(1) PQ symmetry [82][83][84][85]. If the primordial plasma temperature initially exceeds the symmetry breaking scale, then as the universe cools the symmetry is broken in a cosmological phase transition which can lead to the formation of topological defects such as axion strings and domain walls. These can leave distinctive imprints on cosmological observables [86]. The symmetry breaking scale is model-dependent, and some of the most compelling scenarios have ∼ 10 10 − 10 12 GeV [87], but even if inflation occurs at a high energy scale the inflaton decay products may not be hot enough to restore such symmetries if the inflaton decay is slow [19]. The decay of the Higgs condensate, on the other hand, leads to a high temperature plasma irrespective of the reheating temperature and may facilitate such symmetry restorations and the accompanying phase transitions and topological defect formation. However, unless the Peccei-Quinn scale is sufficiently low, i.e.
T h-max from Eq. (89), even the Higgs condensate decay will be insufficient to restore the symmetry.
Note that the large scale inhomogeneity of T h-max does not directly impact relic defects from symmetry breaking except in rare regions where T h-max fluctuates to such low values that the given symmetry is not restored.

VI. CONCLUSION
We have studied the role of the Standard Model Higgs field during the epoch of reheating after inflation. The Higgs forms a condensate during inflation, and we have focused on its behavior in the scenario where the inflaton decays slowly, having a small perturbative decay rate Γ φ m 3 φ /M 2 Pl . We have calculated the maximum temperature achieved by the primordial plasma due to energy transfer from the Higgs and inflaton condensates, and the time evolution of this temperature between the end of inflation and the start of radiation domination.
To do so, we built upon well-established results in the literature and applied them to the Higgs, including the formation of a Higgs condensate during inflation from its quantum fluctuations as a light spectator field, the dynamics of the Higgs condensate after inflation as it oscillates on its quartic potential, and the parametric resonance in the electroweak gauge fields that induces an energy transfer from the Higgs condensate into an effectively thermal plasma which can be understood qualitatively analytically and quantitatively by lattice simulations.
We contrasted the resonant decay of the Higgs condensate with the inflaton's perturbative decay into rare hard particles which thermalize by an in-medium splitting suppressed by the LPM effect, in doing so drawing from the extensive literature on thermalization in nonabelian plasmas in the effective kinetic theory.
From this synthesis of constituent ideas, we are able to develop a comprehensive understanding of the Higgs condensate's role in reheating. The central conclusions of our work are: 1. Regardless of the inflationary history, the Higgs condensate has a typical energy density after inflation which lies within a narrow window between ∼ λ −2/3 H 4 end e −4N and ∼ λ −1 H 4 end e −4N controlled by the Hubble rate at the end of inflation.
2. The Higgs condensate is important to the thermal history of the universe when the inflaton decay rate is sufficiently small (67) 3. The maximum temperature of the primordial plasma is then obtained at the time of the Higgs condensate's fragmentation and decay, typically a few e-folds after the end of inflation, and for typical parameters in Eq. (89) is T h-max ∼ 10 7 GeV H end 10 10 GeV .
In the absence of a Higgs condensate, by contrast, the maximum temperature could be much lower, see Eq. (88). Fig. 3, the presence of the Higgs' decay products changes the thermalization history of the hard primaries produced as the inflaton decays by increasing the scattering targets in the plasma. This enhances the energy transfer from the hard primaries to the thermal sector, but once the inflaton contribution dominates the thermal energy of the universe the subsequent maximum temperature is set by the Higgs-less result above.

As shown in
5. The maximum temperature of our universe will inherit the Higgs' large scale stochastic fluctuations, which are uncorrelated with the curvature fluctuations in our universe. If relics are produced in substantial abundance from the plasma while it is dominated by the Higgs' decay products, they will therefore lead to unacceptable isocurvature fluctuations in the CMB. We show that for production by either thermal freeze out (83) or freeze in (87) this cannot occur so long as the inflationary scale is below H end 10 11 GeV.
6. Even for low reheating temperatures, the electroweak symmetry in our universe is restored by the decay of the Higgs condensate after inflation if H end 10 5 GeV, which complements symmetry restoration from the inflaton decay products. Peccei-Quinn axion symmetry may also be restored if its symmetry breaking scale is below T h-max .