The Decay of the Standard Model Higgs after Inflation

We study the nonperturbative dynamics of the Standard Model (SM) after inflation, in the regime where the SM is decoupled from (or weakly coupled to) the inflationary sector. We use classical lattice simulations in an expanding box in (3+1) dimensions, modeling the SM gauge interactions with both global and Abelian-Higgs analogue scenarios. We consider different post-inflationary expansion rates. During inflation, the Higgs forms a condensate, which starts oscillating soon after inflation ends. Via nonperturbative effects, the oscillations lead to a fast decay of the Higgs into the SM species, transferring most of the energy into $Z$ and $W^{\pm}$ bosons. All species are initially excited far away from equilibrium, but their interactions lead them into a stationary stage, with exact equipartition among the different energy components. From there on the system eventually reaches equilibrium. We have characterized in detail, in the different expansion histories considered, the evolution of the Higgs and of its dominant decay products, until equipartition is established. We provide a useful mapping between simulations with different parameters, from where we derive a master formula for the Higgs decay time, as a function of the coupling constants, Higgs initial amplitude and postinflationary expansion rate.


I. INTRODUCTION
Inflation, an early period of accelerated expansion, is the leading framework to explain the initial conditions of the Universe. The concrete particle physics realization of inflation has eluded any clear identification so far, so the inflationary dynamics is often described in terms of a scalar field, the inflaton, with a vacuum-like energy density. Furthermore, the confirmed discovery of the Standard Model (SM) Higgs in the Large Hadron Collider (LHC) [1,2] has initiated the quest for understanding the cosmological implications of the Higgs, and in particular, its possible role during and after inflation. Intriguingly, the SM Higgs could play the role of the inflaton, if a largish non-minimal coupling to gravity is introduced, appropriately fixed to fit the observed amplitude of the cosmological perturbations [3]. This model, known as Higgs-Inflation, constitutes undoubtedly one of the most attractive and economical scenarios for realizing inflation, though is not free of criticism [4,5]; see however [6].
In this paper, we will rather explore a different route for the possible role of the Higgs during and after inflation. We will merely assume that inflation was driven by a very slowly evolving energy density, without specifying the nature of the field responsible for it. Inflation can then be seen effectively as a quasi-de Sitter background with a slowly changing Hubble rate, departing significantly from de Sitter only at the very end of Inflation. We will assume that the SM Higgs is not coupled directly to the inflationary sector. Under these circumstances, the Higgs behaves simply as a spectator field living in a (quasi-)de Sitter background, with the effective potential of the Higgs ultimately dictating its behaviour. Let us note that even if there is no direct coupling, it is likely that effective operators will connect the Higgs with the inflaton, via some possible mediator field(s). Moreover, the need to reheat the Universe after inflation requires somehow the presence of such coupling, though there is no particular constraints on this. As we will see, the Higgs decays very fast after inflation into all SM species, so one can safely assume that the effect of an inflaton-Higgs effective coupling is negligible, unless that coupling is significantly large. Therefore, even if such coupling is present, we will consider it weak enough, so that any interaction Higgs-inflaton does not affect the dynamics of the latter during or after inflation.
The improved renormalized Higgs potential has been computed at next-to-next-to-leading order [7,8]. It is characterised by the running of the Higgs self-coupling λ(µ), which decreases with energy dλ/dµ < 0, and becomes negative above certain critical scale µ 0 , λ(µ ≥ µ 0 ) ≤ 0. Equivalently, the effective potential develops a barrier at large field amplitudes, reaching a maximum height at some scale µ + < µ 0 , so that at higher energies µ > µ + the effective potential goes down, crosses zero at µ = µ 0 and becomes rapidly negative, possibly reaching eventually a (negative) minimum at some scale µ − µ 0 . This can be seen in Fig. 1. These scales depend sensitively on the Higgs mass m H , the strong coupling constant α s , and specially on the top Yukawa coupling y t . For the SM central values, α s = 0.1184, m H = 125.5 GeV, and the most recent measurement of the top quark mass by CMS, m t = 172.38 GeV [9], one finds µ + 2 × 10 11 GeV and µ 0 3 × 10 11 GeV. If one takes the world average top quark mass m t = 173. 36 GeV [10], then µ + , µ 0 are reduced by a factor ∼ 30. However, by considering a value of y t merely 2-3 sigma smaller than its central one, we can push the critical scales to µ + , µ 0 5 × 10 16 GeV. Besides, minimal additions to the SM such as scalar singlet coupled to the Higgs [11], or even a small non-minimal coupling of the Higgs to gravity [12], can also modify the running of λ(µ) and stabilize the effective potential. In such a case, the Higgs self-coupling may remain always positive λ(µ) > 0.
We will consider that the Higgs amplitude during in- flation remains always in the 'safe' side of the effective potential, where λ(µ) is positive. This can be guaranteed if µ + is sufficiently large (compared to the inflationary scale), or alternatively, if beyond the SM physics stabilizes the potential at high energies. With this considerations, the Higgs fluctuates during inflation, like any light degree of freedom. The fluctuations then pile up at super-Hubble scales, creating a condensate. The amplitude of the Higgs condensate, however, will not grow unbounded with the numbers of e-folds, as it happens in the case of a massless free field. On the contrary, the Higgs self-interactions provide an effective (sub-Hubble) mass to the fluctuations, which eventually saturates the growth of the condensate amplitude. The distribution of the Higgs amplitude at super-Hubble scales enters very fast, within a few-efolds, into a self-similar regime, which continues until the end of Inflation. The Higgs condensate acquires this way a fixed physical correlation length (exponentially larger than the Hubble radius) and a fairly large amplitude. This will set up the initial condition for the behaviour of the Higgs after the end of inflation. Notice however that our analysis will not really depend on the condition λ(µ) > 0 during inflation. The possible implication of the Higgs self-coupling becoming negative λ(µ) < 0 during inflation has indeed been analysed in quite some detail [13][14][15][16]. In this case, if the scale of inflation is sufficiently high, the second minimum can be reached and anti-de Sitter bubbles are formed during inflation. The consequence of this does not need to be catastrophic, but rather indicative that either the condition that our Universe is in the electroweak (EW) vacuum is something very special (very improbable), or either some new physics beyond the SM is necessary to stabilize the Higgs potential. The crucial ingredient for our analysis is, therefore, not that the Higgs self-coupling λ(µ) remains positive during inflation, but the fact that the Higgs develops a vacuum expectation value (VEV) during inflation much larger than the electroweak (EW) scale ∼ O(10 2 ) GeV. The way such condensate is attained is mostly irrelevant. However, for the case we consider in this paper, λ(µ) > 0 all through inflation, the formation of a large VEV during inflation is unavoidable, and its typical amplitude is easily predicted.
In this paper, we investigate in detail the Higgs dynamics during the immediate stages following the end of inflation. The paper is divided as follows. In Section II we present a brief analysis of the behaviour of the Higgs after inflation, ignoring its coupling to the rest of the SM species. In Section III, we switch on the coupling to the SM fields, but ignore the gauge nature of the interactions of the SM. We then develop analytical estimates for a later comparison with numerical simulations. In Section IV we present the first set of lattice simulations, where we follow the Higgs and its decay products, again under the assumption that the gauge nature of the SM interactions can be neglected. In Section V, we finally incorporate gauge interactions into the simulations by modeling the SM with an Abelian-Higgs gauge model. Although this is just an approximation to the full gauge structure of the SM, the non-Abelian non-linearities can arguably be neglected. Therefore, the outcome of these simulations represent the most precise calculation of the dynamics of the SM after inflation, fully incorporating the non-linear and non-perturbative effects of the SM, while considering the gauge nature of its interactions. In Section VII we conclude and discuss some of the possible cosmological implications of our results.
All through the text = c = 1, m p 2.44 × 10 18 GeV is the reduced Planck mass, and t is conformal time. We take the flat FRW metric as ds 2 = a 2 (t)(−dt 2 + dx i dx i ), with a(t) the scale factor.

II. HIGGS OSCILLATIONS AFTER INFLATION
Let us characterise inflation as a de Sitter period with Hubble rate H * M EW , where M EW ∼ O(10 2 ) GeV is the EW scale. In reality, we know that inflation cannot be exactly a de Sitter background, since inflation must end after a finite number of efolds. The curvature perturbation spectral index n s = 0.968 ± 0.006, very well constrained by Planck to be smaller than unity at more than 7-sigma [17], is indeed interpreted as a signature of the quasi-de Sitter nature of Inflation. For our purposes, however, the distinction between de Sitter and quasi-de Sitter is irrelevant.
With a gauge transformation, the SM Higgs doublet can be parametrised in the unitary gauge by a single scalar real degree of freedom, Φ = ϕ/ √ 2. The renormalized improved potential for large field amplitudes |ϕ| M EW , is just given by the quartic part with λ(µ) the Higgs self-coupling at the renormalization scale µ = ϕ. Radiative corrections to the potential are encoded in the running of λ(µ), which to date has been computed to 3-loops when the Higgs is minimally coupled to gravity [7,8].
We ignore the nature of the sector responsible for inflation, so a priori there is no need for the Higgs to be coupled directly 1 to the inflationary sector. We will just consider that the Higgs field simply 'lives' on the de Sitter background, playing no dynamical role during inflation, and behaving simply as a spectator field [13,18,19]. As mentioned in Section I, the need to reheat the Universe after inflation requires somehow a coupling between the SM and the inflationary sector, though there is no particular constraint on this. Therefore, effective operators are expected to connect the Higgs with the inflaton when integrating out some possible mediator field(s). However, as we will show in the following sections, the Higgs decays very fast after inflation into all SM species. Hence, even if there is an inflaton-Higgs effective coupling, we will assume in practice that its effect is negligible, with the possible Higgs-inflaton interactions not affecting the Higgs dynamics during or after inflation.
Under these circumstances, the Higgs amplitude during inflation 'performs' a random walk at superhorizon scales, reaching very quickly, within few e-folds, the equilibrium distribution [20] P eq (ϕ) = N exp − 2π 2 3 The correlation length, i.e. the physical scale above which the Higgs amplitude ϕ fluctuates according to Eq.
(2), is given by l * ≈ exp{3.8/ √ λ} H −1 * [20], so it is exponentially larger than the inflationary Hubble radius H −1 * . After the equilibrium distribution is reached at some point during inflation, the correlation length remains invariant until the end of the exponential expansion. Hence, immediately after inflation, the Higgs amplitude ϕ can be safely considered homogeneous within any volume of size l l * . The Higgs amplitude varies randomly according to Eq. (2), but only when we compare it at scales l l * , much larger than the correlation length. For convenience, we define so that the distribution probability expressed over this dimensionless variable reads P eq ∝ exp[−2π 2 α 4 /3]. The cumulants of the distribution are then given by where ... denotes statistical average over the equilibrium distribution Eq. (2). One finds c 2 0.132, c 4 0.038, c 6 0.015, .... 1 Here we refer to a particle physics coupling, not the gravitational coupling whereas c 1 = c 3 = c 5 = ... = 0. A typical amplitude of the Higgs at the end of inflation is given by its root mean square (rms) value where we have defined the self-coupling normalized to a canonical value λ c ≡ 0.01, As we explain later, reasonable values of λ are taken within the interval [10 −2 , 10 −5 ]. Hence, for λ = 10 −2 , 10 −3 , 10 −4 , 10 −5 (λ 1/4 001 = 1, 0.562, 0.326, 0.178), we conclude that the typical Higgs amplitudes are of the order ϕ rms ∼ H * , independently of the value of λ. We find α ∈ [0.001, 1] with 99.8% probability, whereas α < 0.001 holds only with a 0.17% probability, and α > 1 is yet further suppressed with a 0.03% probability. We do not know the actual value of ϕ within the 'progenitor' patch from which our visible Universe grew up from. Actually we do not know the value of the Higgs condensate within any patch, we just know that typically ϕ/H * ∼ O(0.01) − O(1) for resonable values of λ. That means that just after inflation, within any patch of size l l * , the Higgs has a non-zero amplitude that could be really large, almost at big as H * depending on its realisation. The most updated upper bound for the inflationary Hubble rate is [17] H * ≤ H ). In order to analyse the dynamics of the Higgs after inflation, it is necessary first to fix the post-inflationary expansion rate. Since we do not specify the nature of the inflationary sector here, we can parametrise the scale factor after inflation like with a * the scale factor at the initial time t = t * (i.e. at the end of inflation), and w the equation of state of the Universe characterising the expansion rate of the period following inflation. For instance, if the inflationary sector is described by an inflaton with a quadratic potential, the Universe expands as in a matter-domination (MD) regime during the inflaton oscillations following the end of inflation, so w = 0 and p = 2. If it is an inflaton with a quartic potential, the Universe expands as in a radiationdomination (RD) regime, with w = 1/3 and p = 1. Since we do not really specify the inflaton sector, we are also free to consider other possibilities, including more 'exotic' scenarios where the background energy density decays faster than relativistic degrees of freedom, i.e. with w > 1/3 and p < 1. The paradigmatic example of this is a kination-domination (KD) regime, with w = 1 and p = 1/2, obtained when an abrupt drop of the inflaton potential takes place at the end of inflation, transferring all the energy into kinetic degrees of freedom [21,22].

A. Higgs Oscillations
The amplitude of the Higgs after the end of inflation is non-zero, and given that the Higgs potential is symmetric, the Higgs condensate is 'forced' to oscillate around its minimum at ϕ = 0. The larger the Higgs amplitude, the sooner the oscillations will start after the end of inflation. The eom (equation of motion) of the Higgs just after inflation is where · ≡ d/dt, and H is the comoving Hubble rate, given by We will consider the evolution of the Higgs in an arbitrary patch, inside which its amplitude ϕ * (randomly drawn from Eq. (2)) can be regarded as homogeneous.
The correlation length is exponentially bigger compared to the Hubble radius, l * e 38.2/λ 1/2 001 H −1 * H −1 * , so if we just follow the Higgs within a causal domain of initial size l ∼ 1/H * l * , then we can drop the Laplacian term in the rhs of Eq. (8). The only scale in the problem is therefore a * H * , so it is convenient to define a dimensionless conformal time z ≡ a * H * (t − t * ). We can then write the scale factor as a(z) = a * (1 + p −1 z) p . Introducing the variable with ϕ * the initial amplitude of the Higgs, we can rewrite the Higgs eom in a more convenient form as where ≡ d/dz, and β characterizes the frequency of oscillations. The term on the rhs scales as a /a ∼ (a * /a) 2/p , and hence it becomes irrelevant very soon, since it decays as a /a ∼ z −2/p 1. The initial condition for the Higgs amplitude in the new variables is, by construction, The initial condition for the derivative h * ≡ dh * /dz = 1+ ϕ(t * )/(a * H * ϕ * ), taking into account that the Higgs was in slow roll during inflation (i.e.φ(t * ) = −λa 2 * ϕ 3 * /2H * ,) reads out The initial velocity of the Higgs and the frequency of its oscillations (in the dimensionless variables) both depend, through β, on the initial random amplitude of the Higgs ϕ * and the actual value of λ. Therefore, at different patches of the Universe, (separated at distances larger than the correlation length l l * ), the Higgs will start oscillating with different amplitudes, and the oscillation frequency will also be different.
Depending on the amplitude of β, the Higgs will start oscillating around the minimum of its potential sooner or later. This can be clearly seen in Eq. (11), where the effective squared frequency of the oscillations of h(z) scales as ∝ β 2 . For the canonical value of λ = λ c = 0.01 (λ 001 = 1), the probability for the Higgs to start oscillating immediately at the end of inflation (i.e. that β ≥ 1) is extremely suppressed as 10 −287 %, being even smaller for λ < λ c (λ 001 < 1).
Therefore, at the end of inflation, the Higgs has an initial velocity in slow-roll and a non-zero amplitude as large as ϕ/H * ∼ O(0.01) − O(1) within any arbitrary patch of size smaller than l * . This amplitude remains 'frozen' for a finite time until the start of the oscillations. Looking at Eq. (8), and denoting as z osc (β) the time at which oscillations start at each patch, we see that the condition for the onset of oscillations is a(z osc ) √ λϕ(z osc ) = H(z osc ). For simplicity, we will set the initial value of the scale factor to unity a * ≡ a(t * ) = 1, so that H * ≡ H * , z ≡ H * (t − t * ), and a(z) = (1 + z/p) p . We will also denote any quantity evaluated at z osc with the suffix osc , so for example a osc ≡ a(z osc ). It follows that a osc √ λϕ osc = a osc H osc = H * /a 1/p osc , from where we find For a given expansion rate [p = 2, 1, 1 2 for MD, RD or KD], the period of oscillations depends sensitively on β, since the period is fixed when the oscillations condition a √ λϕ = H is attained at the time of the onset of oscillations z osc , which is itself a function of β and p. The time scale z M at which h(z) reaches its first maximum, characterized by h (z M ) = 0, also depends consequently on β and p. The period of oscillation can be easily obtained from the case of a field with quartic potential, initial amplitude ϕ * , zero initial velocityφ * = 0, and RD background. In conformal time, when the scale factor at the onset oscillations is set to unity, it is given by T = 7.416/( √ λϕ * ) [23]. In our case, we just need to count the oscillations from the first maximum at z = z M , taking into account that in our convention, a(z M ) = 1. The period, in units of z, is then found to be Let us note that the factor 7.416 is only exact for RD. For MD or KD, one expects a similar though somewhat different number, simply due to the term a /a in Eq. (11), which affects the very early stages of the Higgs dynamics (even if it decays very fast after the onset of oscillations). We have obtained fits for z osc , h osc , h(z M ) and Z T as a function of β and for each post-inflationary expansion rate (characterised by ω = 0, 1 3 and 1 for MD, RD and Evolution of the Higgs field for β = 10 −2 , 2.5 × 10 −2 , 5.0×10 −2 , 7.5×10 −2 and 10 −1 (corresponding to the red solid, orange dotted, blue dotted-dashed, green long-dashed and purple short-dashed lines, respectively). The background is RD, so w = 1/3 and p = 1. Dashed vertical lines mark the time zosc(β) when the oscillation condition is attained, a √ λϕ ≡ H, whereas continuous vertical lines mark the time zM (β) when the first maximum in the oscillations is reached, characterized by the condition h (zM ) ≡ 0. Top: Evolution of h(z). Lower: Evolution of the physical Higgs ϕ/ϕ * , which initially is frozen until the oscillations start, and then decreases as ∝ 1/a afterwards, as it oscillates. Similar plots are obtained for MD and KD backgrounds, whereas for other values of β the scale in the horizontal axis changes quite significantly.
KD, respectively). These fits will turn out to be useful later on. We find at the onset of oscillations, On the other hand, we find for the field amplitude at z = z M and the oscillation period (measured from z = z M onwards) where (A, B) (1.27, 6.28), (1.22, 6.33), (1.16, 6.41) respectively for MD, RD and KD.
At the end of inflation, the Higgs energy density at a given patch is mostly dominated by its potential energy, which represents a very small contribution of the total energy budget at that time, ρ * = 3m 2 p H 2 * . Averaging over realizations, we find At the onset of oscillations, part of the potential energy will become kinetic, with the two contributions -kinetic and potential -becoming of the same magnitude. In order to see this, we can average the Higgs field over its oscillations. The total energy density of the Higgs is given by It can be averaged over oscillations 2 as where the average can be splitted into potential and kinetic contributions, In Fig. 3 we can see the total energy density of the Higgs for different values of β, with and without averaging. Of course, as shown in the figure, the oscillation-averaged expressions are only valid when the Higgs has started oscillating, i.e. when z > z M . The figure also shows very nicely the fact that the averaged components verifȳ E V (β) = 1 3Ē (β) andĒ K (β) = 2 3Ē (β). Possibly, the most relevant aspect to be remarked is the well known fact that the Higgs energy density scales as a −4 with the expansion of the Universe, i.e. it behaves as a relativistic fluid. its oscillation-averaged valueρϕ/V * (dashed lines), for β = 10 −2 , 2.5 × 10 −2 , 5.0 × 10 −2 , 7.5 × 10 −2 , 10 −1 (from right to left, red, orange, blue, green, purple). Vertical grey lines mark zM (β), signaling from when the averaged curves should be considered valid. Bottom: The functions EK (z, β) (dashed, purple), EV (z, β) (dotted-dashed, blue) and E(z, β) ≡ EK (z, β) + EV (z, β) (solid, red), and their averagesĒK (β) (purple),ĒV (β) (blue) andĒ(β) ≡ĒK (β) +ĒV (β) (red), for β = 10 −2 . The figure shows thatĒV (β) = 1 2Ē K (β) = 1 3Ē (β). All plots obtained for a RD background.

III. HIGGS DECAY. ANALYTICAL ESTIMATES
As just explained, the Higgs oscillates everywhere in the Universe, although the time to start the oscillations depend sensitively on the initial condensate amplitude, which varies from patch to patch according to P eq (ϕ). Once the oscillations have began within a given patch, all fields coupled directly to the Higgs are excited every time the Higgs goes through the minimum of its potential. In the case of bosonic species this is known as parametric resonance, since a cumulative effect takes place and gives rise to a resonant growth of the number density of species [14,[25][26][27][28][29][30]. Although there is no parametric resonance in the case of fermionic species, yet an interesting effect occurs, since modes with successive higher momenta are excited as the oscillations carry on [28,[31][32][33][34].
All charged leptons of the SM are directly coupled to the Higgs via a Yukawa interaction, so all fermions of the SM will be excited during the oscillations of the Higgs [35], with the possible exception of neutrinos. Among the SM fermions the top quark has the largest coupling to the Higgs, so most of the energy transferred into fermions goes into top quarks [35]. More importantly, the SU (2) L gauge bosons are also coupled directly to the Higgs, and indeed the strength of their coupling is very similar to that of the Yukawa top quark. When two species, one fermionic and another bosonic, are coupled with the same strength to an oscillatory homogeneous field, the first burst of particle production is actually spin independent, and hence an equal number of bosons and fermions are created [36]. However, the successive particle creation bursts at each Higgs zero-crossing take place on top of an already existing number density of previously created species. The spin-statistics becomes then crucial, differentiating bosons from fermions in a noticeable way: bosonic occupation numbers start growing exponentially as the oscillations accumulate, whereas the fermion occupation numbers are always Pauli-blocked, forcing the transfer of energy into modes with higher and higher momenta. Both bosonic and fermionic excitations represent a sizable transfer of energy from the Higgs condensate. However, for equal coupling strength (as it is the case between top-quarks and SU (2) L -gauge bosons), the transfer of energy is much more efficient into the bosonic species [32]. Besides, in the context under study herethe decay of the Higgs after inflation -, the subdominant production of the SM charged leptons has been already adressed in [35]. Therefore, in this paper we will only focus on the production of the most energetically dominant species among the Higgs decay products, the W ± and Z gauge bosons.
In order to study the dynamics after inflation of the Higgs and its most energetic decay products, one should in principle consider the full SU (2) × U (1) gauge structure of the SM electroweak sector. However, one can make reasonable approximations both for analytical and computational purposes. In this work we have considered the following approximate schemes, mimicking the structure of the SM interactions: i) Abelian model. This consists on modeling the interactions between the SU (2) gauge bosons and the Higgs with an Abelian-Higgs analogue. Since gauge fields are initially excited by the Higgs from the vacuum, it is clear that non-linearities due to the truly non-Abelian nature of SU (2) are expected to be negligible during the initial growth of the gauge field occupation numbers [37]. The authors of Ref. [38] have shown that using the Hartreeapproximation, the effective contribution induced by the created gauge bosons onto themselves (due to the non-Abelian non-linearities) can be neglected as long as the backreaction from the gauge fields onto the Higgs does not become significant. In principle, this fact fully justifies to ignore the non-Abelian structure of the SM interactions, while maintaining only the abelian dominant part. ii) Global model. A more crude approximation can yet be done, by ignoring the gauge structure of the interactions. This does not mean that we ignore the interactions themselves, but rather that we consider them as if they were dictated by a global symmetry instead of a gauge one. In this scenario, one simply solves the mode equations of various scalar fields coupled to the Higgs with a quadratic interaction. Each of these scalar fields mimic a component of the gauge fields, with the quadratic interactions reproducing the coupling of the gauge bosons and the Higgs obtained from the SM gauge covariant derivative terms. This way, one can presumably capture the initial stages of the parametric resonance of W ± and Z bosons.
The approach i) is our most precise modeling of the SM interactions, but also the most involved one. We thus postpone its implementation for later on in Section V. The approach ii), though less accurate, has a clear advantage versus the gauge case: it allows not only a lattice implementation (which we introduce in Section IV), but also an analytical treatment (which we present in the remaining of this section). The analytical estimates represent only an approximation to the system described by the scenario ii), but yet provide a valuable insight into the understanding of the dynamics. The order of presentation in the paper of our different approaches is thus based on increasing progresively the degree of proximity to the real system. First, in the remaining of this section, we start with the analytical treatment of the global modeling, ignoring all non-linearities of the system. In Section IV we implement the global scenario ii) on the lattice. This way, we fully capture all non-linearities within this modeling, even if we yet neglect the gauge nature of the interactions. Finally, in Section V, we present a lattice implementation of an Abelian modeling of the system. This fully captures the non-linearities within such modeling, while preserving at the same time the gauge invariant nature of the interactions.

A. Analytical approach to the Higgs decay
In principle, we can follow the initial stages of the parametric resonance of the W ± , Z bosons by simply solving the mode equation for a scalar field χ, coupled to the Higgs with an interaction term of the form g 2 2 χ 2 ϕ 2 . Analytical results following this approach have indeed been presented in [19], so our work in this section should be understood only as complementary to such reference. We develope some formulas which will be useful later on, in order to assess the reliability of this analytical approximation when compared to the fully non-linear numerical lattice simulations.
Following Ref. [23], the equation for the Fourier modes χ k , after an appropriate conformal redefinition of the field χ, can be mapped into with q the resonance parameter, κ ≡ k/( √ λϕ osc ), ≡ d/dz, and z ≡ H osc t, a (dimensionless) time variable. This equation has a well understood structure of resonances: when q ∈ 1 2 [n(n+1), (n+1)(n+2)], n = 1, 3, 5, ... (i.e. q ∈ [1,3], [6,10], ...), there is an infrared band of resonance 0 ≤ k k * ≡ 1 √ 2π q 1/4 H osc [23], for which χ k ∝ e µ k z , with µ k bigger the smaller the k (therefore maximum at k = 0). If the resonance parameter q > 1 is not within one of the resonant bands, but lies in between two adjacent bands, then there is still a resonance of the type χ k ∝ e µ k z , but within a shorter range of momenta 0 < k min ≤ k k * , and therefore with a smaller Floquet index µ k . There is a theoretical maximum value for the Floquet index given by µ (max) k ≡ 0.2377... [23], so that any µ k is always contrained as µ k ≤ µ For simplicity, we will consider until the end of this section that the resonant parameter q = e 2 /λ always falls within one of the resonant bands, q ∈ [1, 3], [6,10], [15,21], .... As a matter of fact, in order to identify e 2 with the gauge coupling g 2 of the Higgs with the gauge boson (g 2 representing either g 2 Z or g 2 W ), we need to make the identification e 2 → g 2 /4. This matches correctly the interaction derived from the SU (2) gauge covariant derivative. The gauge couplings of the Z and W ± gauge bosons verify g 2 Z 2g 2 W ≈ 0.6 at very high energies, so q Z ≈ 2q W . Due to this relation, it is likely that either q W ≡ g 2 W /4λ or q Z ≡ g 2 Z /4λ fall within one of the instability bands. Let us note however that we cannot predict this, since the value of the Higgs self-coupling λ at high-energies is quite sensitive to the uncertainties in the Higgs mass m H , the top quark mass m t , and the strong coupling constant α s . Consequently, we cannot really know the exact value of these resonance parameters. However, taking into account first that the inflationary Hubble rate is constrained from above as H * 10 14 GeV [17], we know that in order to guarantee during inflation that the Higgs fluctuations remain below the critical scale µ + (above which the self-coupling starts decreasing, dλ/dµ ≤ 0), then the Higgs self-coupling value at high energies can only be within the range 10 −2 λ 10 −5 . Pushing λ to smaller values is in principle possible, but it represents a fine-tunning and requires some of the parameters m t , m H , α s to be more than 3-sigma away from their central values. We will consider therefore the range 10 −2 λ 10 −5 as the only acceptable one (with λ ∼ 10 −5 already only marginally valid). If beyond the SM physics affects the running, say stabilizing the Higgs potential at high energies, then λ remains positive and typically of the order λ ∼ 10 −2 − 10 −3 . Considering the range 10 −2 λ 10 −5 , and taking into account the strength of the W ± , Z gauge couplings at high energies, we obtain that the resonant parameters can only possibly be within the range O(10) q O(10 3 ). In particular, since g 2 = g 2 W 0.3 for W gauge bosons, we obtain q = 7.5 for λ = 10 −2 , and q = 3000 for λ = 2.5 × 10 −5 . For Z bosons we obtain similar resonance parameters, but twice as big. For completeness, we have sampled resonance parameters within the interval q ∈ [5, 3000], which corresponds to a range λ = 1.5 × 10 −2 − 2.5 × 10 −5 for W bosons and λ = 3.0 × 10 −2 − 5.0 × 10 −5 for Z bosons.
Let us then consider just one particle species A µ , representing either the Z µ or W ± µ gauge bosons, that will be parametrically excited during the Higgs oscillations. Let g 2 be the coupling stregth to the Higgs and let us represent the gauge field as if it was simply a collection of 3 scalar fields (one for each spatial component), all coupled with the same stregth to the Higgs. The growth of the fluctuations in the initial stages of resonance is described by the linearized Eq. (27). As long as the linear regime holds, even if the amplitude of the fluctuations grows exponentially, the use of 3 scalars should represent a good mapping of the real problem of gauge field excitation. Of course, one is ignoring this way the backreaction of the created bosons into the Higgs, as well as certain contributions in the gauge fields eom, which would be present if the gauge symmetry was restored.
The energy density of the created particles due to the resonance is then given by with the factor 3 accounting for the three spatial components of a gauge field, and where we have introduced an effective mass for the gauge boson, averaged over oscillations as For q 1, the maximum (comoving) momentum possibly excited in broad resonance, is given by from where we see (taking into account that h 2 ∼ h 2 osc ) that In other words, in broad resonance q 1, the decay products are always non-relativistic, and correspondingly we can approximate the effective mode frequency as ω k m A ∼ g  27) for several resonance parameters ranging between q = 5 and q = 3000. In each panel, we plot the corresponding the Floquet index µκ (where χκ ∝ e µκz ) as a function of momentum κ. We can divide the different q in two groups: those which contain a resonance at κ = 0 (blue lines) and those which not (purple dashed lines).
out that h rms h osc independently of β. 3 If q is within a resonant band, then all modes with momenta 0 ≤ k k * are excited with some Floquet index varying within [0, µ (max) k (q)]. This corresponds to the cases with blue solid lines in Fig. 4. We can therefore model the occupation number of the excited modes simply as a step-function n k = e 2µ k y where we have used that This is how the energy density of the gauge bosonsthose which are fully within a resonant band -will grow, at least as long as their backreaction into themselves and the Higgs is negligible. Using this linear approximation we can estimate the moment z eff at which an efficient trasfer of energy has taken place from the Higgs into the gauge bosons, characterized by ρ A (z eff ) = ρ ϕ (z eff ). This will be just a crude estimate of the time scale of the Higgs decay, since by then backreaction and rescattering effects start becoming important (hence invalidating our linear approach). However, the non-linear effects due to backreation/rescattering of the decay products simply tend to shut off the resonance, so using the linear regime for inferring the Higgs time scale should provide, at least, a reasonable estimate of the order of magnitude. More importantly, it provides the parametric dependences of both the time when the resonance is switched off, and the moment when the energy has been efficiently transferred into the gauge bosons.
The energy of the Higgs since the onset of the oscillations decays as where (h/h osc ) 4 ∼ O(1). We can now find z eff by simply equating Eqs. (32), (31), so that Let us recall that g 2 0.3, 0.6 and q ≡ g 2 /(4λ) ∼ O(10)−O(10 3 ), depending on the value of λ. Taking this into account, we find that the first term in the brackets of the rhs is always irrelevant , the second term is constant and of the order 9, and the last term is of order ∼ −1.
Therefore, we can approximate the above expression, using p √ a osc = (1 + 1 p z osc ), as Looking at Fig. 4, we see that the Floquet index of the modes within the resonant band 0 ≤ k k * (i.e. q's within a resonant band, corresponding to the blue solid lines of the figure) can be well approximated by a simple step function µ k ≈ µ k Θ(1 − k/k * ), with a mean Floquet index µ k 0.2. Taking this into account and using the fit of Eq. (16) for the time scale at the onset of oscillations z osc (β), we find The scale factor at z = z eff is then given by It is clear that depending on how small the initial value of β is within a given path of the Universe, the longer it takes for the Higgs to transfer energy efficiently into the gauge bosons. Since β rms ∼ O(0.1), we see that typically the Higgs decays at a time z eff (β rms ) ∼ O(10 2 ). Although the time varies from patch to patch depending on the values of β (with some dependence as well on the expansion rate through w), it is clear that the Higgs tends to decay really fast after inflation, within few oscillations. Later on, we will evaluate the validity of this estimate by comparating it with the outcome obtained directly from the lattice simulations.

IV. LATTICE SIMULATIONS, PART I. GLOBAL MODELING
In this section, we continue modeling the SM interactions with a set of scalar fields. We consider the Lagrangian with i = 1, 2, 3. Varying the action S = d 4 x L leads to the classical eom The term e 2 ϕ 2 χ i , under the identification e 2 = g 2 /4, mimics precisely the interaction term from the covariant derivative of the EW gauge bosons, g 2 2 Φ † ΦA µ , where A µ stands for either Z µ or W ± µ and Φ is the Higgs doublet. More concretely, choosing the unitary gauge for the Higgs Φ = (0, ϕ/ √ 2), and fixing A 0 = 0, we can identify χ i with each spatial component of the gauge boson A i , and the ϕ with the unitary representation of the Higgs. This way, by solving the system of scalar field equations Eqs (39)-(40), we can study the properties of the Higgs interactions with gauge bosons in an approximative way.
In Section III A we studied this scenario, following the fluctuations of the fields χ i with the help of the analytical solutions of the Lamé equation Eq. (27). We exploited the band structure of this equation and used some approximations in order to arrive at our analytical results, summarized in Eqs. (34)- (37). In reality, the scalar fields χ i follow the Lamé equation only initially, in the regime when the non-linearities (due to their small backreaction onto the Higgs) can be neglected. The fluctuations of the χ i fields grow exponentially during the linear regime, and as we will show, it does not take long until they start to impact onto the Higgs dynamics. At that moment, the system becomes non-linear, and only by following in parallel the coupled eom of the Higgs and χ i fields, we can really understand the field dynamics within this scheme. The aim of this section is, precisely, to solve numerically in a three-dimensional lattice the system of equations Eqs. (39), (40). Only in this way, we can fully capture the non-linear behaviour of this system beyond the simpler linear regime of the Lamé equation.
We now present the main results of the lattice simulations carried out for this scenario. We start with the following change of field variables, It is also convenient to redefine new spacetime coordinates z µ = (z 0 , z i ) with respect to the conformal ones To simplify notation, from now on and until the end of this section, all spatial derivatives will be considered to be respect this new set of variables. This change eliminates the friction terms in Eqs. (39), (40), and produces the equivalent set of equations with ≡ d/dz. A lattice version of these equations is presented in appendix A. As already mentioned, we will identify e 2 → g 2 /4, with g 2 either g 2 W or g 2 Z , i.e. the W ± and Z gauge couplings. The resonance parameter that appears naturally in Eq. (44), q ≡ e 2 λ , should therefore be interpreted as q ≡ g 2 4λ .
We have solved Eqs. (43) and (44) in three-dimensional lattices with periodic boundary conditions. We consider initial conditions given by a homogeneous Higgs mode (as described in Section II) and a null zero mode for the scalar fields coupled to the Higgs, We add, on top of the homogeneous contributions, a set of Fourier modes with a probability distribution such that |f k | 2 = 1 2a 2 ω k (in physical variables), which mimic the quantum vacuum fluctuations of the ground state of a scalar field in a FRW background. The Higgs is frozen in slow-roll until the oscillation condition Eq. (14) is attained at z = z osc (see bottom panel of Fig. 2). Hence, during the time 0 ≤ z < z osc , we only evolve in the lattice Eq. (43), corresponding to the slow rolling of the Higgs field (the homogeneous mode of the χ i fields is kept to zero). At z = z osc , we add the small inhomogenous Fourier modes to all fields, and from then on, we evolve together Eqs. (43) and (44). The reader can find more details about our methodology for introducing the initial conditions in Appendix B.
The actual value of λ is quite uncertain, since it depends on the energy scale of inflation. Besides, for a given Hubble rate H * , it can still vary significantly given the uncertainties in m H , α s and m t (mostly in the latter). Due to this, for each value of β, we have chosen a set of 26 resonance parameters q ≡ g 2 4λ , logarithmically spaced between q = 5 and q = 3000 (which corresponds to sampling the Higgs self-coupling from λ ∼ 10 −5 to λ ∼ 10 −2 ). Scanning this way β and q led us to characterise the behaviour of the system, scrutinizing all possible different outcomes depending on λ and ϕ * . In table IV, we list the values of all the resonance parameters q that we have considered. Since we cannot predict the value of the resonance parameters q, we have guaranteed that by sampling different values, we include both the cases in which q is within a resonance band of the Lamé equation, or in the middle of two bands (see Section III). Note that we have run simulations for three different expansion rates, corresponding to a MD universe, a RD universe, and a KD universe, given by ω = 0, 1 3 and 1, respectively in Eq. (7). The following results in this section will be presented only for a RD background, whereas the generalization of the results to other expansion rates will be considered in Section VI.
Our simulations depend only on two parameters, q and β. For each pair of values (q, β), we have run simulations on a lattice with N = 128 points per dimension, with periodic boundary conditions. The minimum momentum captured in each run is k m = 2π N dx , with dx the lattice spacing. The maximum momentum sampled in the lat- 2 k m . The length of the lattice box side is L = N dx. For each value of β and q, we have made sure that our results are not sensitive to the lattice spacing dx and/or the lattice size L. More details about these issues are given in appendix A.
In Fig. 5 we plot, as a function of time, the volumeaverage of the modulus of the (conformally transformed) Higgs field |h(z)|. In this figure, we show the outcome corresponding to β = 0.01, and four different resonance parameters: q = 8, 14, 101 and 354. The values q = 8 and 101 are centered close to the middle of a resonance band of the Lamé equation, while q = 14 and 354 are between adjacent bands. In this figure we also show the corresponding envelope curve of the Higgs oscillations. One conclusion is immediately clear: the time scale of the Higgs amplitude decay depends noticeably on q. By running simulations for each of the q values displayed in Table I, we have fully characterized the q-dependence of the Higgs decay. We have sampled all possible cases, including those where q is in the middle of a resonant band, close to the extremes of a resonant band, or simply outside of any band (between adjacent bands). Note that in Table I   menta range corresponds to the band with the highest Floquet index µ max , which coincides in all cases with the most IR one, see Fig 4. This index was obtained by solving the Lamé equation for a given q parameter, and finding the range of momenta such that µ(k) > 0. The band structure can be well appreciated in Fig. 4, where we plot µ(k) for each of the values of q listed in Table  I. As mentioned above, our range of resonance parameters sample both cases: when q is within a resonance band, and hence k min = 0, or when q is outside, and then k min > 0.
Before examining in more detail the general behaviour of all the fields in the system, we can make some comments about the Higgs behaviour. First of all, let us note that h oscillates with a period T which is, as expected, independent of the value of q. Even if it cannot be really appreciated in Fig. 5, we have checked that the period coincides initially with the analytical expression given by Eq. (15), until it becomes slightly modulated due to the interactions with the χ fields (though it does not change significantly). Looking at the different panels of Fig. 5, it seems that the Higgs decay is slower the greater the resonance parameter q is. This is very opposite to the intuition gained by the study of the Lamé equation in Sec- At later times z > zi, the amplitude of the Higgs drops abruply, due to its decay into the χ fields. In the lower panel we plot the physical Higgs |ϕ|/ϕ * = |h|/a, where we can appreciate that the plateau for h translates into a ∝ 1/a dilution for ϕ, due to the expansion of the universe. The decay of the Higgs into the other fields at later times is manisfested by a significant decrement of |ϕ| well below the 1/a decaying envelope.
tion III, which dictates that the larger the q, the shorter the decay time of the Higgs. 4 We thus see on this the first difference between the simplified study of the system of scalar fields in the linear regime (Section III A), and the real outcome when non-linearities are incorporated in lattice simulations. We will further comment on this issue later on. One can distinguish two different stages in each decay process. Let us look for instance at the upper panel of Fig. 6, where the Higgs modulus is plotted for q = 8, and where we also include the envelope curve of the oscillations (purple line in Figs. 5, 6). One can clearly appreciate that initially, and for some time, the envelope is approximately constant, reducing its amplitude only 4 In reality (and contrary to 'popular wisdom' about parametric resonance), the time scale x eff , identified with the oscillatory field decay time in the linear anaylitical approximation, is in practice mostly independent of q. It is true that the larger the q ther shorter the decay, and viceverse, but the dependence is logarithmically, so the number of oscillations does not change much.
slightly. This is observed as a plateau feature in the upper panel of Fig. 6. The vertical dashed line in the figure indicates the end of this initial behaviour, after which a second stage of rapid decay follows. Let us note that when we talk about the decay of the Higgs amplitude, we refer to the conformally transformed one h. The amplitude of the physical Higgs ϕ/ϕ * = h/a(t) is always decaying with the scale factor, no matter what its coupling to other species is. Before the second stage starts, the physical Higgs amplitude decays mostly due to the expansion of the universe, and not because of an efficient transfer of energy into the scalars. However, both effects are combined afterwards, producing an even more sharp decay in the Higgs amplitude. This is clearly seen in the lower panel of Fig. 6.
In order to understand better this two-stage behaviour, we plot the different contributions to the total energy of the system as a function of time. The energy density can be conveniently written as where, for our choice of variables, the Higgs and χ fields contributions to the kinetic (K) energy are given by (˙≡ the gradient (G) contributions are given by and finally, the Higgs potential (V) energy and the intereaction (int) term are given by In Fig. 7 we have plotted the different contributions to E t (z) for the parameters β = 0.01 and q = 8. Initially, the system is dominated by the kinetic and potential energy densities of the Higgs. This corresponds to the regime of anharmonic oscillations of the Higgs condensate described in Section II, for when the coupling to other fields was ignored (g 2 → 0). However, in reality, as soon as the Higgs starts to oscillate, there is an energy transfer into any species coupled to the Higgs. Each time the Higgs crosses zero, a small fraction of its energy goes into the χ fields. Initially, the amount of energy transferred at each zero crossing is small relative to the total energy stored in the Higgs. Therefore, it takes a time until the transfer becomes noticeable. The Higgs energy components represent the dominant contribution to the total energy during the initial oscillations, and the Higgs evolves without really noticing the presence of the other fields. Eventually, at the time z = z i , the energy transferred into the fields χ becomes significant enough as compared to the Higgs energy itself, say a fraction ρ χ /ρ ϕ = δ < 1. The Higgs condensate becomes affected by the transfer of energy into the other fields when δ ≥ δ s ≡ δ(z i ) 0.1. From then onwards, the Higgs continues pumping energy into the other fields at z > z i , but the amount of energy transferred at each zero crossing is no longer a small fraction of the energy available in the Higgs condensate itself. Therefore, soon after backreaction becomes noticeable at z = z i , the previously exponential growth of the χ fields energy densities stops, eventually saturating to a given amplitude. This is clearly seen in Fig. 7, where the gradient and kinetic energy densities of the χ fields saturate to an almost constant amplitude. This happens because the Higgs has not enough energy anymore to accomplish transferring a sizable fraction of energy into the χ fields. At the same time, immediately after z = z i , the Higgs energy density drops abruptly. This is so because the amount of energy transferred from the Higgs into the other fields, even if not significant anymore compared to the energy stored in the χ fields (hence the saturation of their growth), represents a significant fraction of the energy available in the Higgs at that moment. Therefore, the energy of the Higgs (mostly dominated by the kinetic contribution) drops abruptly, as it can be clearly seen, for instance, from z = z i ≈ 175 till z ≈ 900 in the case q = 8 depicted in Fig. 7.
Note that the Higgs energy density starts decreasing significantly at the moment in which the Higgs amplitude also starts decreasing. This is the moment z = z i , characterized by the fact that sufficient energy has been already transferred into the χ fields. However, the Higgs energy density eventually stops decaying, and its amplitude saturates to an almost constant value, as can be seen in the example of Fig. 7. The amplitude of the Higgs, instead, continues decreasing during a much longer time, as shown in Figs. 5, 6. The long lasting decay of the Higgs amplitude induces the decrease of the potential energy of the Higgs, even long after the dominant energy components of the Higgs have saturated. This can be clearly appreciated in Fig. 7, and is simply due to the fact that the Higgs keeps on oscillating, transferring energy into the χ fields. Since soon after z = z i the energy in the Higgs becomes smaller than the energy in the χ fields, the continuous transfer of energy represents only a marginal fraction of the energy already accumulated in the latter. Hence, the amplitude reached by the gra-dient and kinetic energy terms E χ K , E χ G is not affected anymore, whereas the amplitude of the Higgs potential energy continues decreasing. Eventually, the transfer of energy from the Higgs becomes inefficient and E V also saturates to an approximately constant value. By then, however, the Higgs potential energy is completely irrelevant compared to the gradient and kinetic counterparts.
A very relevant aspect to note is that when all the energy contributions stop growing or decreasing abruptly (with the exception of the Higgs potential energy, which keeps on falling for a long time), the energy components reach equipartition. In particular, somewhen at z > z i , the kinetic energy of the Higgs E ϕ K becomes equal to the sum of the Higgs gradient energy plus the interaction energy, E ϕ G + E int , see the lower planel of Fig. 7. In other words, equipartition in the Higgs sector holds as E ϕ K = E ϕ G +E int , and not as E ϕ K = E ϕ G . This could be expected to some extent, in analogy with a gauge invariant theory, as the one we analize in Section V. There we find that equipartition holds between the Higgs kinetic energy and the Higgs covariant gradient energy, which includes the interaction terms between the Higgs and the gauge fields (see for instance Fig. 13). In a gauge theory, contrary to what happens in the present scenario with only scalar fields, the gauge interactions cannot be isolated as individual terms in a gauge-invariant fashion. Therefore, in the gauge case, it is natural to consider together the interaction term + gradient term of the Higgs, since they are combined in a gauge-invariant manner through the gauge covariant derivatives. However, in the present case, there is no reason to consider the interaction term E int together with the gradient energy of the Higgs, unlike with the gradient energy of the χ fields. As a matter of fact, we have noticed that considering the sum of the gradient energy of the χ fields plus the interaction term with the Higgs, also "equipartitionates" with the kinetic energy of the χ fields. In other words, equipartition also works as E χ K = E χ G + E int . If, on the contrary, we only compare E χ K versus E χ G , we find that the former is slightly larger than the latter, as it can be appreciated in Fig. 7.
All features described so far are, of course, not specific to the particular case q = 8, β = 0.01 and RD shown in Fig. 7. A similar behaviour is observed in the outcome of the field distribution for other choices of β and q. That is, there is always first a plateau-like stage during which the Higgs (conformal) amplitude remains almost constant (or changes amplitude only marginally) for a few oscillations. Then, the amplitude decreases fast when the back-reaction from the χ fields becomes noticeable, which causes at the same time the cease of the exponential growth of the χ fields energy density. Eventually, all fields relax into a stationary distribution with exact equipartition E ϕ K = E ϕ G +E int and E χ K = E χ G +E int . On the other hand, E ϕ V becomes completely negligible as compared to any other energy term E i , in correspondence with the decay of the Higgs amplitude, which continues well after equipartition is set. The duration of the different stages, for a given expansion rate, is directly re-  We also add two new lines: the pink one, which corresponds to the sum of the Higgs gradient energy and the interaction energy, and the light blue line, which is the sum of the scalar gradient energy and the interaction energy. We see that the decay time indicates approximately both the time in which the Higgs kinetic energy has decayed, and the time in which equipartition is set. lated to the specific values of the parameters β and q. In particular, the duration of the initial plateau is directly dependent on the band structure of the Lamé equation of the corresponding q parameter. In the example shown in this section, parameters q = 8 and 101 have broad bands, with larger µ k and modes excited down to k = 0, whereas q = 14 and 354 have narrower bands, exciting modes only from k ≥ k min > 0, and with smaller Floquet indices.
We have characterized the dependence of z i with the resonant parameter q, see Fig. 8 and Eq. (56). Let us recall that this corresponds to the moment when the energy transferred into the χ fields is sufficiently large so that the Higgs amplitude and energy density starts to decrease. Therefore, this is the moment that should be compared to the analytical estimate of the Higgs decay time scale, Eq. (36), derived in Section III A. In light of the results of this section, we know however that the Higgs decay should rather be identified with the abrupt drop of the Higgs energy density, which takes place in reality just afterwards at z z i . After the drop, the kinetic contribution E ϕ K (which is the dominant energy component of the Higgs) enters into a stationary regime, equipartitioned with the gradient energy E ϕ G (as explained, in reality with E ϕ G + E int ). The onset of this regime signals the end of the decrease of the Higgs kinetic energy. We therefore define the time of the Higgs decay z e as the moment when equipartition within the Higgs sector holds. In other words, we define the Higgs decay into the other fields when the equality E ϕ K E ϕ G + E int holds better than 1%. This also provides an operational definition of the time z = z e , see Fig. 8. Although defining the Higgs decay like this might seems arbitrary, looking at the evolution of the different energy components we see that the end of the drop of the Higgs kinetic energy E ϕ K coincides practically with the onset of its equipartition with E ϕ G + E int . From then onwards, i.e. for z > z e , all energy If a given q is within a resonant band, z i (q) is almost independent of q, as appreciated in Fig. 8. As mentioned above, it is this scale z i (q) that should be really compared to the anaylitcal estimate z eff of Eq. (36), even if we know that in reality the time scale of the Higgs decay is rather z e (q), and not z i (q). For RD and β = 0.01, our anyaltical estimate Eq. (36) predicts z eff 200, which is reasonably similar to the fit to the numerical outcome z i (q) ≈ 160 shown in Fig. 8. The analytical estimates are only an approximation to the real dynamics and one cannot expect anything more than a reasonable order of magnitude prediction, as it is indeed the case. More importantly, the analytical calculation predicts that z eff should be only dependent on q logaritmically, which in practice implies that for mildly borad resonance parameters q ∼ O(10) − O(10 3 ) (within a resonance band), z eff is essentially independent of q, as indeed well appreciated in Fig. 8. The dependence of z i (q) with q's outside resonance bands is only logarithmic, but with a big coefficient. As it can be appreciated in the upper panel of Fig. 8, for q 10 2 it is a factor ∼2-4 larger than the analytical prediction Eq. (36), but becomes of the same order for q ∼ 10 2 − 10 3 , modulo a factor ∼1-2. Possibly, for q 10 3 , z i (q) will become smaller, but as said before, such regime is never valid in our case of study.
A fit to the decay time scale z e (q), plotted in the lower panel of Fig. 8, is given by This behaviour is independent of whether q is within or outside a resonance band. We believe this is because z e (q) is set by the non-linearities of the theory (as opposed to z i (q), which is determined by the linear regime of the system). Let us describe here the Higgs decay in the spectral domain. During the initial stages, the modes that are excited correspond to those in the band structure of the Lamé equation for each of the scalar fields. We clearly see this for z < z i (q) in Fig. 9, where we plot both the field spectra k 3 |X k | 2 and its occupation number k 3 n k . We also indicate with dashed lines the resonance bands. As the amplitude of the modes within the resonance bands grows, the system becomes more a more non-linear. Rescattering among modes occur and the bands become wider. Due to the coupling of the modes through Eqs. (43), (44), the initial parametric resonance of the χ k modes within the resonance bands excite at the same time Higgs modes ϕ k , which then rescatter off other modes of χ k , and so on. As a consequence, the power spectrum of the fields grows exponentially and widens, with a typical width 0 ≤ k O(10)k * . At late times, the fields enter into a quasi-stationary turbulent regime, characterized by a very slow evolution of the energy densities [39,40]. This corresponds in the spectral domain to the distributions entering into a slowly selfsimilar regime, with momenta gradually increasing. Note that total thermalization is still a long way off, since this turbulent regime is not very efficient in transferring energy from the long wave modes to the high momentum region.
Our analysis in this section will prove useful for the understanding of some of the properties of the Higgs decay in the gauge simulations of Section V. It will also serve to ascertain the similarities and differences between the two approaches. Let us recall again that all our results of Section IV, correspond to RD and were obtained for a fixed value β = 0.01. We will discuss in Section VI how the outcome changes when varying the Higgs initial amplitude (i.e. β) and the background expansion rate (i.e. w).

V. LATTICE SIMULATIONS, PART 2: ABELIAN-HIGGS MODELING
In this section, we study the properties of the Higgs decay by modeling it with an Abelian-Higgs model. In this approach, and in contrast with the global scenario studied in the last section, we introduce for the first time a gauge structure to the interactions of the Higgs with the gauge bosons. The differences and similarities in the results of these two scenarios will be scrutinized. In this section, we approximate the action of the electroweak sector of the Standard Model, invariant under the local SU(2)×U(1) symmetry group, by a local U(1) gauge theory. This is justified in principle because, as it is shown explicitely in Section V A, we do not expect the corrections due to the truly non-Abelian nature of the SM interactions to play any significant role in the dynamics of the Higgs decay, at least on the onset of the characteristic time scales.
The Abelian-Higgs model is described by the action where the covariant derivative and field strength are defined as usual, as Here, e is the Abelian coupling strength representing the coupling of either the W ± or the Z gauge bosons. As in the last section, in order to mimic correctly the Higgsgauge interactions, we need to take e 2 = g 2 /4, with g 2 = g 2 W or g 2 = g 2 Z for the W ± and Z bosons respectively. Note that we are working with a system with local U (1) gauge symmetry, and consequently, we must take the Higgs field as complex. We write it in terms of its components as From Eq. (57) we derive the following equations of motion As we are dealing with a gauge invariant theory, we have a gauge freedom in the choice of our equations. This allows us to set, from now on, the condition A 0 = 0. In this case, the gauge eom, Eq. (61), can be written in terms of its components as Eq. (63) is the Gauss law, which represents a dynamical constraint that the solution to Eqs. (60) and (62) must preserve at all times. When solving them in a 3dimensional lattice, one must of course check that the Gauss constraint Eq. (63) (or more specifically, its equivalent discrete equation) is indeed preserved during the whole evolution of the system. We also define the gaugeinvariant electric and magnetic fields as usual As in the global scenario, it is really useful to do a change to other spacetime and field variables, which we call "natural" as before. On the one hand, we change to the same set of spacetime coordinates z µ = (z 0 , z i ) introduced in Section IV, which are On the other hand, it is also convenient to define new Higgs and gauge field variables as (with j = 1, 2; i = 1, 2, 3) where ϕ * ≡ |ϕ(t * )| is the initial absolute value of the Higgs field at the end of inflation.
To distinguish between different times, we use a dot or a comma to denote differentiation with respect conformal or natural variables respectively (˙≡ d/dt, ≡ d/dz). Also, and from now on, all spatial derivatives will be with respect the new variables unless otherwise stated. We also define, from now on a dimensionless covariant derivative as With these changes, Eqs. (60), (62), (63) can be written as where the current j µ (x) is defined as The convenience of this new set of variables Eq. (66) is that the a(t)-dependence in the gauge fields' eom disappears. Finally, we define new dimensionless electric and magnetic fields as In this work, we have solved the system of Eqs. (67)-(70) in different 3-dimensional lattices. More specifically, we have solved a gauge-invariant set of analogous equations in a discrete spacetime. In all simulations, we have ensured that the lattice analogue of the Gauss conservation law Eq. (70) is only marginally broken by the timeevolution of the system. The reader can find more details of our lattice formulation in appendix A. We have considered the following initial conditions for the homogeneous modes of the fields. First, let us define |h| ≡ h 2 1 + h 2 2 . From Eq. (66), we have by construction |h * | ≡ |h(t * )| = h 2 1 * + h 2 2 * = 1 at the end of inflation. As long as this condition is satisfied, we can freely distribute this initial value between the components h i * ≡ h i (t * ) due to the symmetries of the Higgs potential. A convenient choice is On the other hand, as we are evolving the system of equations since the end of inflation, the Higgs initial velocity must obey the slow-roll condition (φ i (t * ) = −λa 2 * ϕ 2 * ϕ i /2H * ). With the choice of Eq. (73), the slowroll condition reads, in natural units, We also set the homogeneous mode of the gauge bosons to zero, V i * = V i * = 0, until the onset of the oscillations at z = z osc . The system is solved in the following way. First, for the times 0 < z < z osc , we only evolve the Higgs field equations with the gauge bosons set to zero. At z = z osc , we add fluctuations on top of the homogeneous modes of the different fields, allowing the gauge boson production to take place. Over the homogeneous mode of the Higgs components, we add Fourier modes with a spectrum such that |f k | 2 = 1 2a 2 ω k (in physical variables), which represent the vacuum fluctuations of the ground state of a scalar field in a FRW background. However, due to the gauge nature of the system, the initialization of the gauge fields is more subtle and delicate than in the case of scalar fields. In this case, the fluctuations we add to the gauge fields must preserve the Gauss constraint Eq. (70) initially at every lattice point. Thus, given the spectrum of Higgs fluctuations, we fix the gauge fluctuations as given by the right hand side of Eq. (70). More specifically, we fix the gauge fields amplitude in momentum space as where in the lattice this is done with the corresponding lattice momenta (see [41,42] for a discussion) corresponding to the choice of lattice finite difference operators to mimic continuous derivatives. The implementation of these initial conditions is described in more detail in appendix B. In particular, we discuss the importance of setting appropriately the Higgs fluctuations so that we can impose correctly Eq. (75). From z ≥ z osc onwards, the Gauss law is then preserved by the gauge invariant evolution of the system. How this is checked is discussed in more detail in section A. We now present the main results of the lattice simulations carried out for the Abelian-Higgs model. We will compare the outcome of these simulations with the results of the global scenario presented in Section IV. Like in the case with only scalar fields, we have run simulations for the resonance parameters given in table IV, which range from q = 5 to q = 3000. These values correspond to λ values between 2.5 × 10 −5 and 1.5 × 10 −2 for the W boson, and 5 × 10 −5 and 3 × 10 −2 for the Z boson. We have also runned different simulations for β = 10 −4 , 10 −3 , 10 −2 , 10 −1 and 0.5. The justification of this choice of parameters has been explained in detail in Section IV, so we do not need to enter in further detail here. Again, all results presented in this section will be obtained for a RD background (w = 1/3) and for the β = 0.01 value. We explain in Section VI how these results can be easily extrapolated to other ω and β parameters.
One of the main differences of the Abelian-Higgs model with respect to the global scenario, is that now the Higgs field is described by a set of two real components h 1 and h 2 , combined in a complex variable h = h 1 + ih 2 . This means that the quantity of interest that we must study is the average value of the Higgs absolute value, |h| ≡ h 2 1 + h 2 2 . Note that in the global case of Section IV, we analyzed the analogous quantity |h|, corresponding to the absolute value of the only component of the Higgs. Therefore, the results presented here about the decay of the conformal Higgs amplitude can be easily compared with those of the global scenario. We have plotted in Fig. 10 the volume-average of the (rescaled) Higgs modulus |h| as a function of time and for two specific resonance parameters, q = 9 and q = 23. We have plotted in purple the corresponding oscillations' envelope curve by joining all local maxima with a smooth line. Remember that, according to what we discussed in Section III, all resonance parameters can be divided in two groups: those placed in the middle of a resonance band of the Lamé equation, with an interval of excited momenta of the type 0 ≤ k k max , and those which have a smaller band of the type k min k k max . Examples of these two groups are drawn in blue and purple respectively in Fig. 4. In this regard, q = 9 belongs to the first group, and q = 23 to the second. The initial period of oscillations, before the amplitude of the Higgs drops significantly, fits well the analytical estimate of Eq. (15). This is expected even in the present case with a complex field, since before the Higgs notices the presence of the gauge fields, the dynamics is still effectively equivalent to that of a single real degree of freedom. When the Higgs amplitude (and energy) starts decreasing due to its transfer of energy into the gauge bosons, the period of oscillations is slightly modulated, but never significantly.
For relatively mild broad resonance parameters q 200, we find that the Higgs amplitude behaves in a qualitatively similar way as in the global scenario. This can be rapidly seen by comparing Fig. 10 to the equivalent   6 from the global scenario. In both scenarios, there is first a stage of few oscillations during which the (conformal) Higgs amplitude does not decay, corresponding to the plateau in the envelope functions. After that, at a time z z i (q), the Higgs amplitude starts to decay strongly. This time is indicated in both panels of Fig. 10 with a red, dashed line. Finally, the rescaled Higgs amplitude approaches a constant at late times, |h| → h f (q), which is indicated in both figures with an orange, dashed line. It is important to emphasize that the plateau is only manisfest for the conformal Higgs field |h|, because the physical Higgs field |ϕ|/ϕ * decays as ∼ 1/a(t) due to the expansion of the Universe. The key observation is that, for z z i (q), it decays as the inverse of the scale factor ∝ 1/a, while for z z i (q) the decay is much faster due to the energy transfer to the gauge fields. This can be clearly seen in Fig. 11, shown for q = 23. Therefore, we conclude that for mild resonance parameters q 200, the qualitative behavior of the system is very similar, almost identical, to the global scenario.
For q 200, a new phenomenon takes place, due to the complex nature of the Higgs. We have plotted in Fig. 12 the value of |h| as a function of time, for the resonance parameters q = 214 and q = 584. One can observe, first of all, that the Higgs oscillations do not present a clear initial plateau like for smaller resonance parameters. It is like if, in these cases, z i ≈ 0, with the Higgs decaying into the other species from the onset of oscillations. At the same time, we observe that the amplitude of the oscillations is modulated by another oscillatory function. This behaviour is particularly well observed in two of the specific cases we have studied, q = 214 and q = 275 (which correspond to λ = 3.5 × 10 −4 and λ = 2.7 × 10 −4 for the W ± bosons respectively). This behaviour was not observed in the simulations carried out in the global modeling, so it represents a clear distinction between the global and gauge lattice approaches. It indicates that, in order to capture all the relevant phenomenology of the Higgs decay, we really need to work in a gauge-invariant framework. However, it is true that this particular modulation has no practical consequences for the Higgs decay, since z e (q) still remains very similar to its global counterpart, independently of whether this modulation is present (q 200) or not (q 200), see lower panel of Fig. 14. Besides, this modulation is characteristic of the Abelian-Higgs scenario, so in the non-Abelian real case the modulation could be different. 5 The time scale z i (q) signals, as in the global modeling, the moment at which the decay products (in this case, gauge bosons) have accumulated sufficient energy to start affecting the dynamics of the Higgs condensate. As before, this is understood better if we plot the different contributions to the energy, as a function of time. The energy density of the Abelian-Higgs model Eq. (57) is found to be where V * is the value of the Higgs potential at the end of inflation. The function E t (z) is formed by the sum of 5 Because of this, we do not really attempt to model or explain the details of the modulation. After all, this type of effects are generically expected in 2-dimensional oscillators, as the one we have in the Higgs-Abelian model. the following contributions, Here, E K and E V are the kinetic and potential energies of the Higgs field, E GD is a gauge invariant term formed by the product of two covariant derivatives of the Higgs field (hence containing the spatial Higgs gradients plus the interaction terms), and E E and E M are the electric and magnetic energy densities We have plotted in the upper panel of Fig. 13 these quantities as a function of time for the resonance parameter q = 9, which corresponds to a case within a resonance band. We also show in the lower panel of Fig. 13 the contribution of each energy component to the total, E i /E t . In this lower figure we have also removed the oscillations of each component, showing only the corresponding envelope functions. We see that initially the dominant contributions come from the kinetic and potential energies of the Higgs field, corresponding to the oscillations of the condensate around the minimum of its potential, before it 'feels' the gauge fields. Meanwhile, the other components of the energy, E E , E M and E GD , grow really fast, due to the energy transfer from the Higgs into the gauge fields. Note that for the whole evolution of the system (until equipartition is reached), the electric energy clearly dominates over the magnetic energy.
As in the global analogue, although gauge bosons are being strongly created, the Higgs condensate is at first unaffected. Again, at z ≈ z i (q) (indicated by a dashed red vertical line in the figures) the gauge energy has grown enough to start affecting significantly the Higgs condensate. The consequence of this is a sharp decrease of both the Higgs potential and kinetic energies. Physically, this happens when δ ≡ E E /E t becomes sizeable. When δ 0.1, we can clearly see the correspondence between the backreaction of the gauge fields over the Higgs field and the decrease in the Higgs amplitude.
As in the global scenario, for z z i (q) the Higgs kinetic and potential energies start decreasing sharply. The potential energy becomes very soon irrelevant compared to the other energy contributions, while the kinetic energy approaches an almost constant amplitude. Simultaneously, E GT and E E stop their decay, and also saturate to almost constant values. However, the magnetic energy continues to grow even after E GT and E E have estabilized. Finally, at z = z e , the system arrives again to a stationary regime, in which equipartition between different components is clearly achieved. In this regime, 30% of the total energy goes to the Higgs kinetic energy, 30% to E GD , 20% to electric energy E E , and 20% to magnetic energy E M . The potential energy E V also saturates to a constant, but very subdominant with respect to the other contributions. Quite remarkably, these numerical percentages are independent of the values q and β taken in our simulations. In other words, the final fraction of energies are universal within the Abelian-Higgs formulation. We also expect this to be the case in the real non-Abelian model.
Let us analyze this equipartition regime in more detail. In the gauge scenario, the kinetic energy of the Higgs field E K eventually becomes equal to E V + E GD . Note that since E GD is gauge invariant, it contains both the Higgs gradient terms plus the Higgs interactions with the gauge fields. One can then naturally identify this quantity in the global scenario as the sum of the interaction term plus the Higgs gradients, E int + E ϕ G . Interestingly, the potential energy keeps decaying even after equipartition has been established, similarly as in the global scenario. In principle, therefore, we could think of using just the equipartition relation E K E GD , neglecting the contribution of the potential energy, as we did in the global case. However, in the moment when the rest of energy contributions are stabilized, the potential energy represents still a 1% of the total. So, although the potential energy becomes eventually subdominant, its rate of decay is slower than in the global scenario. This percentage, although small, is still significant at the moment in which equipartition is achieved. Therefore, it is better to follow the equipartition condition E k E V + E GD . The evolution of the different energy components and the achievement of equipartition can be well appreciated in Fig. 13.
It is useful to define the Higgs decay time as the moment at when the Higgs kinetic energy (which dominates over the potential energy) results stabilized at the onset of the stationary regime. As in the global scenario, we will call this quantity, z e (q). Naturally, there is again some degree of arbitrariness in this definition. In the global scenario, we observed that a good operative criterium for defining the z e decay time was based on the degree of equipartition achieved. In our present gauge context, we have observed that an appropriate criterium is to take the moment when the relative difference between E K and the sum E GD + E V becomes less than 1%. We have indicated this time in Fig. 13 with a dashed, vertical line.  . Note that we have extracted a global scale factor term in (76), which means that the total energy of the system is constantly diluting even if these functions seem to go to constants.
As we mentioned above, in the global scenario we did not consider the contribution E V of the Higgs potential energy in the equipartition equalities, simply because its contribution was already irrelevant when equipartition was reached. However, in the Abelian-Higgs scenario, the addition of this contribution to the covariant-gradient one E GD is crucial. Even though E V is also marginal in this case, if we were to consider just E GD in the equipartition analysis, it would achieve equipartition with E K (say to better than 1%) long after the Higgs kinetic energy density has started to saturate. As we can observe in Fig. 13, our criterium E K E GD + E V holding better than 1%, coincides very well with the moment at which all relevant energy densities have stopped either growing or decreasing. Hence itdefines very well what we mean by the end of the Higgs decay.
We have then defined two characteristic time scales, z i and z e , similarly as in the global case. The time z i indicates the moment when the backreaction of the gauge bosons into the Higgs condensate causes the start of its decay (both in amplitude and energy). The time z e characterizes the ending of the decay, which signals also the onset of the equipartitioned stationary regime. As before, we have characterized the dependence of z i and z e with the different q's considered. Let us start with z i . We show in the upper panel of Fig. 14 the different z i (q) obtained for the different resonance parameters with q 200 (remember that for larger values q > 200, the plateau is not present and the Higgs amplitude and energy density start decaying from the very onset of the Higgs oscillations). In the figure, blue squares correspond to q values within a resonance band, and purple circles correspond to values outside bands. We see a clear trend, such that simulations with q within resonance bands have a smaller z i (q) than those with q between adjacent bands. This is, however, only a trend, in the sense that it is not a robust fact as it was in the global scenario (recall Fig. 8). In fact, as seen in Fig. 14, there are a couple of purple points which do not fulfill this property. Given the disparity in the behaviour of z i versus q, we do not attempt to fit it. We simply content ourselves with the observation that, for q (q 200) within a resonance band, the analytical prediction from Section III A still provides a correct estimate of the order of magnitude of z i .
In the lower panel of Fig. 14, we also plot z e as obtained for our different resonance parameters q. We have obtained the following phenomenological fit z e (q) = 588q 0.42 , indicated in the figure with a red line. Note that we have plotted as well the corresponding fit obtained from the global simulations, Eq. (56). Both fits coincide pretty well, indicating that the Higgs decay time z e (q) obtained in the global scenario constitutes already a very good estimation. To some extent this is surprising, since one could expect that the extra terms in the gauge fields eom could play some role, like for example modulating the decay time z e (q) differently than in the pure case of scalar fields. However, our results prove that this is not the case.
In fact, they imply that the interaction term g 2 A µ A µ ϕ 2 (which is the only one kept in the global scenario) is the most relevant one when determining the end of the Higgs decay and the onset of the stationary regime. An alternative source of information about the Abelian-Higgs system comes from the spectra of the different fields. Since we are dealing with a gauge theory, all quantities of physical interest must be gauge invariant. We then plot in Fig. 15 the spectra of the electric and magnetic fields k 3 |E k | 2 and k 3 |B k | 2 at different times. In order to see the dependence of the spectra evolution on the analytical properties of the Lamé equation, we plot both spectra for two different resonance parameters, q = 5 and q = 9. The latter is placed in the middle of a resonance band, while the former is between the first and second resonance bands. The dashed vertical lines in the figures indicate the location of the respective resonance bands. In the case q = 5, one can clearly see that both spectra grow with time, as a consequence of the resonant excitation of gauge bosons. At initial times, there clearly appears a peak in both spectra, centered in the corresponding main resonance band. This confirms that the behaviour derived from the Lamé equation describes well enough the real dynamics during the initial stages, even for the gauge theory. When the gauge bosons start to affect significantly the Higgs condensate, i.e. for z z i (q), both spectra start to displace to the right, populating modes of higher momenta. In this process, new subdominant peaks appear. As times goes on, the peaks disappear, and when the Higgs condensate has decayed (i.e for z z e (q)) the stationary state is established. For the case q = 8, the time scale z i (q) ≈ 150 is much smaller than in the previous case, and the resonance band is much wider. This is expected, as we include modes down to k = 0. In this case, we see that the population of higher modes is much faster than in the latter case, and we do not observe additional subdominant peaks in the spectra.

A. Beyond the Abelian-Higgs
The real nature of the SM interactions is non-Abelian, since the EW sector of the SM is SU (2) × U (1) gauge in- variant. In the eom of the gauge bosons there are therefore non-linear terms 6 of the form ∼ g 2 A 3 , gA∂A, g∂A 2 , where we omit charge and Lorentz indices for simplicity. Following [30], one obtains that within the Hartree approximation, the terms ∼ gA∂A, g∂A 2 vanish, so that in principle only the terms ∼ g 2 A 3 contribute effectively to the dynamics of the gauge fields. We can write the effective mass entering in the gauge fields eom, as given by their interactions with the Higgs plus a contribution from their own non-Abelian self-interactions. Symbolically we write this as The Abelian-Higgs simulations capture the first term g 2 ϕ 2 , which is due to the interaction with the Higgs and is responsible for the resonant excitation of the gauge fields. The self-induced mass due to the gauge field selfinteractions is, of course, not present in the Abelian approach. This second term describes the non-linearities of the non-Abelian nature of the SM interactions, and hence only when the gauge fields have been excited with a sufficiently high amplitude A 2 g 2 ϕ 2 , may its presence have any relevance. The question is then, when do the gauge fields reach the critical amplitude A ∼ A c ≡ gϕ?
The answer can be easily found by analyzing the effective mass of the Higgs. The non-Abelian nature of the interactions does not add any extra contribution into the effective mass of the Higgs field, given by These terms are already captured in our simulations, so the only difference in a non-Abelian simulation would come from the fact that A µ is affected by the nonlinearities of its own eom. The gauge fields backreact into the Higgs dynamics at the time z = z i (q), which corresponds physically with the moment when the amplitude of the gauge fields has grown -due to parametric resonance -up to A 2 λϕ 2 /g 2 . This condition corresponds however to a typical amplitude of the gauge fields A ∼ A(z i ) ≡ √ λϕ/g, which is much smaller than A c . In particular, A(zi) Ac ∼ 1 g √ q < 1 for the typical broad resonant parameters q ∼ O(10) − O(10 3 ). The effective mass of the gauge bosons at z ≈ z i is

It is then clear that m 2
A (z i ) ≈ g 2 ϕ 2 , like if there was no effect from the gauge-field selfinteractions. By the time the gauge field resonant production backreacts on the Higgs dynamics, the gauge fields amplitude stop growing, as explained in detail in Section V. Therefore, we do not expect the non-Abelian terms (neglected in the Abelian-Higgs approach) to play any significant role in the dynamics of the system. 7 It is however likely that the presence of the non-Abelian terms will possibly change the details of the achievement of the equipartition regime. We do not expect however the time scale z i (q) to change at all, whereas the time scale might change perhaps slightly. Only non-Abelian lattice simulations beyond our present work can quantify this.
In light of this analysis, we see a posteriori that neglecting the non-linearities due to the non-Abelian nature of the SM interactions is well justified. 7 It is possible though that for the mildest broad resonance parameters, such as q ≈ 10, there might be some effect from the non-Abelian terms, since in this case the 1 g 2 q correction in Eq. (83) is only marginally smaller than unity.

VI. VARYING THE HIGGS INITIAL AMPLITUDE AND THE EXPANSION RATE
All results from sections IV and V have been presented for a scale factor evolving in a RD universe (ω = 1/3), and for β = 0.01. Naturally, in order to fully understand the dynamical properties of the Higgs decay after inflation, we have run the simulations for other β parameters, and we have also considered other expansion rates such as MD (ω = 0) and KD (ω = 1). Fortunately, one can easily extrapolate the results from one particular set of parameters, say (β 1 , ω 1 ), to another set (β 2 , ω 2 ), using the analytical properties of the Higgs equation described in Section II. In other words, from the results obtained for (β 1 , ω 1 ), we can obtain an approximation to the ones for (β 2 , ω 2 ).
More specifically, we saw in Eq. (18) that in the case of no coupling to the gauge bosons, the conformal period Z T and the value of the conformal Higgs field at the first maximum h(z M ), can be approximated as Here c 1 and c 2 are constants, which are approximately independent of ω and β. From these properties we can see that, if for a given set of values (ω 1 , β 1 ), the volume-averaged Higgs field takes the value h(β 1 , ω 1 ) at the time z(β 1 , ω 1 ), then for (ω 2 , β 2 ) the Higgs field at the time should take the value Notably, this property is maintained quite well even in the presence of a Higgs coupling to its decay products (either scalars in the global simulations or gauge bosons in the Abelian-Higgs simulations). This extrapolation is therefore very powerful. In Fig. 16, we have plotted the volume-averaged value of |h| as a function of time, for both global (top figures) and Abelian-Higgs simulations (bottom figures). Let us focus for instance in the top, left figure. We have obtained for q = 8 the behaviour of |h| as a function of time for β = 10 −4 , 10 −3 , 10 −2 , 10 −1 , and 0.5, directly from the simulations. Using the outcome from these simulations with different β parameters, we have then inverted Eqs. (84) and (85), and obtained the (extrapolated) behaviour corresponding to β = 0.01. These are different predictions for the Higgs decay when β = 0.01, but obtained from the real data from simulations with different β values. We see that the four different extrapolated theoretical predictions obtained for β = 10 −4 , 10 −3 , 10 −1 and 0.5 coincide very well with the real simulation for β = 0.01.
The same is done in the top right figure, but changing the scale factor instead of β (which we fix in this figure as β = 0.01). There, we compare the result of the Higgs decay for a RD universe, on one hand obtained directly from simulations with ω = 1/3, and on the other hand the corresponding extrapolated predictions from the lattice simuLations with ω = 0 (MD) and ω = 1 (KD). The three lines also coincide very well. The same analysis is repeated for Abelian-Higgs simulations in the two bottom figures, with identical conclusions.
This property allows us to extrapolate easily the results for a RD universe with β = 0.01 presented in the last two sections, to any other (ω, β) parameters. In particular, the fit for the Higgs decay time found in Eq. (80) is generalized to an arbitrary (β, ω) as

VII. SUMMARY AND DISCUSSION
The recent measurements of the Higgs boson mass imply a relatively slow rise of its effective potential at high energies. In the regime where the EW vacuum is stable with the Higgs self-coupling kept positive, the Higgs develops a large VEV during inflation, representing a classical condensate, homogeneous over scales exponentially larger than the inflationary radius 1/H * . In this paper we have studied the relaxation of the Higgs, i.e. its decay, during the stages following immediately after inflation. In reality, the origin of the VEV during inflation, which sets up the initial condition for the decaying process, is not particularly relevant for our analysis. If another mechanism (different than quantum fluctuations) is responsible for the development of the Higgs VEV during inflation, our calculations and results would be equally applicable. The case considered in the paper, with the initial amplitude of the Higgs condensate dictated by the equilibrium distribution Eq. (2) due to the stretching of its quantum vacuum fluctuations, simply serves as a starting and practical point, to assess the typical Higgs amplitudes at the end of inflation.
The decay of the Higgs condensate during the early post-inflationary stages constitutes an important event in the evolution of the Universe, which might have interesting cosmological consequences. In this article we have focused on the details of the Higgs decay process itself. We have used different methods of progressive complexity, accuracy and proximity to the real case of the SM. We have modeled the SM interactions in a two-step manner. First, considering a global scenario, ignoring the gauge structure of the SM, representing the gauge fields as a collection of scalar fields appropriately coupled to the Higgs. Secondly, we have considered an Abelian gauge scenario, with the gauge fields and the Higgs embedded within an Abelian-Higgs framework, ignoring the nonlinearities due to the truly non-Abelian nature of the SM. With the global model we have presented both analytical (Section III A) and lattice calculations (Section IV), whereas in the most precise and involved gauge modeling, we have presented only the outcome from lattice simulations (Section V).
The analytical results of the global modeling estimate correctly the right order of magnitude of the Higgs decay time scale. When studying such scenario in the lattice, including all non-linearities within such scheme, we find that the actual Higgs decay takes longer, typically a factor z e /z eff ∼ 3.17q 0.44 larger. This is because the analytical calculations are only capable of estimating the order of magnitude of the time scale when sufficient energy has been transferred into the extra scalar fields mimicking the EW gauge bosons. However, that time only signals the moment when the Higgs condensate really starts noticing that it is coupled to extra species. From then onwards, the Higgs energy density begins to decrease in a noticeable manner, being transferred to the most strongly coupled species, the EW gauge bosons. It is this decrease of the energy of the Higgs that should be interpreted as the decay of the Higgs. Eventually, the Higgs energy density saturates to an approximately constant value. Around the same time, the energy of the species coupled to the Higgs has also stopped growing, and saturates into slowly evolving magnitudes. Very interestingly, the same pattern and time scales are observed in the gauge model, but the final fractions of energies are different. The time scale z e (q) that characterizes the end of the Higgs decay is given by Eq. (80), which represents a factor z e /z eff ∼ 3.68q 0.42 longer with respect to the analytical predictions. We see therefore that, at the end, the differences between the global and gauge modelings are not so relevant, at least in terms of the estimation of the z e (q) time scale.
Interestingly, at the time z ≈ z e (q), both in the global and gauge scenarios, we see that the distributions of fields reach equipartition. In the global model we find that the kinetic energy of the Higgs becomes equal to the sum of the gradient energy of the Higgs plus the interaction with the χ i fields, E ϕ K E ϕ G + E int . This equality holds to better than 1 % from z z e onwards. In the gauge scenario, we find that at the onset of equipartition, also the kinetic energy of the Higgs becomes equal to the sum of the covariant gradient energy (which includes the Higgs-gauge bosons interactions) plus the Higgs potential, E K E GD + E V . This equality also holds to better than 1% from z z e (q) onwards. At the same time, the electric and magnetic energy densities also reach equipartition to better than 1%, E E E M . The distribution of energy in the gauge scenario is actually universal, since the system always reaches equipartition, with E K E GD representing 30% of the total energy, and E E ≈ E M representing 20% each. Let us note that in both global and gauge scenarios, once in the stationary equipartitioned regime, the potential energy becomes gradually more and more irrelevant.
Before we conclude, let us note that the postinflationary decay of the Higgs analyzed here, is very similar to the analogous decay during reheating after Higgs-Inflation [3,[43][44][45][46]. The contexts are however very different. In Higgs-Inflation the Higgs plays the role of the inflaton and dominates the energy budget of the universe, so the decay of the Higgs after inflation truly represents the actual reheating of the universe. In the case we have studied in this paper, the Higgs is simply a spectator field during inflation, and its energy density is only a marginal fraction of the inflationary one, see Eq. (20). In Higgs-Inflation, a non-minimal coupling ξϕ 2 R to gravity is required, with ξ ∼ O(10 4 ) √ λ. The resonance in both Higgs-Inflation and Higgs-Spectator scenarios is dominated by the decay into the gauge bosons W ± , Z. The resonance parameter, however, scales as q ∼ g 2 λ ξ in Higgs-Inflation, versus q ∼ g 2 λ in our Higgs spectator scenario. Therefore, the resonance is ∼ 10 3 √ λ 001 times broader in Higgs-Inflation than in the Higgs-spectator case. However, in Higgs-Inflation the non-perturbatively produced gauge bosons (at each Higgs zero-crossing), decay very fast into the SM fermions via perturbative decays. So for around ∼ 100 oscillations of the Higgs, the resonance is blocked in Higgs-Inflation, simply because the occupation numbers of the gauge bosons do not pile up [44]. This phenomenon is called combined-preheating, and it is absent (or at most is only a marginal effect) in the Higgs spectator scenario.
To conclude, let us note that our paper is intended to be the first one of a series, where we plan to analyze further the details of the Higgs decay (i) and its cosmological consequences (ii). In particular, (i) The results obtained here have gone far beyond the analytical ones available in the literature [19,30]. We have presented different approaches to the nonperturbative and non-linear dynamics of the decay process. Our most precise results are the outcome from our simulations in Section V, corresponding to an Abeliangauge model mimicking the structure of the SM interactions. Even though there is a good motivation to neglect (on the initial stages of the process) the truly non-Abelian nature of the interactions [30], only lattice simulations which fully incorporate the non-Abelian SU (2) × U (1) structure of the SM, will really tell about the (un)importance of the corrections due to the nonlinearities in the gauge sector (absent in the abelian approach). Besides, the details of the stationary stage might very well (and indeed most likely will) change when the full non-Abelian structure of the SM is restored. Therefore, even if the time scales of the start of the Higgs decay and onset of stationary regime may (expectedly) not change much, the fine details can only be quantified in light of such non-Abelian simulations, which are currently beyond our present work. Moreover, in order to assess with even a higher degree of realism the final outcome of the energy distribution among fields, thermal corrections [47][48][49] and fermions [50][51][52] should be effectively incorporated into such simulations.
(ii) The post-inflationary decay of the SM Higgs may have several observable consequences. It has been recently proposed [53] the possibility of realizing baryogenesis via leptogenesis, thanks to the Higgs oscillatory behaviour. The time dependence of the Higgs condensate oscillations can create an effective chemical potential for the lepton number, which could lead to the generation of a lepton asymmetry in the presence of righthanded majorana fermions with sufficiently large masses. The electroweak sphalerons would then redistribute such asymmetry among leptons and baryons. Secondly, the fields excited from the decay of the Higgs may act as a source of gravitational waves [54][55][56][57][58][59]. The case of the charged fermions of the SM was considered [35], but it is expected that the background of gravitational waves from the EW gauge bosons contributes to a much larger signal [35]. Besides, the fact that the Higgs is a condensate varying at superhorizon scales, may give rise to interesting anisotropic effects [41,42] in the amplitude of such background of gravitational waves. Thirdly, it is indeed possible that the gauge field production that we have described in this paper could provide the necessary conditions for primordial magnetogenesis. Although it is likely challenging to obtain a sufficiently large correlation length, it is conceivable that an inverse cascade process provides the appropriate mechanism for the growth of an initially small correlation length [60,61]. Finally, if dark matter is a gauge singlet field coupled to the Higgs, it is also possible that the Higgs oscillations could produce the right amount of dark matter, such that its distribution could account for the correct relic abundance [62].
where we have defined Jn at the lattice pointn as and the lattice covariant derivatives as D + µ φ = 1 dµ (U µ φ(n+μ)−φ) and D − µ φ = 1 dµ (φ−U * µ (n−μ)φ(n−μ)). One needs to check that for all times, the discrete Gauss law (A12) is conserved. In particular, we require that for all times ∆ G 1, where We have checked that this is indeed the case. In particular, we find that depending on the simulation, at the end of the running time the Gauss law is in fact only marginally broken, with ∆ G 10 −13 − 10 −15 .
All results presented in this work have been obtained for N = 128 points for both global and Abelian-Higgs simulations. Apart from N , we also need to fit the range of momenta that we want to cover in our simulations. This is a crucial step, as this range must be chosen carefully in order to capture all the relevant phenomenology of the Higgs decay. Let us call p min the minimum momentum covered by the lattice. We can then fix the length of the cube L and the maximum momentum covered by the lattice, p max , in terms of N and p min as Note that in this appendix, k refers to physical momentum and p to lattice momentum. As discussed in section III, the Higgs eom possess a well known structure of resonance bands, which can be either of the form 0 < k < k * or k min < k < k * . We expect these momenta to be physically excited, at least at the first stages of the Higgs decay. Therefore, we must have a good coverage of this range of momenta. Let us define the coefficient α c The larger α c is, the better the infrared coverage of the resonance band, but the worse the ultraviolet scales are captured. In order to probe well the posterior displacement of the spectra to higher momenta when the system becomes non-linear, we need to choose α c judiciously. With this idea in mind, we have determined for each q, the α c parameter that ensures a good infrared coverage without spoiling the ultraviolet part of the spectra. For each simulation, we have fixed α c , typically within the range 4 α c 11, depending on N and the appropriate dynamical range. We show in Fig. 17 two particular spectra obtained from simulations of the global scenarios at two different times, for different (N, α c ) parameters. Apart from a better or worse coverage of the ultraviolet or infrared regimes, the main physical results are well captured in all simulations, and are also consistent between each other. The same consistency is observed in the Abelian-Higgs simulations, making our results robust versus lattice artefacts.
state of a scalar field in a FLRW universe where we have |f k | 2 = 1 2a 2 osc ω k,osc . (B3) Here, ω k,osc ≡ k 2 + a 2 osc m 2 osc is the frequency of the field at the time z osc , and m osc is the mass at this same time, m 2 osc = (∂ 2 V /∂f 2 )(z osc ) with V the potential. The mode f k also contains an arbitrary random constant phase ∀Arg(f k ). To mantain isotropy properties, we add both left-moving and right-moving waves, so that we take where θ 1 and θ 2 are constants with θ i ∈ [0, 2π). In the discrete lattice, we set the fluctuations in momentum space so that, from lattice point to lattice point, |f k | varies accordingly to Eq. (B2), and the phases θ 1 and θ 2 vary randomly within the interval θ i ∈ [0, 2π). From the properties of the Lamé equation discussed in Section III, we know that depending on the value of the resonance parameter q ≡ g 2 /(4λ), the Higgs equation has a certain structure of resonance bands. As we see in Fig. 4, the most infrared band is always the one with the greatest Floquet index, and we hence expect that the Higgs decay will be dominated by this band, at least at initial times. It has a maximum at a given momentum, which we call k max . Therefore, this allows us to set a cutoff to the probability spectrum (B2), such that for k > k max , |f k | = 0. As it should be, we have confirmed that changing this cutoff within a wide range of values does not significantly modify our results.

Abelian-Higgs model: Gauss conservation law
We now discuss how we set the initial quantum fluctuations in the Abelian-Higgs model. Caution must be taken in this step, because as we will see, we must ensure that the Gauss condition (Eq. (A13) in the discrete) holds at the beginning of the simulations.
Let us come back to natural variables. In this section, we define h j ( z, z osc ) with j = 1, 2 to be the fluctuations of the two components of the conformally rescaled Higgs field at the time z osc . Let us also define h j (k) ≡ h j ( k, z osc ) to be their corresponding Fourier transforms. Following very closely our discussion of the initial conditions in the global modelling, we impose, over the two components of the Higgs, the spectra Here, A and B are quantities that change, from point to point of the lattice in momentum space, according to the probability distribution function P (y)dy = 2y y 2 e − y 2 y 2 dy (B7) where y 2 ≡ (2ω j ) −1 . The frequency is ω j ≡ | k| 2 + m 2 j , with the (rescaled) masses defined as m 2 1 ≡ (3h 2 1osc + h 2 2osc )β 2 = 3h 2 1osc β 2 , m 2 2 ≡ (h 2 1osc + 3h 2 2osc )β 2 = h 2 1osc β 2 , where h 1osc ≡ h 1 (z osc ) and h 2osc ≡ h 2 (z osc ). For the last equality, we have used that, for the initial conditions of the Higgs homogeneous mode given in Eq. (74), we have h 2osc = 0. From (B6), the fluctuations of the Higgs derivatives are Also, the four different phases vary, in momentum space, from lattice point to lattice point. These phases would vary in principle randomly within the interval θ i ∈ [0, 2π), but as we are working also with gauge bosons, we need to preserve the Gauss law initially. Due to this, we thus may need to impose one simple constraint to the phases.
Let us discuss this in more detail. As mentioned before, we must ensure at the initial time the Gauss law (70), Here, j 0 (z) is j 0 (z) ≡ qβ 2 Im[(h 1 −ih 2 )(h 1 +ih 2 )]. Therefore, the quantum fluctuations we impose to the gauge fields at z osc must preserve this condition. Let us write the Gauss law (B10) in momentum space as where j 0 (k) is the Fourier transform of j 0 (z) at the time z = z osc , and p min is the minimum momentum of the lattice. This allows us to set fluctuations to the gauge fields in the following way: First, for a given lattice point in momentum space, we produce the Higgs fluctuations according to Eqs. (B6) and (B9). With these Higgs fluctautions, we obtain the correspondent fluctuations of j 0 (z) and its corresponding Fourier transform j 0 (k). Finally, we fix V i ( k, z osc ) according to Eq. (B11). We have then obtained a spectrum of initial gauge fluctuations. However, in order this procedure to be valid, we must ensure that our current j 0 (k) does not possess a zero mode, i.e. j 0 ( k = 0) = 0. This requirement can be clearly seen in Eq. (B11). This is equivalent to saying that there must not be a total electric charge in our lattice box. However, from the spectrum of Higgs fluctuations described above, we obtain j( k = 0) = d 3 zj 0 (z) = (B12) with Re[h 1 (k)h 2 (k) − h 2 (k)h 1 (k)] = e 2 λ β 2 × cos θ3+θ4−θ1−θ2 (B13) This quantity is not zero in general. There does not seem to be a particular reason why we should have a total electric charge in our box, so we should find a way of making Eq. (B13) null. We have found two different ways of modifying slightly the initial quantum fluctuations of the Higgs field to make the integrand of Eq. (B13) is zero, which do not modify significantly the amplitude of the fluctuations with respect to the approach used in the global model. The first one is to impose, at each lattice point, the following constraint to the four arbitrary phases of the Higgs fluctuations θ 4 = θ 1 + θ 2 − θ 3 + π (B14) so that the phases θ 1 , θ 2 and θ 3 are randomly generated within the interval θ i ∈ [0, 2π), and θ 4 is fixed through Eq. (B14). The second one is to leave the four phases totally random, but to do, at each lattice point, the following shift to the Higgs fluctuationṡ is a sum over all lattice points. One can easily confirm that both methods make zero the integrand of Eq. (B13). In practice, we have confirmed that both methods produce almost identical results. This is normal, as in order to trust our lattice simulations, the way in which we set the initial fluctuations must not play any relevant role, as long as their amplitude does not significantly change.