Preheating after Higgs Inflation: Self-Resonance and Gauge boson production

We perform an extensive analysis of linear fluctuations during preheating in Higgs inflation in the Einstein frame, where the fields are minimally coupled to gravity, but the field-space metric is nontrivial. The self-resonance of the Higgs and the Higgsed gauge bosons are governed by effective masses that scale differently with the nonminimal couplings and evolve differently in time. Coupled metric perturbations enhance Higgs self-resonance and make it possible for Higgs inflation to preheat solely through this channel. For $\xi\gtrsim 100$ the total energy of the Higgs-inflaton condensate can be transferred to Higgs particles within $3$ $e$-folds after the end of inflation. For smaller values of the nonminimal coupling preheating takes longer, completely shutting off at around $\xi\simeq 30$. The production of gauge bosons is dominated by the gauge boson mass and the field space curvature. For large values of the nonminimal coupling $\xi \gtrsim 1000$, it is possible for the Higgs condensate to transfer the entirety of its energy into gauge fields within one oscillation. For smaller values of the nonminimal coupling gauge bosons decay very quickly into fermions, thereby shutting off Bose enhancement. Estimates of non-Abelian interactions indicate that they will not suppress preheating into gauge bosons for $\xi \gtrsim 1000$.

We perform an extensive analysis of linear fluctuations during preheating in Higgs inflation in the Einstein frame, where the fields are minimally coupled to gravity, but the field-space metric is nontrivial. The self-resonance of the Higgs and the Higgsed gauge bosons are governed by effective masses that scale differently with the nonminimal couplings and evolve differently in time. Coupled metric perturbations enhance Higgs self-resonance and make it possible for Higgs inflation to preheat solely through this channel. For ξ ≳ 100 the total energy of the Higgs-inflaton condensate can be transferred to Higgs particles within 3 e-folds after the end of inflation. For smaller values of the nonminimal coupling preheating takes longer, completely shutting off at around ξ ≃ 30. The production of gauge bosons is dominated by the gauge boson mass and the field-space curvature. For large values of the nonminimal coupling ξ ≳ 1000, it is possible for the Higgs condensate to transfer the entirety of its energy into gauge fields within one oscillation. For smaller values of the nonminimal coupling gauge bosons decay very quickly into fermions, thereby shutting off Bose enhancement. Estimates of non-Abelian interactions indicate that they will not suppress preheating into gauge bosons for ξ ≳ 1000.

I. INTRODUCTION
While the discovery of the Higgs boson at CERN [1] solidified our understanding of the Standard Model (SM), its behavior in the early Universe, above the electroweak symmetry-breaking scale, remains unclear. An intriguing possibility is the identification of the Higgs boson with the scalar field(s) necessary for driving inflation, the rapid acceleration phase of the Universe required to solve both the horizon and flatness problems, as well as seeding primordial fluctuations necessary for structure formation [2][3][4].
The original attempt to use the Higgs or a Higgs-like sector to drive inflation resulted in an inconsistently large amplitude of fluctuations [5], because of the value of the Higgs self-coupling λ in the Standard Model. However, the introduction of a nonminimal coupling between the Higgs field and the Ricci scalar can remedy this [6]. Such nonminimal couplings are not only generic, since they arise as necessary renormalization counterterms for scalar fields in curved spacetime [7][8][9][10][11][12][13][14][15][16], but they also grow without a UV fixed point under renormalization group flow, at least below the Planck scale [10]. The inherent ambiguity in the running of the Higgs self-coupling λ at high energies, due to our incomplete knowledge of possible new physics between the TeV and inflationary scales, leads to an ambiguity in the exact value of the required nonminimal coupling [17][18][19]. While simple estimates like λ infl ¼ Oð0.01Þ lead to the requirement of ξ ¼ Oð10 4 Þ, smaller values of λ can allow for much smaller nonminimal couplings. We will remain agnostic about the exact running of the Standard Model couplings at high energies and instead explore a broad parameter range 1 covering 10 ≲ ξ ≲ 10 4 .
A basic feature of inflationary models with nonminimal couplings is that they provide universal predictions for the spectral observables n s and r, largely independent of the exact model parameters and initial conditions [23,24]. These observables fall in line with the Starobinsky model [25] as well as with the large family of α attractors [26]. 2 1 For inflation on a flat plateau one should consider ξ ≳ 440 (e.g., Ref. [20]). In models of hilltop or inflection-point inflation, smaller values of ξ are possible, although UV corrections are expected to be larger. A discussion of the effects of UV corrections (as well as initial conditions) in Higgs inflation can be found, e.g., in Refs. [21,22]. In order to provide as complete a treatment of Higgs inflation as possible without referring to specific unknown physics, we choose to consider a broad range of nonminimal couplings that go below ξ ≈ 400. 2 See Ref. [27] for a way to alter the predictions of α-attractor models through multifield effects.
Even after the latest Planck release [28], these models, which predict n s ¼ 1-2=N Ã and r ¼ Oð1=N 2 Ã Þ, continue to be compatible with the data for modes that exit the horizon at N Ã ≃ 55 e-folds before the end of inflation.
While inflation provides a robust framework for computing the evolution of the Universe and the generation of fluctuations [3,[29][30][31][32][33][34][35], the transition from an inflating universe to a radiation bath (as required for big bang nucleosynthesis [36][37][38]), known as reheating, remains a weakly constrained era in the cosmic evolution. Despite the difficulty of directly observing reheating due to the very short length scales involved, knowledge of how the equation of state of the Universe transitioned from w ≃ −1 to w ¼ 1=3 is crucial, since it affects how one relates the observed cosmic microwave background (CMB) modes to the time during inflation when they exited the horizon [39][40][41][42][43][44][45][46]. This becomes increasingly relevant as new data shrink the experimental bounds on primordial observables.
The transfer of energy from the inflaton [which carries (almost) the entirety of the energy density of the Universe during inflation] to radiation degrees of freedom (d.o.f.) can occur either through perturbative decays or through nonperturbative processes. The latter case, denoted as preheating, includes parametric and tachyonic resonances (see Ref. [47] for a review). The end state of any (p)reheating scenario must be a universe filled with SM and dark matter (DM) particles, or at least intermediary particles that decay into the SM and DM sectors. Preheating therefore has the potential to address other long-standing challenges in cosmological theory (such as generating the observed baryon-antibaryon asymmetry [48][49][50][51][52]) or leave behind cosmological relics (such as cosmological magnetic fields [53,54] or primordial black holes [55][56][57]).
Higgs inflation provides a unique opportunity to study the transition from inflation to radiation domination, since the couplings of the Higgs inflaton to the rest of the SM are known. Detailed analyses of reheating in Higgs inflation were first performed in Refs. [58,59] and extended in Ref. [60] using lattice simulations. However, as discussed later in Refs. [61][62][63] and independently in Ref. [64], multifield models of inflation with nonminimal couplings to gravity can exhibit more efficient preheating behavior than previously thought, due to the contribution of the fieldspace structure to the effective mass of the fluctuations. Furthermore, it was shown in Refs. [61][62][63] that, in nonminimally coupled models, preheating efficiency can be vastly different for different values of the nonminimal coupling, even if these values lead to otherwise identical predictions for CMB observables. We will thus perform a detailed study of preheating in Higgs inflation, extending the results of Refs. [58,59,[61][62][63][64], in order to distinguish between Higgs inflation models with different values of the nonminimal coupling.
Because of the appeal of Higgs inflation as an economical model of realizing inflation within the particle content of the Standard Model, the unitarity cutoff scale has been extensively studied [19,[65][66][67][68] (see also Ref. [69] for a recent review). For large values of the Higgs vacuum expectation value (VEV) (like the ones appearing during inflation) the appropriate unitarity cutoff scale is M pl = ffiffi ffi ξ p , while for small values of the Higgs VEV it must be substituted by M pl =ξ, where M pl ≡ ð8πGÞ −1=2 is the reduced Planck mass.
In Sec. II of this work, we introduce a simplified model of a complex Higgs field coupled to an Abelian gauge field. Section III describes the generalization of this model to the full electroweak sector of the Standard Model. In Sec. IV we study the self-resonance of Higgs modes. Section V deals with the evolution of the gauge fields during and after Higgs inflation. At the end of this section we also address the unitarity scale. The decays and scattering processes that involve the produced Higgs and gauge bosons are described in Sec. VI, and observational consequences are described in Sec. VII. Concluding remarks follow in Sec. VIII.

II. ABELIAN MODEL AND FORMALISM
We build on the formalism of Ref. [70] for the evolution of nonminimally coupled multifield models, as it was applied in Refs. [23,71,72] during inflation and in Refs. [61][62][63] during preheating. The electroweak sector consists of a complex Higgs doublet, expressed using four real-valued scalar fields in 3 þ 1 spacetime dimensions: where φ is the background value of the Higgs field, h denotes the Higgs fluctuations, and θ; ϕ 3 , and ϕ 4 are the Goldstone modes. We also add the SUð2Þ L and Uð1Þ Y gauge sectors. We will start by closely examining an Abelian simplified model of the full electroweak sector, consisting of the complex scalar field and a Uð1Þ gauge field only. The full equations of the Higgsed electroweak sector are given in Sec. III, where we also discuss their relation to the Abelian simplified model. In order to connect our notation to that of Ref. [70], we identify ϕ 1 ¼ φ þ h and ϕ 2 ¼ θ. We will start by deriving the equations of motion for general ϕ I fields for notational simplicity. We use uppercase latin letters to label fieldspace indices, I; J ¼ 1, 2, 3, 4 (or just I; J ¼ 1, 2 in the Abelian case); greek letters to label spacetime indices, μ, ν ¼ 0, 1, 2, 3; and lowercase latin letters to label spatial indices, i, j ¼ 1, 2, 3. The spacetime metric has signature ð−; þ; þ; þÞ.
We first consider Uð1Þ symmetry with the corresponding gauge field B μ . The Lagrangian in the Jordan frame is given by The covariant derivative∇ μ is given bỹ whereD μ is a covariant derivative with respect to the spacetime metricg μν and e is the coupling constant. The corresponding field strength tensor 3 is By performing a conformal transformatioñ the action in the Einstein frame becomes and as in Refs. [70,71]. The potential in the Jordan frame is the usual Standard Model Higgs potential where the Higgs vacuum expectation value v ¼ 246 GeV can be safely neglected at field values that arise during inflation and preheating. Hence, the Higgs potential can be adequately modeled by a pure quartic term. For the sake of readability, we will drop the arguments of G, V, and f from now on. Varying the action with respect to the scalar fields ϕ I , the corresponding equation of motion for ϕ I is We work to first order in fluctuations, in both the scalar fields and spacetime metric. The gauge fields have no background component, and thus we only treat them as first-order perturbations. We consider scalar metric perturbations around a spatially flat Friedmann-Lemaître-Robertson-Walker metric, where aðtÞ is the scale factor. We may always choose a coordinate transformation and eliminate two of the four scalar metric functions that appear in Eq. (12). We work in the longitudinal gauge, where BðxÞ ¼ EðxÞ ¼ 0. Furthermore, in the absence of anisotropic pressure perturbations, the remaining two functions are equal, AðxÞ ¼ ψðxÞ.
We also expand the fields, Note that for Higgs inflation only ϕ 1 has a background value, φðtÞ, whereas the background value of ϕ 2 is zero. We may then construct generalizations of the Mukhanov-Sasaki variable that are invariant with respect to spacetime gauge transformations up to first order in the perturbations (see Ref. [61] and references therein): where we used Eq. (2) and we stress again that F νμ is defined using covariant derivatives. Until now we have worked in full generality and not chosen a gauge. Hence, we are in principle working with more d.o.f. than needed. We will distinguish two frequently used gauges: unitary and Coulomb gauge. The equation of motion of X h is unaffected by the gauge choice.

A. Unitary gauge
In the unitary gauge θ ¼ 0 ¼ X θ . Equation (23) thus becomes a constraint equation, The equations of motion for the gauge fields are rewritten as Separating the time and space components, the equation for Performing the analysis in Fourier space, with the convention fðxÞ ¼ R d 3 k ð2πÞ 3=2 f k e −ik·x , we derive an algebraic equation for B 0;k , The equation of motion for the spatial components B i is Using the constraint (26), going to Fourier space, and multiplying by a 2 , the equation of motion becomes̈B which is somewhat simplified in conformal time, We now distinguish between transverse (B AE k ) and longitudinal (B L k ) modes: The equations of motion for the transverse and longitudinal modes become

B. Coulomb gauge
In the Coulomb gauge (∂ i B i ¼ 0), the Goldstone mode θ remains an explicit dynamical d.o.f., and thus the relevant equations of motion are Going to Fourier space, we can solve for B 0;k in terms of θ k , similarly to the situation in the unitary gauge, By plugging the longitudinal mode into the last equation of Eq. (36) and demanding that it is zero, we get the additional constraint Substituting the two previous relations into the equation of motion for X θ gives We must demand that physical observables are identical in the two gauges, and derive a relation between θ k in the Coulomb gauge and B L k in the unitary gauge. X h k and B AE k are already identical in the two gauges. The longitudinal component of the electric field 6 is given by In the unitary and Coulomb gauge we get Since E L should not depend on the gauge, we can use these expressions to solve for B L in terms of θ. We obtain It is a straightforward algebraic exercise to show that by using Eq. (42), the equations of motion for B L k and θ k can be transformed into each other, providing a useful check for our derivation.
During preheating, when the background inflaton field oscillates, the unitary gauge becomes ill defined at the times where φðtÞ ¼ 0, as can be seen for example in the transformation relation of Eq. (42). We will compute gauge boson production during preheating by solving the corresponding linearized equations of motion in the Coulomb gauge, which is always well defined. The evolution of the background Higgs field φðtÞ will be numerically computed in an expanding background according to Eqs. (15) and (16), neglecting backreaction from fluctuations and produced particles.

C. Single-field attractor and parameter choices
For Higgs inflation, the function fðΦ; Φ † Þ is given by [6] fðΦ; For typical values of Higgs inflation λ ¼ Oð0.01Þ, and correspondingly ξ ∼ 10 4 . If we consider a different renormalization group flow for the self-coupling λ through the introduction of unknown physics before the inflationary scale, λ will become smaller or larger at inflationary energies. As we will show below, since the combination λ=ξ 2 is fixed by the amplitude of the scalar power spectrum, a larger or smaller value of λ during inflation will lead to a correspondingly larger or smaller value of the nonminimal coupling ξ. We will consider values of ξ in the range 10 ≤ ξ ≤ 10 4 . The inflationary predictions for the scalar and tensor modes for nonminimally coupled models with ξ ≥ 10 fall into the large-ξ single-field attractor regime, as described, e.g., in Ref. [23]. This results in very simple expressions for the scalar spectral index n s , the tensor-toscalar ratio r, and the running of the spectral index α as a function of the number of e-folds at horizon crossing N Ã The values for the spectral observables given in Eq. (44) correspond to single-field background motion. Multifield nonminimally coupled models of inflation at large ξ show a very strong single-field attractor behavior. The strength of the attractor was analyzed in Ref. [71] for the case of an SOðNÞ-symmetric model, similar to Higgs inflation without gauge fields. The more general case of two-field inflation with generic potential parameters was studied in Refs. [61,72], where it was shown that the single-field attractor becomes stronger for larger ξ and that it persists not only during inflation but also during the (p)reheating era. For generic initial conditions, the isocurvature fraction β iso is exponentially small for random potentials, while for a symmetric potential β iso ¼ Oð10 −5 Þ, as was shown in Ref. [72]. As discussed in Sec. VA, during inflation the gauge bosons are very massive compared to the Hubble scale, making the single-field attractor behavior of Higgs inflation stronger than the one described in Ref. [71] for the scalar symmetric case. Hence, the use of a single-field motion φðtÞ for the background is well justified during and after Higgs inflation. The dimensionless power spectrum of the (scalar) density perturbations is measured to be Using the tensor-to-scalar ratio from Eq. (44) with N Ã ¼ 55 yields r ≃ 3.3 × 10 −3 , and hence the tensor power spectrum becomes Given that the Hubble scale during inflation is approx- the Higgs self-coupling and nonminimal coupling must obey the relation We keep the value of the Hubble scale fixed and determine the value of λ that corresponds to each ξ through Eq. (48).

III. ELECTROWEAK SECTOR
We now consider the full SUð2Þ × Uð1Þ gauge symmetry, as it exists in the electroweak sector of the SM. The Lagrangian in the Jordan frame is given by with the Higgs doublet The covariant derivative∇ μ is given bỹ where Y is the generator of hypercharge Uð1Þ and B μ is the corresponding gauge field. The Higgs doublet has hypercharge þ1. We have also introduced the vector notation The A μ are the gauge fields corresponding to SUð2Þ and τ i are the Pauli matrices. The corresponding field-strength tensors are Defining the fields W μ , W † μ , Z μ , and A μ as the components of the covariant derivative of Φ are given bỹ The structure of the equations is almost identical to the one studied in the earlier parts of this work, where we focused on the Abelian case. For θ, the Goldstone mode that becomes the longitudinal polarization of the Z boson, we substitute At quadratic order, the field-strength term for the electroweak case is no more complicated than the Abelian case; it simply contains more fields, with Comparing Eqs. (56) and (62) with the Abelian case, we can easily find the equations of motion for the gauge fields.
The photon A μ does not couple to the Higgs: The Z boson obeys and correspondingly for the W AE bosons

A. Unitary gauge
In the unitary gauge, the equations of motion for the three Goldstone d.o.f. θ and ϕ 3 ; ϕ 4 give the constraints The equations of motion for Z μ and W μ are These equations are identical to the equations in the Abelian case, albeit with different couplings. The equations for the longitudinal and transverse modes are thus given by where W AE denotes the AE polarization of the field W (so the AE does not distinguish W or W † ).

B. Coulomb gauge
The Coulomb gauge for the two types of bosons, Z and W AE , is defined through the conditions In Fourier space, we can express Z 0;k in terms of θ and W 0;k in terms of ϕ 3 and ϕ 4 : From the decoupling of the longitudinal modes from the equations for the corresponding transverse ones we get the constraints Substituting into the equation for X θ gives and substituting into the equation for X ϕ 3 gives and likewise The equations of motion of the transverse modes of the Z and W arë

IV. HIGGS SELF-RESONANCE
We now focus on the Higgs fluctuations, neglecting the effects of Goldstone modes and gauge fields. In our linear analysis the Higgs fluctuations do not couple to the gauge field. The equation of motion for the rescaled fluctuations where the effective frequency is defined as For notational simplicity and connection to earlier work [61][62][63] we define the various contributions to the effective mass of the Higgs fluctuations, where M h h was defined in Eq. (19) and For the case of fluctuations along the straight background trajectory (like Higgs fluctuations), the Riemann contribution m 2 2;h vanishes identically. As described in Ref. [75] and further utilized in Ref. [61], the mode functions can be decomposed using the vielbeins of the field-space metric. In the single-field attractor-which exists in the nonminimally coupled models that include Higgs inflation both during [23] and after inflation [61]-the decomposition of X h k into creation and annihilation operators is trivial, Since the vielbeins obey the parallel transport equation D τ e h 1 ¼ 0, the equation of motion for the mode function v k becomes We solve the equation in cosmic rather than conformal time, which is better suited for computations after inflation,̈v where the frequency is defined in Eq. (80). We examine the two dominant terms of the effective mass: the one arising from the potential (m 2 1;h ) and the one arising from the coupled metric perturbations (m 2 3;h ). The latter is often overlooked in studies of preheating, perhaps because it is vastly subdominant during inflation. It arises by combining the equation of motion for δϕ and the metric perturbation ψ, defined through Eq. (12), in conjunction with the definition of the Mukhanov-Sasaki variables, given in Eq. (14).
The expression for m 2 1;h is where we used ξ ≫ 1 in expressions such as ð6ξ þ 1Þ ≃ 6ξ. Furthermore, since we are initially interested in studying the behavior during inflation where analytic progress can be made, we use ξφ 2 ≫ M 2 Pl as an approximation. As we will see, this works reasonably well even close to the end of inflation. We normalize the effective mass by the Hubble scale, We can use the single-field slow-roll results where we went beyond lowest order in ξφ 2 and we measure the number of e-folds from the end of inflation, meaning that negative values correspond to the inflationary era. 7 This leads to If we minimize m 2 1 as a function of δ ¼ ffiffi ffi ξ p φ, the field amplitude that minimizes the mass is 7 We neglected the contributions coming from the lower end of the integral leading to Eq. (91).
or, equivalently, N min ≃ −1.5. For the minimization we used the full expression for the effective mass and only took the Taylor expansion for large ξ at the end. We can see that for ξ ≫ 1 the minimum of m 2 1 is independent of ξ and thus occurs at the same value of δ (which will also be the same value of N) in the approximation of Eq. (91). In general, the function m 2 1;h ðNÞ=H 2 shows no appreciable difference for different values of ξ ≫ 1 during inflation. This can be easily seen by substituting Eq. (91) into Eq. (90). As shown in Ref. [62], this behavior persists during the time of coherent inflaton oscillations.
The mass component arising from the metric fluctuations is where the last approximation holds during inflation. Using the slow-roll expression for _ φ we get that during inflation This contribution is clearly subdominant to m 2 1;h , and hence it can be safely neglected during inflation. However, jm 2 3;h j grows near the end of inflation, since it is proportional to _ φ 2 , which at the end of inflation is given by It has been numerically shown in Ref. [61] that the field value at the end of inflation is ffiffi ffi Numerically, we get m 2 3;h =H 2 ðtÞ ≃ −11 at the end of inflation, in rough agreement with the approximate expressions given above.
The numerical results for ξ ¼ 10 are shown in Fig. 1, along with the approximate analytical expressions that we derived. We only show the ξ ¼ 10 case, since all cases with higher values of the nonminimal coupling exhibit visually identical results. After the end of inflation the two dominant components of the effective mass of the Higgs fluctuations evolve differently for different values of ξ. In Ref. [62] the behavior of m 2 1;h was analyzed in the static universe approximation. It was shown that for ξ ≳ 100 the effective mass component m 2 1;h quickly approaches a uniform shape regardless of the value of ξ. The consequence of this is that the Floquet chart for the inflaton self-resonance also approaches a common form for ξ ≳ 100. This can be seen in the left panel of Fig. 2, where m 2 1;h is very similar between ξ ¼ 100 and ξ ¼ 10 3 , but different for ξ ¼ 10. The coupled metric fluctuations component of the effective mass has a similar shape for ξ ¼ 100 and ξ ¼ 10 3 , but for ξ ¼ 10 it is significantly less pronounced, as seen in the right panel of Fig. 2.

A. Superhorizon evolution and thermalization
An important notion when dealing with (p)reheating is the transfer of energy from the inflaton condensate to the radiation d.o.f. Naively, one must compute all the power concentrated in the wave numbers that are excited above the vacuum state (meaning that they are different than the adiabatic vacuum at any time) and compare that to the energy density stored in the condensate. However, when dealing with inflationary perturbations, one must keep in mind that computations should refer to modes, whose length scales are relevant to the dynamics being studied. For curvature perturbations, the use of a finite box was described in Ref. [76]. For preheating, since thermalization proceeds through particle interactions, the relevant length scales are those that allow for particle interactions, meaning subhorizon scales or short wavelengths.
While UV modes encounter the complication of possibly being excited for wave numbers that exceed the unitarity bound (this does not occur for Higgs modes), the IR modes have a different conceptual difficulty: since thermalization occurs when particles interact and exchange energy, in order to lead to a thermal distribution modes that are superhorizon are "frozen in" and hence cannot take part in such processes. 8 Hence, it is normal to only consider modes that have large enough physical wave numbers that place them inside the horizon at the instant in time that we are considering. Modes that have longer wavelengths are frozen outside the horizon and do not contribute to the thermalization process. They should be summed over and added to the local background energy density. We will skip this last step, as their contribution is subdominant compared to the energy density stored in the inflaton condensate. In Fig. 3 we see the evolution of the comoving Hubble radius, which shrinks during inflation and grows afterwards. We also see that different values of ξ lead to different post-inflationary evolutions, which is expected, since the effective equation of state of the background dynamics after inflation depends strongly on ξ, as shown in Ref. [61]. More specifically, large nonminimal couplings ξ ≳ 100 lead to a prolonged period of matter-dominationlike expansion, which can last for several e-folds in the absence of backreaction. As we will see in the next sections, the majority of the parametric resonance effects occur for N ≲ 3 e-folds, placing the entirety of the reheating dynamics inside the matter-dominated background era for large values of ξ. In order to consistently take into account the relevant wave numbers, we use an adaptive code that only sums up the contributions of modes that are inside the horizon at the point in time when we compute the energy density of the Higgs field fluctuations.

B. Preheating
We now move to the computation of the energy density in the Higgs particles that are produced during preheating. A detailed analysis was performed in Ref. [63]. However, all computations were initialized at the end of inflation, thereby neglecting the amplification of long-wavelength modes during the last e-folds of inflation. We initialize all computations at 4.5 e-folds before the end of inflation, in order to ensure that all relevant modes are well described by We see in the right panels of Fig. 4 that at early times (before the end of inflation) the energy density in Higgs modes (indicated by the solid blue line) decays as a −4 (indicated by the dotted line), in keeping with the expectation for modes in the BD state. However, approximately one e-fold before the end of inflation, the evolution of the energy density in Higgs modes departs from a −4 , because the low-k modes are enhanced with respect to the BD spectrum. This enhancement occurs because m 2 eff;h < 0, an early tachyonic amplification phase driven largely by the effect of coupled metric perturbations. An immediate consequence of this fact is that one would underestimate the true amount of growth by starting the computation in a BD-like vacuum state at the end of inflation. FIG. 4. Left: The effective mass squared (black-dotted curve), along with the contributions from the potential (blue curve) and the coupled metric perturbations (red curve). Right: The energy density in the background Higgs condensate (orange curve) and the Higgs fluctuations (blue curve) for ξ ¼ 10; 10 2 ; 10 3 (top to bottom). The green curve shows 10% of the background energy density, which is used as a proxy for the limit of our linear analysis. The orange-dashed curve is ρ 0 a −4 , corresponding to the redshifting of the background energy density during radiation-dominated expansion. The energy density is plotted in units of M 4 Pl .
The right panels of Fig. 4 present the results for the energy transfer into Higgs particles for ξ ¼ 10; 10 2 ; 10 3 . Preheating is completed when the energy density in the Higgs fluctuations (blue line) becomes equal to the energy density of the background field (orange line). However, the linear analysis is expected to break down when the energy density of the Higgs fluctuations becomes comparable to that of the inflaton field. As an indicator of the validity of the linear theory, which neglects backreaction of the excited modes onto the background, the green line shows 10% of the energy density of the inflation field.
For all values of ξ studied, the system exhibits an amplification of inflaton (Higgs) fluctuations. This is mainly caused by the periodic negative contribution of m 2 3;h to the effective mass squared m 2 eff;h , which is plotted in the left panels of Fig. 4. This is the term arising from considering the effect of the coupled metric perturbations at linear order. As shown in Ref. [63] and further reiterated in Fig. 4, the amplification driven by m 2 3;h lasts longer for larger values of ξ. Specifically, the time at which the tachyonic resonance regime stops scales as t ∼ ffiffi ffi ξ p H −1 end , as shown in Ref. [63]. However, for ξ > 100 the differences are irrelevant (in the simplified linear treatment), since the Universe will have preheated already by N ≃ 3 e-folds. Hence, for ξ > 100 the self-resonance of the Higgs field leads to predictions for the duration of preheating that are almost independent of the exact value of ξ.
After the tachyonic resonance has ended (and if preheating has not completed yet), the modes undergo parametric resonance, driven by the oscillating effective mass term m 2 1;h . However, for very long-wavelength modes k ≃ 0 the Floquet exponent vanishes [62] and the amplification is polynomial in time rather than exponential, and hence significantly weaker. As shown in Ref. [62], the maximum Floquet exponent in the static universe approximation is μ k;max T ≈ 0.3, where T is the background period. Using the relation ω=H ≃ 4, which was derived in Ref. [61] for ω ¼ 2π=T, the maximum Floquet exponent is expressed as μ k ∼ 0.5H. Hence, the Floquet exponent is too small to lead to an efficient amplification of Higgs fluctuations in an expanding universe. Thus the early-time tachyonic resonance driven by the coupled metric fluctuation is crucial for preheating the Universe through Higgs particle production.
For ξ ¼ 10 the situation is significantly different. Both tachyonic resonance (due to the coupled metric fluctuations encoded in m 2 3;h ) as well as parametric resonance (due to the potential term m 2 1;h ) become inefficient earlier, leading to a slower growth of the fluctuations and the energy density that they carry and an incomplete preheating. However, for smaller values of the nonminimal coupling ξ ¼ Oð10Þ one must take into account another important feature, namely, the evolution of the background. As shown in Ref. [61], larger values of ξ put the Universe into a prolonged matter-dominated state (w ¼ 0). This means that the energy density of the background condensate redshifts as a −3 ¼ e −3N . For small values of ξ, however, the Universe passes briefly through the background (average) equation of state w ¼ 0 and after the first e-fold approaches w ≃ 1=3. Figure 5 shows the evolution of the energy density in Higgs modes for the marginal case of ξ ¼ 30. We see that the fluctuation energy density in the Higgs modes would always be smaller than the background if the background evolved with w ≃ 0, as indicated by the orange dashed line. However, the fact that the background energy density redshifts faster (w ≃ 1=3) allows for complete preheating. Simply put, nonminimal couplings in the "intermediate" regime of ξ ¼ Oð10Þ exhibit a shorter period of tachyonicparametric amplification, while at the same time following a background evolution of ρ ϕ ∼ e −4N .
We distinguish two time points relevant for preheating: N reh is the time at which the energy density in the linear fluctuations equals the background energy density (which we take as the time of complete preheating), and N br is the time at which the energy density in the linear fluctuations equals 10% of the background energy density (which is the point at which backreaction effects may become important). We have numerically found that the self-resonance of the Higgs field becomes insufficient to preheat the Universe at ξ < 30. In particular, the results for N reh ðξÞ can be fitted by a simple analytical function, as shown in Fig. 6: for ξ ≳ 30, where complete preheating is possible, at least in the linear approximation that we used. For ξ > 100, N reh becomes largely independent of ξ, as expected from the results of Fig. 4. As a final note, we must say that the results were insensitive to the exact value of the maximum wave number considered. This is due to the fact that the small (but subhorizon) wave numbers k ¼ OðH end Þ are exponentially amplified and dominate the fluctuation energy density shortly after the end of inflation. Hence, we do not need to implement any scheme to subtract the vacuum contribution from large-k modes, since it is vastly subdominant for any reasonable UV cutoff.

A. Evolution during inflation and initial conditions for preheating
We use the equations of motion derived in the Abelian model in the unitary gauge, in order to study the evolution of gauge fields during inflation. The unitary gauge is well defined in this period, since φðtÞ does not vanish. The values of B AE;L k at the end of inflation serve as initial conditions for preheating. Accurate knowledge of the spectrum of gauge fields at the end of inflation is essential, especially when initializing lattice simulations which are increasingly expensive to start deeper within inflation. During preheating, the unitary gauge is not well defined at moments when φðtÞ ¼ 0, so we use the Coulomb gauge. In order to determine the initial condition for θ k , we will use Eq. (42), which relates B L k in the unitary gauge to θ k in the Coulomb gauge.
The equations of motion for the longitudinal and transverse modes in the unitary gauge in conformal time τ are where e is the Uð1Þ gauge coupling. These equations are of the form with I denoting either L or AE polarization, and b I and ω 2 I are given by After integrating by parts, we rewrite the quadratic action in Fourier space as and follow the same quantization procedure as the one appearing in Ref. [87]. This is the standard method used to quantize models with noncanonical kinetic terms, which include nonminimally coupled models in the Einstein frame. The canonical momentum is and the commutator relation of the operatorB I k ðτÞ is We decompose the field operatorB I k ðτÞ in terms of creation and annihilation operators, where the mode function u I k ðτÞ satisfies the same equation of motion as the field operatorB I k ðτÞ [Eq. (100)]. As long as the adiabaticity condition j ∂ τ ω ω 2 j ≪ 1 holds [87], the modes can be described by the WKB approximation, The behavior is different for modes with jkτj > x c (early times/short wavelengths) and jkτj < x c (late times/long wavelengths), with x c given by corresponding to the ratio of the gauge boson mass to the Hubble scale during inflation. Both cases can be described by the WKB approximation, but they exhibit different behavior, which we describe below. At early times and for large values of the wave number k the wave function u I k ðτÞ must match the Bunch-Davies vacuum solution. We focus on the longitudinal mode first. We can take the limit of early times (or subhorizon modes) analytically, when 1 ≪ jkτj ≃ k=ðaHÞ, resulting in and b L ðk; τÞ → M 2 Pl a 2 e 2 φ 2 k 2 2f : ð111Þ Putting everything together, the mode function for jkτj > x c becomes where we used 2f ≃ ξφ 2 and x 2 c ¼ M 2 Pl e 2 =ðξH 2 Þ. The transverse modes are canonically normalized and furthermore conformally coupled at early times, and hence their mode function becomes Overall, α L;AE ¼ 1 and β L;AE ¼ 0 in Eq. (108).

Single-field attractor strength from gauge interactions
The superhorizon evolution (k ≪ aH) of isocurvature fluctuations is an indicator of the (in)stability of the classical background trajectory (see, e.g., Ref. [88]). We will analyze the behavior of the gauge fields and the possible effects on the stability of the single-field attractor.
We will mainly focus on the longitudinal mode, since it will be amplified most efficiently during preheating.
During inflation we can rewrite the equations of motion using x ¼ −kτ as the time variable. If we further make use of the de Sitter approximation (τ ¼ −1=aH) and take φðtÞ as a constant, the equations of motion become 9 with x c as defined in Eq. (109). As expected, we recover the solution of Eq. (112) in the limit jkτj ≫ x c . By using the relation between ξ and λ given in Eq. (48) which is required by the normalization of the power spectrum, Eq. (109) gives where we took e ≃ 1 for Standard Model gauge couplings during inflation. For ξ ≳ 1, where the CMB observables and the inflationary dynamics fall into the "large-ξ" attractor, x c ≫ 1 for all values of interest. Thus, the first of the two cases that were examined in Ref. [87], x c < 1 and x c > 1, does not apply for nonminimally coupled models of inflation with large ξ unless one takes a very weakly coupled gauge field e ≪ 1, making such a value very different to gauge couplings found in the SM. For the longitudinal mode B L the presence of a firstderivative term is important for x < x c , leading to The details of the derivation are given in the Appendix. Figure 7 shows the evolution of certain wave numbers from kjτj ≫ x c ≫ 1 until the end of inflation. It is evident that the simple scalings of ju L ðkjτj ≫ x c Þj ∝ jτj and ju L ðkjτj ≪ x c Þj ∝ ffiffiffiffiffi jτj p agree very well with the full numerical evolution across a wide range of wave numbers. While ξ ¼ 1000 was chosen for Fig. 7, different values of a nonminimal coupling ξ ≫ 1 lead to similar results.
Following Eq. (116), gauge fields during Higgs inflation become very massive, when compared to the Hubble scale. This further reinforces the single-field description of the 9 It is worth noting that the equations of motion for the gauge fields during inflation look very similar in structure to the ones derived for a minimally coupled charged inflaton in Ref. [87]. As shown in Ref. [70], the field space is asymptotically flat for large field values, and hence all covariant derivatives can be substituted for partial derivatives during inflation, at lowest order in 1=ξ. background trajectory discussed in Sec. II C, since the orthogonal direction(s)-described equally well through the scalar d.o.f. θ or through the longitudinal polarization of the gauge boson B L -are very massive, making the background single-field trajectory a stable one. 10 It is interesting to note that, due to the relation between ξ and λ arising from the normalization of the power spectrum, the gauge fields become less massive for larger ξ, meaning that the ratio of the gauge field mass to the Hubble scale becomes smaller. Hence, the single-field attractor, at least in the linearized analysis, becomes weaker for larger ξ. This is opposite to the case of a scalar-only multifield model with a nonsymmetric potential, where the attractor strength increases with ξ, as shown in Ref. [61]. While for the SM the gauge couplings are large enough to make the gauge field much heavier than the Hubble scale, one can construct more general inflationary models involving a Higgs-like field and the associated gauge sector. In this case weakly coupled gauge sectors might leave observational imprints through oscillations of the background during inflation. A search for "primordial clocks" [89,90] in these models is beyond the scope of the present work because they do not arise in SM Higgs inflation, but they could provide a useful tool for exploring gauge-field phenomena in broader classes of nonminimally coupled inflation.
The transverse modes are significantly easier to analyze, since Eq. (114) makes clear that the B AE are conformally coupled at early times and will become massive (and thus be suppressed) for x < x c . In the de Sitter approximation, Eq. (115) can be solved exactly using Hankel functions, resulting in where z ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi 1 4 − x 2 c q , as described, e.g., in Ref. [87]. The analysis of the transverse modes is essentially identical to the minimally coupled case of Ref. [87]. Since they will not be significantly amplified during preheating, we will not discuss them further.

Initial conditions in the Coulomb gauge
Having explored in detail the behavior of the longitudinal gauge fields during inflation, we focus on their form close to the end of inflation and the start of the (p)reheating era. Since we are interested in the details of the vacuum (as will be evident later), we compare the adiabatic vacuum during inflation [given in Eq. (108)] to the approximate analytic expressions derived for u L ðk; τÞ, as well as to the numerically derived values. It is a straightforward exercise to expand b L ðk; τÞ and ω L ðk; τÞ in the two limiting cases of kjτj to see that the WKB expression given in Eq. (108) with α ¼ 1 and β ¼ 0 matches Eqs. (112) and (117) in the appropriate limits. Figure 8 shows the comparison of the WKB solution of Eq. (108), the approximate expressions of Eqs. (112) and (A8), as well as the numerical results from the modes shown in Fig. 7. We see an excellent agreement between all three, with the exception of the modes around kjτj ∼ x c , where the approximate expressions fail, since they were derived using the limits kjτj ≫ x c or kjτj ≪ x c . We must also note that we used the approximation τ ¼ −1=aH for the analytically derived expressions, and hence we expect some discrepancy close to the end of inflation. This agreement has a significant physical meaning: since the adiabatic vacuum follows the evolution of the mode functions, there is no particle production during inflation. We can thus begin our numerical computations at the end of inflation, unlike the case of Higgs self-resonance where we needed to initialize our simulations several e-folds before the end of inflation in order to capture nontrivial dynamics that took place during the last stages of inflation itself.
The initial conditions in the Coulomb gauge can be easily read off from the unitary gauge solutions using Eq. (42). It is interesting to note that there is no ξ-dependent term in the relation of θ k to B L k . The initial conditions that we will use for the computations in the Coulomb gauge are For the brown curve this does not occur during inflation. 10 We must note here, that our analysis only shows the linearized stability of the single-field trajectory, not the approach towards it from generic initial conditions fΦ; ∂ t Φg. The latter was performed in Ref. [71] for an SOð4Þ-symmetric model meant to describe Higgs inflation without gauge couplings.
Before we conclude the analysis of the gauge-field evolution during inflation, let us focus on the case of kjτj ≫ x c , where the initial conditions for preheating are It is reassuring that for large wave numbers the coupling constant e drops out of the initial conditions for the θ field (since x c ∝ e), and hence the decoupling limit is trivially obtained. For kjτj < x c it is slightly more complicated to see this, since for e → 0 we get x c → 0, and hence that region shrinks into nonexistence as we take the decoupling limit. Also, we would have to compute the expressions for x c ≪ 1 before we send e → 0 in that case. Since the case of e ≪ 1 does not apply to Higgs inflation, we will not pursue it further.

B. Preheating
We start by rewriting Eq. (39) in a somewhat more compact way, where we defined the gauge-field mass and X θ ¼ aðtÞ · θ. We normalize the scale factor as a ≡ 1 at the end of inflation. The effective mass of the Goldstone mode θ in the absence of gauge fields is The numerical solution of Eq. (123) was obtained in cosmic rather than conformal time, since this is more convenient for numerical simulations after the end of inflation. The computations were initialized at the end of inflation, according to Eqs. (119) and (120).
We can follow the quantization method described in Ref. [61] and utilized in Sec. IV for the study of the Higgs self-resonance, where e θ 2 ¼ ffiffiffiffiffiffi ffi G θθ p . Using the vielbein decomposition, the covariant derivatives are effectively substituted by partial ones, In order to eliminate the first-derivative term we can use the rescaled variablez k , defined as leading to where wherem 2 B is larger than m 2 1;θ and m 2 4;θ . As discussed extensively in Refs. [61][62][63] for the case of a purely scalar multifield model with large nonminimal couplings to gravity, the field-space manifold is asymptotically flat for large field values and exhibits a curvature "spike" at the origin φðtÞ ≃ 0. This "Riemann spike" is exhibited in the effective mass of the isocurvature modes m 2 eff;θ , more specifically in the m 2 2;θ component, which is subdominant for all times away from the zero crossings of the background value of the inflaton field φðtÞ. We will not reproduce the entirety of the Floquet structure of this model, both because we do not wish to repeat the analysis of Ref. [62], and because (as we will see in the subsequent section) the first zero crossing of φðtÞ is the only relevant one for preheating through gauge modes.
In order to estimate the maximum excited wave number k max , we consider the following approximation containing only the dominant terms: wherem 2 B dominates over all subsequent terms in Eq. (134) for large k. Figure 9 shows the three contributions to ω 2 z;approx for ξ ¼ 10 3 ; 10 4 . As shown in Ref. [63], the scaling of the spike in the effective mass is where hHðtÞi is a time-averaged version of the Hubble scale over the early oscillatory behavior. The range of excited wave numbers is given by the relation assuming that the spike of m 2 2;θ dominates overm 2 B near φðtÞ ¼ 0. Each subsequent inflaton zero crossing affects a smaller range of wave numbers, since m 2 2;θ ∝ H 2 ∝ ρ infl: ∝ a −3 , where we assumed w ¼ 0 for the averaged background evolution. Altogether k 2 max ∝ a −1 , and hence the maximum excited wave number shrinks for every subsequent inflaton oscillation. The maximum comoving wave number after the first inflaton zero crossing, where aðtÞ ≈ 1, is where we used Eq. (47) and H end ≈ 0.5H infl . This is in agreement with Ref. [64]. We focus primarily on the first inflaton zero crossing, since the produced gauge bosons will decay into fermions between two subsequent background zero crossings, and hence Bose enhancement is lost. This was shown in Refs. [58,59] and will be discussed in detail in Sec. VI. The second dominant component of the gauge-field effective frequency squared ism 2 B , which scales simply as where the λ − ξ relation given in Eq. (48) was used in the last step. We can see that for ξ ¼ 10 3 the maxima of the two contributionsm 2 B and m 3 2;θ are similar, as shown in Fig. 9. Computing the energy density transferred from the inflaton condensate into the gauge-field modes requires more attention than the corresponding computation of Sec. IV for the Higgs self-resonance. In the case of Higgs self-resonance, the range of excited wave numbers is k h max ∼ H. A naive computation of the energy density in the local adiabatic (WKB) vacuum for the same modes gives ρ BD ∼ k 4 max ∼ H 4 , which is 10 orders of magnitude smaller than the background energy density. 11 In that case, we do not need to subtract this unphysical vacuum contribution from the energy density of the Higgs modes, since the energy density in the parametrically amplified modes is exponentially larger. 11 Any computation that does not involve vacuum subtraction, including lattice simulations (such as in Refs. [53,54]), deals with classical quantities and computes the energy density of the vacuum modes as if they were physical. Such a computation is valid as long as the unphysical energy density of the vacuum modes is vastly subdominant.
For the case of gauge fields the maximum wave number up to which modes can be excited is given in Eq. (138). The vacuum energy density in these modes, naively computed, is ρ BD ∼ k 4 max ∼ λ 2 M 4 Pl . The total energy density in the inflaton field is ρ infl ¼ 3H 2 M 2 Pl , leading to ρ BD =ρ infl ∼ λξ 2 ∼ 10 −10 ξ 4 . This is much greater than unity for large values of the nonminimal coupling. We thus need to remove the unphysical vacuum contribution to the energy density using the adiabatic subtraction scheme [9]. In this scheme we compare the wave function of the gauge fields to the instantaneous adiabatic vacuum, computed in the WKB approximation, isolating the particle number for each wave number k. The particle number corresponding to a mode v k is given by A drawback of this method is that the particle number is only well defined when the adiabaticity condition holds, _ ω k =ω 2 k ≪ 1, and thus we cannot define the particle number in the vicinity of the "Riemann spike" when φðtÞ ¼ 0. 12 The energy density is easily computed through the particle number as Both the particle number and the energy density can be computed equally well using the field θ k or B L k , since the only moment for which the longitudinal gauge fields are not defined is when φðtÞ ¼ 0. At this instant we cannot define the particle number either way, since there is no well-defined adiabatic vacuum. Figure 10 shows the evolution of the particle number density for a few values of the comoving wave number after the first few inflaton zero crossings, neglecting the effect of particle decays, as described in Sec. VI. The left panel of Fig. 11 shows the particle number density per k mode for ξ ¼ 10; 10 2 ; 10 3 after the first inflaton zero crossing. The condition of Eq. (138) for the maximum excited wave number k max is evident.
At this point, it is worth performing a simple estimate of the energy density that can be transferred to the gauge-field modes away from the first point φðtÞ ¼ 0, where hni is the average occupation number. The background inflaton energy density is ρ infl ¼ 3H 2 M 2 Pl ∼ 10 −11 M 2 Pl , and hence for ξ ≳ 10 3 the transfer of energy is enough to completely drain the inflaton condensate within one zero crossing of φðtÞ if we take the particle number shown in Fig. 11 into account. The right panel of Fig. 11 shows the ratio of the energy density in gauge fields to the background energy density of the inflaton after the first zero crossing. Obviously, values of ρ gauge =ρ infl > 1 are not physical but signal the possibility of complete preheating. FIG. 9. Dominant components of the effective frequency squared for ξ ¼ 10 3 (left) and ξ ¼ 10 4 (right). Color coding is as follows: m 2 B =a 2 (red), m 2 2;θ (blue), and k 2 =a 2 (black) for the maximum excited wave number k max . The orange-dotted curve shows the scaling a −3 . 12 Reference [64] computed the particle number, working in the Jordan frame, and arrived at similar results. The energy of the gauge fields was subsequently computed using the value of the gauge-field mass directly on the "Riemann spike." We refrain from using m 2 2;θ j max as an indicator of the gauge-field mass, since the particle number is not a well-defined quantity there. For ξ ≈ 10 3 , the two contributions to the gauge-field mass m 2 2;θ andm 2 B are comparable, as shown in Fig. 9, which does not hold for other values of ξ.

C. Unitarity scale cutoff
So far we have computed the excitation of gauge-field modes of arbitrary wave number k < M Pl . However, the unitarity scale sets a limit above which no analytical (perturbative) treatment can be trusted. The unitarity scale for Higgs inflation, and more generally for nonminimally coupled models, has received extensive attention in the literature. We will follow the analysis of Ref. [67], where a field-dependent unitarity scale was derived in both the Jordan and Einstein frames.
The unitarity scale at the end of inflation is k UV;1 ≡ M Pl = ffiffi ffi ξ p , which becomes k UV;2 ≡ M Pl =ξ for even smaller values of the background Higgs field. It is straightforward to estimate the relation of the unitarity scale to the maximum excited wave number, We see that, depending on the value of the nonminimal coupling ξ, the wave number of the produced gauge bosons can exceed the field-dependent unitarity scale. New physics is needed above the unitarity scale and it is not clear how this new physics will change particle production for such large wave numbers. We do not wish to propose any UV completion of the Standard Model in order to address the dynamics above the unitarity scale. We will instead provide a conservative estimate of the energy density in gauge bosons in the presence of unknown UV physics that suppresses particle production with large wave numbers (above the unitarity scale). Simply put, we will compute the energy density by introducing a UV cutoff at k UV;1 or k UV;2 . If we consider the UV cutoff at k UV;1 , both ξ ¼ 10 3 and ξ ¼ 10 4 preheat entirely after one inflaton zero crossing, FIG. 10. The particle number density for k=H end ¼ 1; 150; 550; 2600; 28000 (blue, black, green, red, and purple, respectively). From left to right: ξ ¼ 10 2 ; 10 3 ; 10 4 . If a colored curve is missing from a panel, the corresponding wave number is not excited.
FIG. 11. Left: The particle number density after the first inflaton zero crossing for ξ ¼ 10 2 ; 10 3 ; 10 4 (blue, orange, and green, respectively) Right: The ratio of the energy density in gauge fields to the background inflaton energy density as a function of the nonminimal coupling ξ after the first zero crossing. We see that for ξ ≳ 10 3 gauge boson production can preheat the Universe after one background inflaton zero crossing, and hence it is much more efficient than Higgs self-resonance. since k UV;1 ≳ k max ðξ ¼ 1000Þ, as can be seen from Fig. 11. If instead we place the UV cutoff at k UV;2 , the gauge fields do not carry enough energy to completely preheat the Universe after one inflaton zero crossing, regardless of the value of the nonminimal coupling ξ. We thus conclude that preheating into gauge fields is very sensitive to unknown UV physics, since the majority of the energy density is carried by high-k modes, whose number density in a UVcomplete model can be much different than the one computed here. It is worth noting that the excitation of Higgs fluctuations occurs entirely below the unitarity scale, and hence it is not UV sensitive. We will not consider any UV cutoff for the remainder of this work, unless explicitly stated.

VI. SCATTERING, DECAY, AND BACKREACTION
So far we have computed the parametric excitation of particles, either Higgs or gauge bosons, from the oscillating Higgs condensate during preheating. With the exception of the brief discussion in Sec. IVA, the interactions of the resulting particles have been completely ignored. However, as discussed in Refs. [58,59], certain types of decays of the produced particles can suppress Bose enhancement and thus effectively shut off preheating. We will discuss in turn (A) the decay of Higgs particles into gauge bosons and fermions, (B) the scattering of Higgs particles into gauge bosons and fermions, (C) the decay of parametrically produced gauge bosons, (D) the scattering of gauge bosons into fermions and Higgs bosons, and (E) possible effects arising from non-Abelian interactions of the produced W and Z bosons. Any of the above-mentioned processes can suppress or shut off the resonances. Due to their inherent differences, we will explore them separately.

A. Higgs decay
In the Standard Model, Higgs particles can decay into pairs of fermions or gauge bosons. The fermion masses are while the gauge boson masses were extensively studied in Sec. V B. For now, it is enough to consider the part of the gauge-field mass analogous to m f in Eq. (145) with the Yukawa coupling substituted by the gauge coupling. We start with the process of a Higgs particle decaying into two gauge bosons. In order for this to be kinematically allowed, the following relation must hold: m h > 2m gauge . It is straightforward to see that m h ≪ 2m gauge , at least for φðtÞ ≠ 0. When φðtÞ ¼ 0, the Riemann contribution to the gauge-field mass (the "Riemann spike") dominates, keeping the relation m h ≪ 2m gauge valid at all times. Hence, the Higgs field cannot decay into gauge bosons, as long as the background Higgs condensate follows the evolution that is derived neglecting backreaction.
The decays of Higgs bosons to fermions deserve closer attention, due to the fact that small Yukawa couplings for some fermions (like electrons and positrons) can make them much lighter than the Higgs particles, and hence kinematically open the decay channel. Furthermore, fermion masses do not have a Riemann component, and hence when φðtÞ crosses zero fermions become instantaneously massless, making the decay even easier. A similar analysis of kinematical blocking of perturbative decays during reheating was performed in Ref. [91], where the Higgs field was a light spectator field during inflation rather than playing the role of the inflaton itself.
We will compute each component of the Higgs field m h;1 and m h;3 separately. We begin with the potential contribution where δ ¼ ξφ 2 and δ ≃ 0.8 at the end of inflation, as discussed in Refs. [61][62][63]. The value in Eq. (146) holds at the start of preheating and until the crossover time t cross ∼ ffiffi ffi ξ p H −1 end . The expression for t cross was derived in Ref. [63]. For t < t cross metric perturbations dominate the effective mass, resulting in tachyonic amplification. For t > t cross the Higgs particle mass m 2 h;1 decreases slowly with time. For ξ ¼ 10; 10 2 ; 10 3 the crossover time occurs at N cross ≃ 1.5; 2.5; 3.2, respectively. Figure 12 shows the ratio of the fermion to the Higgs mass at the end of inflation as a function of the Yukawa coupling, for different values of the nonminimal coupling. We see that the decay is kinematically possible for small Yukawa couplings. Furthermore, the decay channel is less constrained for later times and larger nonminimal coupling. The perturbative decay rate of Higgs particles to fermions is given by Figure 12 shows the ratio Γ=H, which must be greater than unity in order for the decay to be efficient. It is clear that, in the parameter range where the decay is kinematically allowed, it is very inefficient. This can be intuitively understood since m h ∼ H, m f is proportional to y f and Γ is proportional to y 2 f , and hence Γ=H is suppressed by an extra factor of the Yukawa coupling compared to m f =m h . This conclusion does not change, even if one considers the short increase in the mass of the Higgs modes due to the coupled metric fluctuations term m 2 3;h . Even though m 2 3;h has a large positive spike, its duration is too small to allow for a significant decay of Higgs particles into fermions. Before we conclude this section, we will make one further note regarding the evolution of fermion masses. Equation (145) shows that fermions become massless when φðtÞ ¼ 0. The distinction between computing the fermion mass during reheating by either using an averaged quantity for the Higgs VEVor by using the full time dependence was explored in Ref. [91]. In order to explore possible effects of the time dependence of the fermion mass, we focus on the case of ξ ¼ 10 3 and choose a large Yukawa coupling y w ¼ 1, since that provides the largest decay rate to Hubble scale ratio Γ=H ≃ 10, as shown in Fig. 12. The time per oscillation that m f < m h is Δt Ã H ≃ 10 −3 . Hence, Γ=Δt ≪ 1, meaning that the time when fermions are massless is too small to significantly deplete the Higgs boson population.

B. Higgs scattering
While we saw that Higgs decays to both gauge bosons and fermions are either kinematically blocked or extremely weak during preheating, the same might not be true for Higgs scatterings, due to the large occupation number, close to the time of complete preheating. The kinematical blocking arguments still apply, since the relation m h > 2m f;A is replaced by m h > m f;A , and hence is weakened only by a factor of 2. As we saw, the kinematical constraints are significant, and hence we will only consider scattering of Higgs particles into the lightest fermions (electronpositron pairs), with y e ¼ Oð10 −6 Þ. The relevant rate is The Higgs particles are heavy (m h > H) and have small wave numbers (k=a ≲ H), and hence will be nonrelativistic. We will take v ¼ c ≡ 1 as an upper limit. The number density of Higgs particles is approximately where ρ h ¼ ρ infl at the point of complete preheating. The cross section is Putting everything together, we arrive at It is easy to see that Γ=H ≪ 1 since y 4 e ≃ 10 −24 , M 2 Pl =H 2 ≃ 10 10 , and H=m h < 1.
It is also worth briefly noting other scattering diagrams leading to the depletion of the Higgs population. Two examples are shown in Fig. 13, which are the inverse of gluon-fusion processes considered for the LHC. In general, they suffer from the same suppression factors as the treelevel scattering: light fermions come with small Yukawa couplings, while heavy ones will lead to suppression factors from the fermion loops. We will not discuss these processes further.

C. Gauge decay
Following Refs. [58,59], the decay widths of the W and Z bosons to fermions are given by where the decay widths are obtained by summing over all allowed decay channels into SM fermions. The decay of the Z boson to a pair of Higgs particles proceeds similarly. Using the gauge boson mass given in Eq. (139), we see that Γ W;Z =H ≫ 1, and hence the produced gauge-boson population is depleted within a Hubble time, or between two consecutive inflaton background zero crossings φðtÞ ¼ 0.
There are two issues that need to be addressed: the possible decay of particles during their production close to the Riemann spike at φðtÞ ¼ 0 and the decay away from φðtÞ ¼ 0, when them 2 B ¼ e 2 φ 2 ðM 2 Pl =2fÞ component dominates the gauge-field mass. In both cases, we will approximate the total decay of the particle number as where ΓðtÞ is defined through Eqs. (152) and (153) by considering the time-dependent mass of the gauge bosons.
We will focus only on the cases of ξ ¼ 10 3 and ξ ¼ 10 4 . During the spike, the particle number is not a well-defined quantity, since an adiabatic vacuum cannot be constructed due to the violation of the adiabaticity condition. We will however compute the exponential decay factor of Eq. (154) as an estimate of possible particle decays. We choose the limits of integration to correspond to the times for which adiabaticity is violated, and hence particle production occurs. This is also the time at which the Riemann spike is pronounced. For all cases we get e − R t t 0 Γðt 0 Þdt 0 > 0.5, and hence there is no significant particle decay. We will thus neglect this altogether. However, after the particle production has taken place at φ ¼ 0, the particle number is a constant, if one neglects decays, and the particle mass grows sharply as m 2 W;Z ∼ φ 2 . We rewrite the equation for the energy density in the gauge sector as Γðt 0 Þdt 0 ω k : ð155Þ Figure 14 shows the energy density per particle number of a random excited k mode as ρ ≃ nðtÞm A , with A denoting any gauge field. We see that decays into fermions completely deplete the produced gauge-boson population within far less than a period of background oscillations. Hence, in order for the energy transfer to be able to preheat the Universe, the energy density in the gauge fields must be equal to the energy density in the inflaton condensate as soon as the adiabaticity condition is restored. The fact that the particle decays during the "Riemann spike" are insufficient to suppress gauge boson production shows that this is indeed possible.

D. Gauge scattering
Instead of decaying into fermions, gauge bosons can also scatter into Higgs particles or fermion-antifermion pairs. We will estimate the rate of the Higgs scattering to Higgs bosons. The scattering rate is Γ ¼ nσv, where we take v ¼ c and FIG. 14. Energy density per mode (in arbitrary units) with (blue) and without (red) considering particle decays for ξ ¼ 10 3 (left) and ξ ¼ 10 4 (right). The time is rescaled by the period of background oscillations.
where we computed the number density using the condition of complete preheating and we took the Mandelstam variable s ≃ k 2 max . Altogether, where we used the relation between λ and ξ given in Eq. (48). Since Γ=H ≲ 1 for ξ ≳ 10 3 , gauge field scatterings are not important. This is different from other cases of preheating into gauge bosons, such as that in Ref. [54] where gauge boson scattering is extremely efficient. The difference is that in the present case the number density is not large, but the average energy carried by each gauge boson is, due to the large range of excited wave numbers.

E. Non-Abelian effects
Since we are using an Abelian Uð1Þ gauge field as a proxy for preheating into SM W and Z bosons, we must estimate the possible non-Abelian effects. As long as the linear analysis holds, the electroweak sector can be decomposed into three almost identical Abelian copies. A numerical example of the relation between an SUð2Þ gauge field and its three Abelian copies at low field values was shown in Ref. [92].
However, once the gauge-field modes become sufficiently populated, their true non-Abelian nature cannot be neglected. The relevant term in the non-Abelian Lagrangian is where f abc and f ade are SUð2Þ structure constants. In the equation of motion for the gauge field strength A i , this term in the Lagrangian will induce a term of the form g 2 A j A j A i , which has the form of an effective non-Abelian mass term. Using a Hartree-type approximation we can define the non-Abelian contribution to the gauge-field mass squared as m 2 non-Abelian ∼ g 2 hAAi. We estimate hAAi through the energy density of the gauge fields as ρ ≃ m 2 A hA 2 i. Taking as a maximum value ρ In order for the non-Abelian mass contribution to suppress particle production, it must dominate over m 2 θ;2 . However, we know that m 2 θ;2 ≃ ξ 2 H 2 ≃ ξ 2 10 −12 M 2 Pl , meaning that for ξ ≳ 10 3 the "Riemann spike" dominates over the possible non-Abelian mass contribution. Hence, we expect the explosive transfer of energy from the inflaton to the gauge fields to persist even in the full SUð2Þ × Uð1Þ sector.
A further phenomenon that has been observed during simulations of preheating of a non-Abelian Higgsed sector was described in Ref. [93]. There, the decay of the Higgs condensate through resonant decay of electroweak bosons was simulated. Non-Abelian gauge boson interactions led to an extended momentum distribution. Particles with such high momenta are energetic enough to scatter off the Higgs condensate and fragment it, thereby shutting off any further parametric resonance. In the case of Higgs inflation the gauge fields produced do not survive long before decaying into fermions, due to their large masses. Hence, this is unlikely to be an issue in the present case.

VII. OBSERVATIONAL CONSEQUENCES
Observing reheating is difficult due to the inherently small length scales involved. However, there are two important quantities that can be used to connect reheating to particle physics processes or CMB observables: the reheat temperature T reh , and the number of e-folds of an early matter-dominated epoch in the expansion history of the Universe N matter .
The reheat temperature is computed using the Hubble scale at the instant when ρ infl ¼ ρ rad as where g Ã ¼ 106.75 is the number of relativistic d.o.f. at high energies. For instantaneous reheating from gauge field production, which happens for ξ ≳ 1000, the Hubble scale is H ≃ H end . For ξ ≲ 1000 preheating proceeds through Higgs self-resonance, leading to a smaller value of the energy density, as shown in Fig. 4. The monotonic increase of the reheat temperature T reh as a function of the nonminimal coupling ξ is shown in Fig. 15. It must be noted that Eq. (161) assumes the immediate transition to a thermal state after preheating has ended. For the case of Higgs selfresonance, this will occur through efficient scattering of Higgs bosons to the rest of the SM. For the case of instantaneous preheating to gauge fields, the situation is more complicated. In that case the number density of gauge bosons is not exponentially large, as is usually the case in preheating. On the contrary, the transfer of energy to gauge fields is done primarily through the production of fewer high-momentum modes k max ∼ ffiffi ffi λ p M Pl . A fraction of the produced W and Z bosons will decay to leptons, while another fraction will decay into quarks and antiquarks that will eventually hadronize. The approach to thermal equilibrium will thus be more complicated. We leave the study of the thermalization process for future work and we use Eq. (161) as an estimate of the reheat temperature, under the assumption of efficient thermalization.
However, a high reheat temperature may pose a challenge for any computation that goes beyond the linearized analysis that we presented, due to possible conflicts with the unitarity scale. Since thermalization of the reheating products will result in a blackbody spectrum, we can take the typical momentum involved to be k ∼ 3T reh , which is thus the typical momentum exchange in particle scatterings inside the plasma. Since complete reheating means that the inflaton condensate will have completely decayed, the unitarity scale is k UV;2 ≡ M Pl =ξ. The typical particle momenta are below the unitarity scale for 3T < k UV;2 . As shown in Fig. 15, for ξ ≲ 300, the resulting plasma has a low enough temperature to avoid processes that exceed the unitarity scale, at least neglecting the tail of the thermal spectrum. For ξ ≳ 300, the unitarity scale k UV;2 will be exceeded by the typical wave numbers in the system, provided that thermalization is efficient. Finding out whether the unitarity scale is indeed violated for ξ ≳ 300 requires a more extensive study of the thermalization phase, which lies outside the scope of the present work.
The number of matter-dominated e-folds of postinflationary expansion is a nonmonotonic function of the nonminimal coupling. For ξ ≳ 10 3 , instantaneous reheating leads to a universe filled with gauge-field modes of high wave numbers, and hence the Universe transitions immediately to radiation domination (assuming no UV suppression). We must note that the decay of the inflaton condensate makes the gauge fields light, and hence relativistic. For small values of the nonminimal coupling ξ ¼ Oð10Þ the background evolves as w ≈ 1=3, and hence the evolution of the Universe is that of radiation domination soon after the end of inflaton, even if preheating is not efficient. Hence, N matter ¼ 0 for both large and Oð10Þ values of the nonminimal coupling. There is an intermediate region of ξ ¼ Oð100Þ, where preheating happens through self-resonance and the background evolves following an average equation of state of w ≈ 0 [61] before preheating completes. In that regime of nonminimal couplings N matter ≈ N reh ≈ 3, slightly shifting the predictions of the CMB compared to the approximation of instantaneous reheating [47], where the equation of state is assumed to transition from w ¼ −1=3 at the end of inflation to w ¼ 1=3 immediately afterwards.

VIII. CONCLUSIONS
Higgs inflation is an appealing way to realize inflation within the particle content of the Standard Model, by coupling the Higgs field nonminimally to the gravity sector with a large value of the nonminimal coupling. We analyzed the nonperturbative decay of the Higgs condensate into Higgs bosons and electroweak gauge fields, finding distinct behavior for different ranges of values of the nonminimal coupling ξ.
The self-resonance of the Higgs field leads to preheating after N reh ≃ 4 e-folds for values of the nonminimal coupling ξ ≳ 30. For large values ξ > 100 the inflaton can transfer all of its energy into nonrelativistic Higgs modes within N reh ≈ 3, independent of the exact value of the nonminimal coupling. The dominant contribution to the parametric excitation of Higgs modes is the effect of coupled metric fluctuations. In order to accurately capture the amplitude of the Higgs wave function, the computation must be initiated before the end of inflation.
The excitation of gauge bosons is much more dramatic, reminiscent of the purely scalar case of preheating in multifield inflation with nonminimal couplings [61][62][63]. Gauge fields are excited after the first zero crossing of the inflaton field, up to wave numbers k max ∼ ffiffi ffi λ p M Pl . This leads to the possibility of the inflaton condensate transferring the entirety of its energy density to W and Z bosons immediately after the end of inflation, leading to instantaneous preheating. W and Z bosons will efficiently decay into SM fermions, ultimately filling the Universe with a thermal plasma. Estimates of perturbative decay and non-Abelian effects show that gauge field production is robust against both.
The efficiency of the reheating stage can have observational consequences. The values of the spectral observables n s and r are related to the time N Ã when the CMB-relevant modes exited the horizon during inflation. For Higgs inflation and related models the CMB observables are given by n s ≃ 1-2=N Ã − 3=N 2 Ã and r ≃ 12=N 2 Ã . Depending on the speed of the transition from the end of inflation to radiationdominated expansion of the Universe, the observationally relevant N Ã may vary, shifting the predictions for n s and r.
The use of the Coulomb rather than unitary gauge for our computations allowed us to tie the results to the purely scalar case studied in Refs. [61][62][63], as well as apply the results to other models with curved field-space manifolds. One such example is another version of Higgs and Higgs-like inflation, proposed in Ref. [94]. In that model, the necessary nonminimal coupling is small and negative, accompanied by a minimum of the Higgs potential at a large vacuum expectation value during inflation. The analysis of this model is left for future work and can provide a possible method for probing the Higgs potential during inflation through its effect on the preheating behavior and the reheat temperature, rather than the CMB observables alone.
Another modification of Higgs inflation is based on the possible existence of an inflection point in the Higgs inflation potential [95][96][97]. The possibility of primordial FIG. 15. Reheat temperature in units of M Pl as a function of the nonminimal coupling ξ. The discontinuity at ξ ≃ 10 3 occurs due to the instantaneous preheating to gauge fields. The light red region represents the uncertainty of the exact threshold of instantaneous preheating to gauge fields. The black-dotted line corresponds to the unitarity scale constraint T reh ≲ M Pl =3ξ. The blue-dashed line shows the reheat temperature due entirely to Higgs self-resonance, assuming gauge boson production above the unitarity scale is suppressed due to unknown UV physics. black hole production in critical Higgs inflation has been proposed in Ref. [98]. It was however shown that large curvature fluctuations in Higgs inflation are hard or impossible to produce without violating observational limits on the tensor-to-scalar ratio and the running of the spectral tilt [99,100], unless one postulates a large running of the nonminimal coupling ξ that is not found in the renormalization group flow of the Standard Model nonminimally coupled to gravity. Recent studies of Higgs inflation involving nonminimal couplings in the Palatini formulation of gravity [101,102] can also have different preheating phenomenology. Exploring the preheating phenomenology of these models is interesting and can be performed using the techniques applied here. Such analyses can provide unique methods to probe the Higgs potential at energy scales that are out of reach for the LHC and any future accelerator.
where we dropped an arbitrary pure phase term. The derivative is ∂ τ u L ðk; τÞ u L ðk; τÞ As a side note, the fact that the term proportional to i is negative shows that we correctly chose the right-moving wave. 13 Again, dropping an arbitrary phase, the initial conditions for preheating computations are for wave numbers such that jkτ in j < x c . Since x c ≫ 1, we can drop the 1=2 factor in the above equation. 13 In reality, when solving the full equation of motion in cosmic time with all factors included, the result is not a perfect rightmoving wave. However, this is still a very good approximation to use as an initial condition both for the current linear computation of fluctuations as well as for future lattice simulations.