Gravitational wave production from the decay of the Standard Model Higgs field after inflation

During or towards the end of inflation, the Standard Model (SM) Higgs forms a condensate with a large amplitude. Following inflation, the condensate oscillates, decaying non-perturbatively into the rest of the SM species. The resulting out-of-equilibrium dynamics converts a fraction of the energy available into gravitational waves (GW). We study this process using classical lattice simulations in an expanding box, following the energetically dominant electroweak gauge bosons $W^\pm$ and $Z$. We characterize the GW spectrum as a function of the running couplings, Higgs initial amplitude, and post-inflationary expansion rate. As long as the SM is decoupled from the inflationary sector, the generation of this background is universally expected, independently of the nature of inflation. Our study demonstrates the efficiency of GW emission by gauge fields undergoing parametric resonance. The initial energy of the Higgs condensate represents however, only a tiny fraction of the inflationary energy. Consequently, the resulting background is very suppressed, with an amplitude $h^2 \Omega_{\rm GW}^{(o)} \lesssim 10^{-29}$ today. The amplitude can be boosted to $h^2 \Omega_{\rm GW}^{(o)} \lesssim 10^{-16}$, if following inflation the universe undergoes a kination-domination stage; however the background is shifted in this case to high frequencies $f_p \lesssim 10^{11} {\rm Hz}$. In all cases the signal is out of the range of current or planned GW detectors. This background will therefore remain, most likely, as a curiosity of the SM.


I. INTRODUCTION
Gravitational waves (GW) are ripples of the spacetime which propagate at the speed of light. Until very recently, the existence of GW had only been proven indirectly, through the modulation of the orbital period of binary pulsars [1]. Advanced LIGO [2] has just announced, however, the first direct detection of GW from the coalescence of two massive black holes. This historical event opens a new window into the Universe, which will allow us to probe astrophysical and cosmological environments previously inaccessible. This milestone detection will very likely inaugurate a new era in cosmology.
A compelling possibility is the Higgs-Inflation scenario [48,49], in which the SM Higgs is responsible for inflation, thanks to a large non-minimal coupling to gravity. Interestingly, even if inflation is driven by an inflaton scalar field other than the SM Higgs (typically a singlet under the SM), the Higgs may still play a very relevant role after inflation. In this work we assume that the SM Higgs is not responsible for inflation, and indeed we consider it very weakly coupled to the inflationary sector, even at loop order. Under this circumstance, we can distinguish two possibilities: i) If the Higgs is minimally coupled to gravity, it behaves as a spectator field during inflation. It then forms a condensate with a large vacuum expectation value (VEV), and a correlation length exponentially larger than the Hubble radius [50,51]. If on the contrary, ii) the Higgs is non-minimally coupled to gravity with a sufficiently large coupling, the Higgs is not excited during inflation [29,52], but it is however strongly excited during the transition period at the end of inflation [53]. In this case, the Higgs forms also a condensate with large VEV, but with a correlation only of the order of the Hubble radius.
In either case i) or ii), shortly after the end of inflation, the Higgs condensate starts oscillating around the minimum of its potential. This gives rise to its decay into all the species of the SM, as the latter are non-perturbatively excited through parametric effects [30][31][32][33][34]54] (see also [55][56][57] in the Higgs-Inflation context). All the SM species coupled directly to the Higgs, i.e. the electroweak gauge bosons W ± , Z, and the massive fermions (quarks and charged leptons), are all highly excited. This is a violent non-equilibrium process, creating large timedependent matter density inhomogeneities, which therefore act as a classical source of GW. arXiv:1602.03085v2 [astro-ph.CO] 31 May 2016 The decay of the Higgs into fermions, and the corresponding GW production, was studied in [31], following the formalism of [58,59]. Fermions are excited through parametric effects [60,61], though the growth of their occupation numbers is Pauli blocked. The most energetic fermion species excited is the top quark, since its Yukawa coupling is the largest one within the SM. In this paper we focus instead in the production of GW by the gauge bosons. The gauge field production is expected indeed to be more efficient than that of fermions, as their occupation numbers grow exponentially [30,[32][33][34]. Most of the energy of the Higgs condensate is actually transferred into the electroweak W ± , Z gauge bosons. Therefore, even if the final GW background is contributed by all the Higgs decay product species, the gauge fields we study here represent in fact the dominant contributors.
The decay of the SM Higgs condensate into gauge bosons after inflation, has been studied recently very extensively. It was first studied in [30,32] with analytical methods based on the linear regime, and later extended in [33]. Beyond the linear regime, a full set of lattice simulations of the process was presented in [33], modeling the SM gauge interactions with an Abelian-Higgs set-up. Although this is just an approximation to the gauge structure of the electroweak interactions, the non-Abelian effects can be arguably neglected for a large fraction of the physically motivated values of the Higgs self-coupling. The outcome of these simulations, though neglecting the truly non-Abelian structure, represent a precise calculation of the dynamics of the SM after inflation, fully incorporating the nonlinear and non-perturbative effects of the SM gauge interactions between the Higgs and the W ± and Z gauge bosons.
More recently, lattice simulations presented in [34] have considered the non-Abelian structure of the SM. They have shown interesting effects due to the new nonlinearities introduced. The non-Abelian corrections are however suppressed by the smallness of the Higgs selfcoupling [32,33]. In high-energy inflationary models, the Higgs self-coupling runs in fact into small values [62,63], making the non-Abelian corrections less relevant, the larger the energy scale. In this paper we are mostly interested in scenarios with the highest possible energy scale of inflation, as this enhances the production of GW in the system. Therefore, the use of an Abelian modeling will suffice for our aim to study the GW production from the SM fields after inflation.
The structure of this paper is as follows. In Section II we review the creation of a Higgs condensate during or towards the end of inflation. In Section III we review the post-inflationary dynamics of the Higgs and of its decay products, summarizing the results of [33]. In Section IV we discuss our formalism to study GW production in this process. In Section V we present our results, describing the general features of the GW spectra obtained from our lattice simulations. In Section VI we parametrize the GW spectra as a function of the Higgs initial amplitude, Higgs self-coupling, and post-inflationary expansion rate. In section VII we discuss how the GW background redshifts until today. Finally, in Section VIII we wrap up our results and conclude.
From now on, m 2 p = (8πG) −1 = 2.44 · 10 18 GeV is the reduced Planck mass, and we consider a flat Friedman-Robertson-Walker background metric ds 2 = a 2 (t)(−dt 2 + d x 2 ), with t the conformal time and a(t) the scale factor.

II. HIGGS EXCITATION DURING (OR TOWARDS THE END OF) INFLATION
Relevant properties of the SM Higgs, like the nature of its gravitational coupling, or its coupling to the inflationary sector, are currently unknown. As a consequence, the role played by the Higgs, and in general its dynamics during the early Universe, are uncertain. In this paper we consider the Higgs to be sufficiently weakly coupled to the inflationary sector, so that the Higgs does not develop a super-Hubble mass during inflation. The need to reheat the Universe after inflation requires, however, the presence of a Higgs-inflaton coupling, induced through radiative corrections from some mediator field(s).
The only scale-free renormalizable Higgs-inflaton operator is g 2 φ 2 Φ † Φ, where Φ is the SM Higgs doublet, φ is the inflaton, and g 2 is the coupling strength. In order not to spoil the inflationary predictions, it is required that g 2 g 2 * = 10 −6 [42]. Furthermore, to avoid that the Higgs developing a super-Hubble mass during inflation, one needs g 2 ≤ 10 −4 g 2 * . We will consider this second inequality as valid from now on, providing in this way an operational definition of what we mean by the Higgs being sufficiently weakly coupled to the inflationary sector.
The concrete particle physics realization of inflation has eluded any clear identification so far. The inflationary dynamics is normally described in terms of a scalar field, the inflaton, a singlet under the SM, and with a vacuum-like energy density. For simplicity, we will describe inflation as a de Sitter background with Hubble rate H * MEW, where MEW ∼ O(10 2 ) GeV is the electroweak (EW) scale. The current upper bound of the inflationary Hubble rate is [64] so, in principle, there is 'plenty of room' to fulfill the demand MEW H * ≤ H (max) * . In the unitary gauge, the Higgs doublet can be written as Φ = ϕ/ √ 2, with ϕ a real degree of freedom with renormalized potential at large field values (ϕ MEW) where the effective self-coupling λ(ϕ), encapsulates the radiative corrections to the Higgs potential [65,66]. Under the circumstance that the Higgs is not responsible for Inflation, and is decoupled from (or weakly coupled to) the inflationary sector, we can consider two possibilities: (i) Higgs minimally coupled to gravity -. In this case, the Higgs plays no dynamical role during inflation. It behaves as a light spectator field, independently of its initial amplitude [30,67]. The Higgs performs a random walk at superhorizon scales during inflation, reaching within few e-folds an equilibrium distribution with variance [50] The Higgs forms this way a condensate with a large VEV during inflation. A typical amplitude of the Higgs condensate is then ϕrms ∼ H * , (almost) independently of the Higgs self-coupling for reasonable values of λ. The scale over which the Higgs condensate amplitude fluctuates, i.e., the correlation length of the Higgs condensate, is exponentially larger than the Hubble radius [50], l * ∼ The Higgs condensate is then homogeneous within cosmological scales.
(ii) Higgs non-minimally coupled to gravity -. An interaction ξΦ † ΦR, with R the Ricci scalar, is required by the renormalisation of the SM in curved space.
If the value of ξ at the inflationary scale lies below ξ 0.1, the Higgs is sufficiently light during inflation, and hence we recover scenario i). If on the contrary, ξ > 0.1, the Higgs becomes too heavy during inflation. It develops a super-Hubble mass and, consequently, it is not excited as a condensate. However, the sudden drop of the curvature R, during the transition from the end of inflation to a standard power-law post-inflationary regime, induces a non-adiabatic change in the effective Higgs mass m 2 ϕ = ξR. This translates into a significant excitation of the Higgs modes at the Hubble scale k ∼ a * H * . Following the recent analysis of [53], the Higgs excitation acquires a large variance 1 depending on the post-inflationary equation of state w. Furthermore, if after inflation the inflaton oscillates around the minimum of its potential, the periodic time-dependent behavior of R excites the Higgs modes through parametric resonance. The Higgs amplitude cannot in any case exceed [53] as the Higgs self-interactions prevents any further growth above this value. In summary, in the presence of a large non-minimal coupling, the Higgs forms a condensate immediately after inflation, with 1 Denoting the post-inflationary equation of state as w, in [53] it is found that a large amplitude bounded as ϕrms H * /λ 1/2 ξ 1/4 . The coherent scale is, however, only of the size of the Hubble radius, l * H −1 * , as the fluctuations at different horizon patches are uncorrelated. The running of the Higgs self-coupling λ(µ) has been computed up to three loops [62,63]. The self-coupling decreases with energy (dλ/dµ < 0), becoming negative at a given critical scale µc. Due to this, the Higgs potential Eq. (2) possesses a maximum (a barrier) at a scale µ+ µc, crosses zero at µc, and it (possibly) develops a negative minimum at higher energies. These scales are very sensitive to the Higgs mass mH , the strong coupling constant αs, and the top quark mass mt. For the SM central values for αs and mH , as well as the world average top quark mass mt = 173.34GeV [68], one finds µ+ ≈ 7×10 9 GeV and µc ≈ 10 10 GeV. However, considering a value of mt two/three sigma below its central value, we put the critical scales at µ+, µ0 ≥ 5 × 10 16 GeV H (max) * , which is one way of ensuring the stability. Another way is to consider Higgs portals to scalar fields, changing the running of λ so that it always remains positive [69]. In this work we assume that the Higgs potential is stable up to inflationary energies, so we consider that λ never becomes negative. The higher the energy scale of inflation, the smaller is λ, with reasonable values only within the interval 10 −2 λ < 10 −5 [33] (λ ∼ 10 −5 being only marginally valid).

III. POST-INFLATIONARY DYNAMICS OF THE SM HIGGS AND ITS DECAY PRODUCTS
As we have just discussed, either during inflation [case i)] or just towards its end [case ii)], the Higgs is excited in the form of a condensate with a large VEV. As a consequence, the Higgs condensate starts oscillating around the minimum of its potential, soon after the end of inflation. Each time the Higgs crosses zero, all particles coupled to the Higgs, i.e. the electroweak gauge bosons and the charged fermions of the SM, are created in nonperturbative bursts [30][31][32][33][34]. In the case of gauge bosons this phenomenon is called parametric resonance, and it is very similar to the decay process of an inflaton in preheating with a quartic potential [70]. The main difference with respect to preheating is that, in the present scenario, the Higgs does not dominate the energy budget of the Universe. On the contrary, given the typical Higgs condensate amplitudes ϕrms [Eq. (3) in case i), or Eq. (4) in case ii)], the Higgs energy density is always much smaller than the energy density of the inflationary sector ρ * ≡ 3m 2 The decay of the Higgs and the dynamics of its energetically dominant decay products have been recently studied under these circumstances, with the help of classical lattice simulations [33,34]. The lattice approach takes into account the non-linearities of the system beyond the analytical approach (based on the initially valid linear regime). Strictly speaking, however, the analysis presented in [33,34] only describes the post-inflationary dynamics of the system in the spectator field case i), when the Higgs is excited as a condensate during inflation.
The post-inflationary dynamics in the case ii), when a non-minimal coupling to gravity ξϕ 2 R is present, can however be expected to be (qualitatively) similar to that of case i), at least for moderate values of ξ. The reason for this is that the coupling to the scalar curvature R, introduces an additive term in the Higgs oscillation frequency, proportional to the square of the Hubble rate, R ∝ H 2 . As H 2 ∼ 1/t 2 , the new term decays rapidly in time. Therefore, the Higgs oscillations in the case ii) are expected to be initially modulated by the non-minimal coupling to gravity, but tending rapidly to the oscillations of case i) (characterized by an oscillation frequency ∝ λϕ 2 , given by the Higgs self-interactions).
In what follows, we will describe the details of the postinflationary dynamics, considering only the scenario i) with ξ set to zero. This should capture equally well the dynamics when 0 < ξ 0.1. For the scenario ii) with a non-minimal coupling ξ > 0.1, we expect the dynamics to tend rapidly to the case of ξ < 0.1, unless ξ is extremely large. When ξ 1, a new analysis beyond our current study should be performed. However, the larger the ξ, the smaller the GW production is expected to be 2 . Therefore, we only focus from now on in the details of scenario i) with ξ = 0, which maximizes the GW production. This should still capture qualitatively well the dynamics of the system, even in the presence of a moderate non-minimal coupling to gravity, as long as ξ is not extremely large.

A. Abelian model of the electroweak interactions
In order to describe the production of gauge bosons from the Higgs decay after inflation, we will follow the approach presented in [33]. We will model the electroweak 2 The larger the ξ, the faster the energy transfer from the Higgs into its decay products, as the particle production (due to the Higgs oscillations) is expected to be modulated by a larger oscillation frequency. As we will explain later on in Section IV, the faster the decay of the Higgs proceeds, the higher the frequency of the GW produced during the process. However, the amount of GW produced is mostly determined by the energy stored initially in the Higgs condensate, and not by the rapidity of its decay. In the scenario ii), the larger the ξ, the smaller the initial energy stored in the Higgs [53]. Therefore, the larger the ξ, the larger the frequency of the GW background, but the smaller the amplitude of the background. We expect therefore the GW production in scenario i) with ξ = 0 to be the largest possible one.
sector of the SM with an Abelian-Higgs set-up, ignoring the non-linearities arising from the full non-Abelian structure of the electroweak interactions. Non-abelian corrections over the Abelian dynamics are indeed suppressed as ∝ 1/ √ q [32][33][34], where q ≡ e 2 /λ is the resonance parameter of the gauge field(s), and e 2 represents the Abelian coupling (mimicking either of the W ± or Z gauge couplings). The Abelian approximation to the full Higgs-gauge electroweak interactions works better the larger is the resonance parameter q. For GW production through the Higgs decay products, we are interested in inflationary scenarios with the highest possible energy scales, since this enhances the GW produced, see Section IV. If the inflationary Hubble rate H * is sufficiently high, say of the order of (though somewhat smaller than) its current upper bound H * H (max) * ∼ 10 14 GeV, the value of λ runs into small values. This implies that the resonance parameters are, in fact, rather large, 1. Therefore, taking the Abelian approximation in this regime is well justified 3 , and hence in the present work we will consider the Abelian modeling from now on.
Within the Abelian-Higgs modeling, the interactions of the Higgs with a gauge boson species A µ can be described by the action Here Dµ ≡ ∂µ − iAµ is the covariant gauge derivative, Fµν ≡ ∂µAν − ∂ν Aµ is the (Abelian) field strength, and e is the strength of the Abelian coupling. In the Abelian modeling the Higgs needs to be considered as a complex scalar field Φ ≡ 1 As we are dealing with a gauge theory, we have a gauge freedom to choose the field components. This allows us to set the gauge A0 = 0 from now on. By varying the action, the equations of motion can be derived as Eqs. (7) and (8) of the system describe the dynamics of the Higgs and the gauge boson, while Eq. (9) is the Gauss law, representing a constraint that must be obeyed at all times. The gauge-invariant electric and magnetic fields are Ei ≡ F 0i and Bi = 1 2 ijk F jk .
We model the scale factor as where w is the equation of the state of the universe after inflation. The value of w depends on the inflationary sector, which is not specified here, so we will consider various values representing different expansion rates. For example, for matter-domination (MD), radiation-domination (RD) or kination-domination (KD), the equation of state is w = 0, 1 3 , 1, respectively. From now on we fix a * = a(t * ) ≡ 1.
It is convenient to define dimensionless spacetime vari- It is also convenient to define dimensionless Higgs and gauge field amplitudes as | is the initial modulus of the Higgs field at the end of inflation. To distinguish between variables, we use a dot or a prime to denote differentiation with respect to conformal or natural variables,˙≡ d/dt or , ≡ d/dz. From now on, all spatial derivatives will also be taken with respect the new variables, unless otherwise stated. We also define a dimensionless covariant derivative as Di ≡ ∂ ∂z i − iVi. With these changes, Eqs. (7)-(9) can be written as where q ≡ e 2 λ is the resonance parameter, and we have defined the parameter characterizing the initial Higgs amplitude. Let us denote by W + µ , W − µ and Zµ, the Abelian version of the SM electroweak W, Z gauge bosons. Let us also denote by gW and gZ the SM gauge couplings of the Higgs to such bosons. In order to mimic correctly the real Higgs-gauge interactions, we need to identify e 2 = g 2 /4 in Eqs. (7)-(9), with g 2 = g 2 W for W bosons and g 2 = g 2 Z for Z bosons. This translates into the following resonance parameters in Eqs. (13)-(15), .
In principle, for each of the three gauge bosons W ± and Z, there should be an equation of motion like Eq. (8) [or equivalently Eq. (14)], and a Gauss constraint like Eq. (9) [or equivalently Eq. (15)]. However, we demonstrated in [33] that this system of three gauge fields can be mapped identically into another system with a single gauge boson, defined as with associated resonance parameter (gauge coupling) See Section V.B. in [33] for details. At any moment one can recover the original fields by simply taking At very high energies, the running of the SM gauge couplings show that g 2 The energy density of the Higgs + gauge fields is where V * is the value of the Higgs potential at the end of inflation, and the function E s (z) is formed by the sum of the following energetic contributions: Here, E K and E V are the kinetic and potential energies of the Higgs field, and E E and E M are the electric and magnetic energy densities of the super gauge-boson S µ of Eq. (18) (hence containing the contribution from all the gauge bosons W ± , Z), where we have defined dimensionless electric and magnetic fields as Ei = Ei/H 2 * and Bi = Bi/H 2 * . The last contribution EGD is a gauge-invariant term formed by the product of covariant derivatives of the Higgs field, hence representing the energy stored in the spatial Higgs gradients and the Higgs-gauge interactions. We have plotted in Fig. 1 the different energy contributions as a function of time for the resonance parameter qs = 79. We have normalized each energy component as Ei/Es, and we have removed their oscillations, hence showing only the corresponding envelope functions. We see that initially the dominant contributions come from the kinetic and potential energies of the Higgs field. This corresponds to the oscillations of the Higgs condensate around the minimum of its potential, before the gauge fields backreact onto the Higgs. Meanwhile, the other components of the energy (EE, EM and EGD) grow very fast, due to the energy transfer -via parametric resonance -from the Higgs to the gauge fields. We can understand the time evolution of these energies in light of the context of the next subsection.

B. Summary of the post-inflationary dynamics of the Higgs and its decay products
The lattice simulations presented in [33] represent a precise calculation of the dynamics of the SM after inflation, fully incorporating the nonlinear and nonperturbative effects of the SM gauge interactions between the Higgs and the W ± and Z gauge bosons, though neglecting the truly non-Abelian structure. The advantage of a lattice approach with respect an analytical one is that the first includes the effect of the non-linearities of the system into the Higgs decay dynamics, which become relevant very soon after the end of inflation. Here we just briefly summarize the results obtained in [33].
The production of the gauge bosons from the postinflationary Higgs decay is controlled by three parameters: the amplitude of the Higgs condensate ϕ * at the end of inflation, the equation of state w of the Universe in the period following immediately after inflation, and the resonance parameter q ≡ g 2 /4λ of the gauge fields, with g 2 equal to either of the W ± , Z gauge couplings 4 . Taking the end of inflation as the initial time z * = 0, the dynamics of the system at later times z ≥ z * can be characterized by three time scales: • z = zosc: This time signals the moment after inflation when the Higgs effective mass becomes larger than the Hubble rate. This 'forces' the amplitude of the Higgs condensate to start rolling down its potential. Previous to this moment, during z * ≤ z ≤ zosc, the Higgs amplitude remains frozen in slow-roll 5 . Hence z = zosc signals the onset of the Higgs oscillations around the minimum of its potential. Every time the Higgs crosses around zero, particle creation via parametric effects occur. It is found that zosc O(10), the exact value depending on the particular Higgs amplitude ϕ * and postinflationary equation of state w. Therefore, the Higgs oscillations start, in all cases, shortly after inflation ends. From then on, top quarks and gauge bosons are strongly created every time the Higgs crosses through zero. After few oscillations, the gauge boson production dominates, as the gauge fields develop parametric resonance, whereas the fermions are Pauli blocked [72]. The excitation of the dominant species Z, W ± is very similar to the excitation of preheat fields coupled to the inflaton, as preheating is due to oscillations of an inflaton with quartic potential.
• z = zi: This second time scale signals the moment when the produced gauge bosons have accumulated sufficient energy such that they start backreacting onto the Higgs condensate, starting to affect the latter severely. At z ≥ z i there is a sharp decrease of the amplitude and energy density of the Higgs condensate. This scale depends strongly on the resonance structure of the gauge field dynamics, characterized by the particular value of the resonance parameter qs ≡ (g 2 For the larger resonance parameters (due to small values of λ), zi tends to be shorter, while weakerresonance parameters cases have larger zi's. In practice, zi zosc + (O(0.1) − O(10 3 )). 4 The production of fermions is equally controlled by the same parameters, ϕ * , w, q, but substituting the gauge coupling g 2 in the definition of q by the Yukawa coupling y 2 i , with the i-index referring to the fermion species [31]. 5 In the scenario ii) with ξ > 0.1, the Higgs starts rolling down its potential immediately after inflation, since the effective mass of the Higgs (given by its non-minimal coupling to gravity) is larger than the Hubble already at z * . Hence, zosc = z * = 0 in this case.
• z = ze: This third and final time scale signals the end of the transfer of energy from the Higgs into the SM species, as well as the onset of a stationary regime. The time ze depends on the resonance parameter as a power law ∝ q α s , where a phenomenological fit to the simulations shows α ∼ 0.42 [33]. Therefore, the greater the qs, the longer the time ze. At z ≥ ze, the energies of the different components of the system have established an equipartition regime, with fixed relative ratios of the energies independent of qs. For reasonable values of λ, q s ranges between ∼ O(10) and ∼ O(10 3 ), so ze ranges between ∼ 10 3 and ∼ 10 4 .

IV. GRAVITATIONAL WAVE PRODUCTION
Gravitational waves (GW) are tensor perturbations which propagate following the equation of motion [22] where the source of GW, Π TT ij , is the transverse-traceless (TT) part of the anisotropic stress tensor Πij. In our case, in the presence of both scalar and vector fields, the source is effectively given by [73] where {...} TT represents the TT part of the quantity inside the brackets, Dµ is the gauge covariant derivative in Eq. (6), Ei ≡ F 0i and Bi = 1 2 ijk F jk are the electric and magnetic fields formed from the superfield S µ in Eq. (18), and we have discarded a term ∝ (|E| 2 + |B| 2 )δij because it is actually a pure trace term. Note that here the charge must be identified with that of S µ , so that It is convenient to redefine the tensor mode amplitude through a conformal redefinition like (recall that initially we take a * = 1) so that Eq. (27) can be written in terms of the dimensionless variables Eqs. (11), (12) as with We recall that E i ≡ E i /H 2 * , B i ≡ B i /H 2 * , and qs = (2g 2 W + g 2 Z )/4 is the total resonance parameter Eq. (19). It is clear from Eq. (32) that both the Higgs field and the gauge bosons will contribute as a source of GW.
The spectrum of the GW energy density contained within a volume V , and normalized to the total energy density ρ tot of the Universe (at the time of GW production), can be written in the continuum as where ... 4π ≡ 1 4π dΩ k ..., with dΩ k a solid angle differential in k-space.
In light of the parameters factorized out in the source term of Eq. (30), it is convenient to define a new variablē It is then useful to express Ω GW (k, z) in terms of the natural variables of the problem zµ and wij. We can factorize out this way the dependence with the Hubble scale H * and the background expansion rate as where and In order to derive Eqs. (35)-(37), we have used that the total energy density of the Universe can be expressed as ρtot = 3m 2 p H 2 * a −3(1+w) , with w the post-inflationary equation of state. We recall as well that H ≡ a /a. The factorization Ω GW = δ * w Θ GW in Eq. (35) is indeed very convenient: the dependence on {qs, β, w} of Θ GW (k, z), comes only from the effect of these parameters on the dynamics of the eom of the Higgs + gauge fields system, Eqs. (13)- (15).
Note that the prefactor δ * in Eq. (35), implies a suppression of the GW (energy density) as ∼ (H * /mp) 4 ≪ 1. This effect is related to the fact that the typical initial amplitude of the Higgs condensate is ϕ 2 * ∼ ϕ 2 rms ∼ H 2 * , which is then suppressed by the appearance of a Planck mass factor as 1/m 2 p in the rhs of the GWs' Eq. (27). The scaling ∝ δ * is ultimately responsible for the smallness of the GW background today, as we will emphasize later on in section VII. Note that in standard preheating scenarios, say after chaotic inflation, the inflaton and preheat fields dominate the energy budget of the universe, and have typically much larger field amplitudes. Therefore, there is no such suppression in standard preheating via parametric resonance. The production of GW from subdominant field(s), like inflationary spectator fields as in our case, will however be always suppressed by the smallness of the fields amplitude ϕ ∼ H * mp.
Depending on whether the post-inflationary equation of state is stiff, w > 1/3, or not, w ≤ 1/3, the background energy density of the Universe will correspondingly decrease slower or faster than relativistic species, i.e. d log ρtot/d log a ∝ −3(w + 1) will be < −4 for stiff backgrounds, or ≥ −4 for non-stiff backgrounds. The prefactor w = (a/a * ) 3w−1 in Eq. (37) will, therefore, either suppress the GW background as ∝ w < 1 for w < 1/3 (e.g. w = 0 for MD), or enhance it as ∝ w > 1 for w > 1/3 (e.g. w = 1 for KD). For w = 1/3 the background energy density corresponds to a RD Universe, and hence w = 1, so that there is neither a suppression nor an enhancement. As we will discuss in section VII, in a KD scenario with w = +1, the amplitude of the GW background will be maximally enhanced since w 1. However, even in this case, the large suppression due to δ * 1 will still dominate over this enhancement, so that the overall modulation of the signal is Ω GW ∝ δ * w ∼ (H * /mp) 2 , which still represents a suppression, though a milder one.
In order to solve the eom Eq. (30) for the GW, we have followed the standard procedure first introduced in [8], solving a relativistic wave-like equation in real space sourced by the full Pij, with no TT projection, We can then recover wij at any moment, in Fourier space, through the relation (40) where Λ ij,lk (k) is a geometrical projector that filters out the TT degrees of freedom in Fourier space. Since Λij,pq(k)Λ pq,lm (k) = Λ ij,lm (k), the argument inside the angular-average ... in Eq. (37) can be computed as We have studied the GW creation process in lattices of N = 256 points per dimension. To solve the Higgs + gauge fields eom Eqs. (13), (14), while verifying the constraint Eq. (15) at every time, we have used the non-compact Abelian lattice formulation presented in Ref. [33]. The exact details of the discretization formalism preserving gauge invariance are described in the Appendix A of [33], so we do not repeat them here. In the same way, a discussion of how we set the initial conditions of the different fields can be found in the Appendix B of the same paper. In all simulations we have ensured that the lattice resolution covers well the dynamical range of momenta excited in the process, for both the matter and the GW fields.
The discrete version of the GW eom Eq. (38), contrary to the matter fields eom Eqs. (13)-(15), does not follow from a discretized action. Instead, we simply substituted the continuous derivatives ∂µ in Eq. (38), with standard forward/backward lattice derivatives. In order to introduce a lattice version of the energy density spectrum of GW Eq. (37), we followed the prescription introduced in [74]. In our case, this translates into where dx ≡ H * dx is the dimensionless lattice spacing, κ(ñ) ≡ k(ñ)/H * the dimensionless momenta, k(ñ) ≡ (2π/L)|ñ| the momentum at the Fourier lattice siteñ, L the length of the lattice box, and wij ≡ wij(ñ, z) the discrete Fourier transform of wij(n, z), with n labeling the lattice sites. Note that Λ (L) ij,lm is a discretized version of the TT projector given in Eq. (40), and multiple choices are possible. We have chosen a lattice projector based on forward derivatives, noticing that other choices did not change the GW spectra appreciably, see [74] for a thoughtful discussion on this point.

V. RESULTS FROM LATTICE SIMULATIONS
In this section we present the basic features of the GW spectra produced during the post-inflationary Higgs decay process, obtained from the outcome of our lattice simulations. We leave a detailed parametrization of the spectra for Section VI, and the analysis of the redshift of the GW background until today for Section VII.
Our simulations depend on a series of parameters, some of them being unknown quantities of our system. The first unknown quantity is the initial amplitude of the Higgs field ϕ * . We know that, at the end of inflation, ϕ * changes from patch to patch with variance given by Eq. (3). Our lattice simulations do correspond to a single patch, inside which we consider the random value ϕ * to be homogeneous. The second unknown parameter is the amplitude of the Higgs self-coupling λ at the post-inflationary scales. This determines the exact form/amplitude of the Higgs potential V = λϕ 4 /4 introduced in the lattice. If we fix the strength of the gauge couplings to their value at very high energies, Taking into account the large freedom in the Hubble rate, 10 2 GeV H * 10 14 GeV, and the experimental uncertainty in the top quark mass mt (which affects the running of λ), a good physical range for these parameters is β ∈ [5 · 10 −4 , 0.3] and qs ∈ [20, 3000] (corresponding to λ ∈ [1.5 · 10 −2 , 10 −4 ]).
Note that we have also simulated values within the range 5 ≤ qs ≤ 20 for completeness, although for highenergy inflationary scales with H * of the order of (or somewhat smaller than) H max * ∼ 10 14 GeV, those values correspond to excessively high λ. Only for inflationary Hubble rates H * H (max) * we expect to obtain qs O(10); however this kills completely the GW signal Eq. (35), as the latter scales as ∝ (H * /m p ) 4 1. As we do not consider any particular inflationary model, the post-inflationary expansion rate is also unknown. We can characterize it by the equation of state w, see Eq. (10). For example, if inflation is caused by an inflaton with a quadratic potential, the Universe following inflation expands effectively as MD, with ρtot ∝ 1/a 3 . If the inflaton potential is quartic instead, it behaves as RD with ρtot ∝ 1/a 4 . We can even consider more exotic scenarios, like a KD universe, in which the energy density decays faster than that of relativistic species, with ρtot ∝ 1/a 6 . As we do not make any assumptions on the particular inflationary model, we consider w as a free parameter determining the expansion rate. In conclusion, we parametrize the GW spectra as a function of three independent variables qs, β, and w.
As described in [33], the exact dynamics of the Higgs decay process depend sensitively on the value of qs. Correspondingly, it is expected that the exact details of the GW spectra will also depend sensitively on qs. However, the qualitative aspects of these spectra can be easily understood, without the need to specify the particular value of qs. To see this, let us look at Fig. 2. There we show the temporal evolution of the spectrum Θ GW (k, z; qs, β, w). The plots correspond to the resonance parameters qs = 61 and 750, and to KD, RD and MD post-inflationary expansion rates. Within each plot, each line corresponds to the GW spectra at a particular time, showing its evolution from approximately the start of the Higgs oscillations until well after the production of GW ceases. Note that in these plots we consider the particular value β = 0.01, but a scaling of the results to arbitrary β values will be presented in the next section.
Let us now discuss three qualitative aspects of the Θ GW (k, z; qs, β, w) spectra shown in the figure: the time evolution of the spectra, the amplitude when GW stop growing, and the appearance of peaks.
Let us focus first on the time evolution of the spectra, and its relation with the time scales of the postinflationary Higgs dynamics introduced in the last section: zosc (onset of the Higgs oscillatory regime), zi (time at which the backreaction of the gauge bosons onto the Higgs condensate starts becoming effective), and ze (stabilization of the Higgs energy density and the onset of equipartition). We observe in Fig. 2 that the GW pro-duction begins shortly after the start of the Higgs oscillations, i.e. at the onset of parametric resonance at z zosc. From then on, we observe a significant growth of the GW amplitude during the linear stage z zi. This is due to the initial exponential excitation of the gauge bosons, due to the parametric resonance induced by the Higgs condensate oscillations. However, the final amplitude of the spectra is mostly determined by the non-linear dynamics during some time after the onset of backreaction z > zi, while the Higgs condensate is decaying noticeably. We can define z GW as the time scale at which GW stop being produced, so that Θ GW saturates to a fixed amplitude. In general, one finds that z GW < ze. In other words, the GW stop being produced before the onset of equipartition. This can be clearly observed in Fig. 2. In Ref. [33] we found which in the examples shown in the figure, corresponds to ze ≈ 1520, 3270, 7040 for qs = 67, and to ze ≈ 4350, 9370, 20190 for qs = 750, for KD, RD, MD postinflationary expansion rates, respectively. Note that these times are much longer than the final times displayed in Fig. 2, when the spectra have already saturated.
The fact that z GW < ze is indeed not surprising. The precise moment when GW cease to be produced is better determined when the Higgs energy density stops dropping abruptly, and this happens sometime after z = zi but before z = ze, see Fig. 1. From this time onwards (z > z GW ), even if the Higgs energy density is still decaying until the onset of equipartition at z = z e , the matter fields are only evolving smoothly, adjusting themselves towards equipartition. The time ze simply indicates when the Higgs (comoving) energy density is finally stabilized to a fixed amplitude, coinciding with the onset of equipartition. In conclusion, there is no more GW production after z = z GW . The growth of the GW spectra saturates at that moment, and the GW simply redshift from then on, due to the expansion of the Universe.
Let us now discuss the final amplitude of ΘGW after it has saturated, i.e. for z > z GW . If we focus on the panels where qs = 61, we see that, independently of the chosen post-inflationary expansion rate (either KD, RD or MD), the maximum amplitude of the GW spectra is of the same order of magnitude, ΘGW ∼ O(10 −10 ). Of course, the particular shape of the final spectra is different in each case, but the final amplitude seems to very similar. The same happens if we focus on the qs = 750 case, in which, for the three KD, RD and MD spectra, we have ΘGW ∼ O(10 −8 ). This indicates that the final amplitude of ΘGW at saturation is roughly independent on the postinflationary expansion rate.
This should not be confused, however, with the standard change of amplitude of the GW due to their nature as relativistic species. The prefactor w in Eq. (36), which verifies w 1 > w 2 if w1 > w2, accounts precisely for this effect. Therefore, the final amplitude of the GW is indeed much more affected by their natural redshifting, than by the small dependence of ΘGW on the rate of expansion. Since w has an exponential dependence on w, see Eq. (36), the proportionality Ω GW ∝ w impacts dramatically on the amplitude of the GW background today. It can affect the GW in both ways, either suppressing the amplitude by w < 1 for w < 1/3, or enhancing it by w > 1 for w > 1/3. Furthermore, this modulation of the GW amplitude will continue even after the production of GW has ceased, i.e. at z > z GW . This is because the GW only redshift at the same rate as the expand-ing background when the Universe becomes RD, hence turning the prefactor into unity, w = 1. In summary, the Ω GW ∝ w dependence means that the slower the postinflationary expansion rate (i.e. the larger the equation of state w), the higher the final amplitude of the GW background.
Let us finally discuss the appearance of peaks in the GW spectra. In Fig. 2 we can see that, during the growth of the GW spectra, a structure of peaks develops. The GW are sourced by both the Higgs and gauge fields through the terms Pij of Eq. (32), acting in the rhs of Eq. (30). In momentum space, the spectrum of GW is then sourced by a convolution of the Higgs and gauge fields spectra. Therefore, the position of the peaks is correlated with the appearance of peaks in the spectra of both the Higgs and the gauge fields.
Let us denote by u [g] ij the contribution to the GW sourced only by the gauge fields term P [g] ij (E, B), and by u [h] ij the contribution sourced by the Higgs covariant derivatives P [h] ij (Dh), see Eq. (32). From the linearity of Eq. (38), it follows that Similarly, let us denote as Θ [g] GW and Θ [h] GW the contribution to the GW spectra associated to these fields respectively. Clearly, as the GW spectrum is quadratic in uij, then GW represents an interference contribution from the convolution of a term like ∼ P [g] ij P [h] ij . In Fig. 3 we show, for the case qs = 61 and β = 0.01, both Θ [g] GW and Θ [h] GW as well as the total spectrum Θ GW for two different times. One can see that GW evolve in a similar manner, being almost identical, especially in the infrared regime. In particular, they both show some peaks at certain scales. This is a reflection of the dynamics of the system, which creates similar peaks in the spectra of Ei, Bj and Dih, correspondingly transferring those peaks into P [g] ij and P [h] ij : during the linear regime of parametric resonance, the fast creation of gauge bosons induces a similar growth of the electric and magnetic fields, as well as of the Higgs covariant derivatives. As a consequence, P [g] ij and P [h] ij contribute very similarly to the total spectrum of GW. This has in fact a very interesting effect in the GW spectrum, as it produces a clear destructive interference effect in the infrared, suppressing the total amplitude Θ GW with respect GW . At the same time, this softens (in some cases it almost washes out) the peak structure, which becomes much more smoothed in the final spectrum. This is clearly shown by the continuous curves in Fig. 3, as compared with the dashed and dotted-dashed curves.
The origin of the peaks can be understood by examining the spectra of the matter fields, i.e. of the Higgs and gauge bosons 6 . Looking at the initial stages of the process, a growth in both the Higgs and gauge fields spectra takes place in infrared scales (small k). Peaks are generated in the matter fields spectra, according to the band structure of the Lamé equation. These peaks are created during the initial stages of the process, when parametric resonance starts building up, and the Lamé equation applies. These scales are essentially imprinted in the spectrum of the GW during the excitation of the matter fields. The position of the most-infrared peak in the GW spectra, common to both the qs = 61 and 750 cases in Fig. 2, is indicated with a dotted-dashed line. It corresponds to the initial resonance band in the spectra of the gauge fields. In the qs = 750 case, there is even a second peak in the GW spectrum, indicated with a dotted line. It corresponds to another peak appearing in the spectrum of the Higgs field. When the system becomes fully non-linear, the spectra of both fields show a rescattering effect towards the ultraviolet, populating modes of higher and higher momenta. This generates a characteristic feature in the fields' spectra, which develops a relatively wide peak with a 'hunchback' shape in the ultraviolet scales. This last peak is shifted towards higher momenta according to how large the resonance parameter qs is. Again, this scale is imprinted in the GW spectrum, and it is indicated with a dashed line in both cases qs = 61 and 750 in Fig 2. A systematic study of the origin and correlation of the peaks in the GW and the matter fields spectra would be extremely interesting. Such a study could be used to probe the properties of the SM at high energies, such as the running of the couplings involved. This goes, however, beyond the scope of this paper where, as a first step, we only want to obtain a more modest characterization of the general aspects of the GW background. We will therefore adopt a phenomenological approach in section VI, characterizing the peak structure in the GW spectra, by means of simple fitting formulas.

VI. PARAMETRIZATION OF THE GRAVITATIONAL WAVE SPECTRA
In this section we parametrize the position and amplitude of the final peaks in the GW spectra as a function of q s , w, and β. We will focus first on the particular case β = 0.01, and from this, we will apply the scaling found in [33] to extrapolate and generalize this parametrization to other β values.
Let us start with the position of the peaks. We show in Figure 4 the momenta ki at which the peaks appear as a function of qs, for the β = 0.01 case, and for the different expansion rates we have simulated: KD, RD and MD. The maximum number of peaks we can observe in the spectra is three: one associated to the hunchback, whose position we denote by k3 (red circles), and two associated with the initial parametric resonance dynamics, whose position we denote by k1 (purple squares) and k2 (orange triangles). However, for some values of q s we do not see all three peaks: for qs 200 the k2 peak is not clearly observed, as it overlaps with either of the two. Also, for some qs the peaks k1 and k3 are too near to each other, and hence it is difficult to attribute a particular peak to either of them. This explains why, for some specific values of qs (particularly at low qs), we just show the red circles corresponding to k3.
The key idea is that, except for very low qs, we appreciate a clear separation between the hunchback k3 scale and the other scales k1, k2. This separation is appreciated in all the post-inflationary expansion rates. More specifically, the position of the hunchback peak increases with qs, exhibiting a clear power-law dependence. We find the following fit On the other hand, the position of k1 and k2 are mostly independent 7 on qs. We find these peaks to be well fitted by 7 In reality there is a logarithmic dependence to qs, but we would need to go to very large values qs 10 3 to start noticing it. with parameter values (again for β = 0.01) as These fits are depicted with straight lines in Fig. 4.
On the other hand, we show in Fig. 5, the amplitude of the spectrum evaluated at the highest peak Θ GW (kp), for the different qs considered, and for different postinflationary expansion rates. For β = 0.01, we find the following phenomenological fit where This peak corresponds to the maximum amplitude of the GW at the moment when they stop being actively created, i.e., at z = z GW . However, note that kp does not necessarily correspond always to the same peak k1, k2 or k3; rather, it alternates among these [for KD and RD expansion rates we normally have ΘGW(kp) ΘGW(k3), while for MD we have ΘGW(kp) ΘGW(k1)]. We see in Fig. 5 that the three fits for KD, RD, and MD coincide pretty well, confirming what we pointed out in the last section: the maximum amplitude of ΘGW at saturation time z GW is roughly independent of the post-inflationary expansion rate (the shape, however, is not; see Fig. 2).
These fits have been obtained for the particular β = 0.01 case, but a generic extrapolation to other β values can be easily carried out. We just need to use the rescaling laws that we found in [33], which connect scales and field amplitudes, from one simulation with Higgs initial amplitude and post-inflationary equation of state (β1, w1), to another simulation with different parameters (β2, w2). In particular, where Using these rescaling laws we predict the position of the peaks in the GW spectrum for arbitrary initial Higgs amplitudes β as On the other hand, using the scaling laws Eqs. (52)-(54), we demonstrated in Ref. [33] that we can recover the dynamics of the matter fields for a given set of (β, w) parameters, in terms of the results from an actual simulation done for another set (β , w ). Likewise, rescaling the terms involved in the GW source Eq. (32) by means of Eqs. (52)-(54), we can predict now the scaling of Θ GW [Eq. (43)], and, hence, how the amplitude of the background of GW scales with β. We find that We have confirmed the validity of these predictions by carrying out several lattice simulations with different β and w parameters. As an example, in Fig. 6 we show various spectra of GW for qs = 354, for both RD (w = 1/3) and KD (w = 1). The continuous red, dashed yellow, and dotted-dashed blue lines, show the spectra for β = 0.2, 0.03, 0.004 respectively, obtained directly from lattice simulations. We indicate with arrows the theoretical predictions for β = 0.2, as obtained from the β = 0.03 and β = 0.004 lattice simulations, using the extrapolation laws Eqs. (56), (57). We see that the two extrapolated predictions match quite well the output of the real β = 0.2 lattice simulations within errors.
Using Eqs. (35), (50) and (57), we obtain that the maximum amplitude of the GW background at the end of the production stage, as a function of β, qs, w, is given by where w , δ * are given by Eq. (36), and A GW , α GW by Eq. (50). The amplitude in Eq. (58) constitutes one of the key results of our analysis. However, in order to quantify the amplitude of the signal today, we need to redshift its amplitude and frequency.

VII. THE GRAVITATIONAL WAVE BACKGROUND TODAY: REDSHIFTING THE SPECTRUM THROUGH COSMIC HISTORY
We now compute how the GW background redshifts until today. We first define with a * the scale factor at the end of inflation at z = z * = 0, a RD the scale factor at the onset of the Radiation-Dominated stage of the Universe at z = z RD , and w the effective equation of state between z * and z = z RD . Essentially, RD quantifies our ignorance about the expansion rate between z * and z RD . Let us take as a frequency of reference the one corresponding to the mode kp of the highest peak of the spectrum. The frequency today associated to the peak scale kp is then given by The highest peak of the GW spectrum today is of course characterized by the highest peak of Θ GW , parametrized 8 by Eqs. (45), (50). The frequency and amplitude of highest peak today is then In order to understand what frequencies and amplitudes these expressions really imply, we need to consider specific cases. For instance, let us assume that the universe is RD after inflation, so that RD = 1, and let us consider that the inflationary Hubble rate is close to its upper bound, H * H (max) * . Taking qs = 100 and βrms 0.1, we obtain RD : h 2 Ω (o) GW (f p ) 10 −29 , at f p 3 · 10 8 Hz. (66) This amplitude is tiny, so unfortunately there is not much hope to expect to detect it in the future, unless high-frequency GW detection technology undergoes unforeseen development. The main reason why this signal is so small lies in the suppression ∝ δ * = (H * /mp) 4 4 1.
If the Universe was MD after inflation, the situation becomes even worse, because there is an extra dilution of the signal, as the latter is now proportional to some factor RD 1, which becomes smaller and smaller the longer it takes for the Universe to reach a RD regime at z = z RD . This dilution is simply a consequence of the fact that GW scale with the expansion of the Universe as relativistic species, ρ GW ∝ 1/a 4 , whereas a MD background energy density dilutes slower as ρ ∝ 1/a 3 .
If the Universe is KD after inflation, the GW signal is, however, enhanced significantly. In particular, given the initial ratio of energies ∆ ≡ V * /ρ * ∼ 10 −12 [Eq. (5)], the Universe will sustain a KD expansion rate until the moment when the relativistic SM fields become dominating the energy budget. This implies that the GW signal is enhanced by a factor ∝ RD = 1/∆ ∼ 10 12 . The scaling of the signal also goes as ∝ (β/0.01) 4+v(1) with v(1) = 2/3, instead of with v(1/3) = 0 as in RD. In addition, A KD GW A RD GW . Compared to a RD background, and for β = 0.1, there is therefore another enhancement (however milder) 8 Although the highest-amplitude peak kp is normally k3, this is not always the case. However, when kp is instead associated with k1 or k2 (typically for low qs), the spectral amplitude at the k3 peak is still very similar to that of the highest peak. Therefore, for simplicity, we are going to associate here the amplitude This corresponds yet to a small signal, but its amplitude is in fact comparable 9 to the standard scale-invariant inflationary background h 2 Ω (Inf) GW 5 · 10 −16 (H * /H (max) * ) 2 . Our signal in this case of KD after inflation, however, lies at extremely high frequencies ∼ 10 11 Hz, beyond the range of planned GW detectors.
Before we move into the concluding section, it is perhaps interesting to make a remark about a particular aspect of this GW background, which we have not studied in this paper. Given the condition of the Higgs as a condensate with a finite correlation length at the end of inflation, it is indeed expected that this GW background will have anisotropies on angular scales corresponding to that correlation scale. This is very similar to the case of preheating scenarios with light preheat fields [75,76]. Of course, given the smallness of the background amplitude itself, detecting such anisotropic variation in the sky seems even more chimeric than detecting the background itself. Yet, it is interesting to note that these anisotropies are expected in the present case, as opposed to the situation when the decaying (oscillatory) field is an inflaton, which is considered to be homogeneous 10 over scales much larger than our present horizon.

VIII. SUMMARY AND DISCUSSION
Independently of the nature of the inflationary sector, a stochastic background of GW is expected simply due to the existence of the Standard Model Higgs. This background of GW is always generated after inflation, as long as the Higgs is decoupled from (or sufficiently weakly coupled to) the inflationary sector. In such a case, the Higgs is excited either: i) during inflation, if it is minimally coupled to gravity (or if it has a weak non-minimal coupling), or ii) at the end of inflation, if it has a (sufficiently) large non-minimal coupling to gravity. Either way, we expect the Higgs to be in the form of a condensate after inflation, decaying very rapidly -via non-perturbative effects -into the rest of the SM species. The resulting post-inflationary out-of-equilibrium dynamics of the SM fields, generates necessarily a stochastic background of GW. Both the SM Higgs and the electroweak gauge bosons act as the dominant sources of the GW background.
We have studied the details of the form of the GW spectrum, determining its frequency, amplitude and shape. We have characterized the GW spectrum dependence on the unknown parameters of the system, namely the Higgs initial amplitude at the end of inflation β = √ λϕ * /H * , the equation of state w characterizing the post-inflationary expansion rate of the Universe, and the resonance parameter qs = (g 2 Z + 2g 2 W )/4λ. The running of the Higgs self-coupling at high energies is in fact quite uncertain within the experimental input, so λ can vary within the range 10 −2 λ < 10 −5 . This translates into some uncertainty in the regime of the resonance parameter, which may vary within the range qs ∼ O(10)−O(10 3 ).
We have used real-time classical gauge field lattice simulations in an expanding box in (3 + 1) dimensions. We chose N = 256 points per dimension, ensuring that the relevant modes involved in this process were well captured within the dynamical range of the simulations. Our results have been obtained within an Abelian-Higgs modeling, expected to describe sufficiently well the system when qs 1. Only in the case of the smallest resonance parameters qs ∼ O(10), does one expect the dynamics of the system to be affected by the presence of the full non-Abelian structure of the SM.
From our lattice simulations, we have obtained Eq. (58), which is a phenomenological fit of the amplitude of the GW spectra as a function of the different unknown parameters described above. We also obtain a parametrization of the observed redshifted amplitude until today in Eq. (65). However, the GW signal is suppressed by the inflationary Hubble rate as ∝ (H * /mp) 4 . The largest amplitudes for the GW background are therefore obtained when H * is only somewhat smaller than (but of the order of) its current upper bound H (max) * ∼ 10 14 GeV. This implies that λ runs to small values λ < 10 −2 , hence making the resonance parameter large, qs > 10. In light of this, the use of the Abelian approach is fully justified. In any case, the basic features of the fields dynamics and GW production, i.e. its dependence on qs, β and w, are not expected to change drastically in the full non-Abelian scenario. Our study can be considered therefore as a good indicator of the GW amplitudes to expect in general, even if non-Abelian corrections were to be considered.
If the Universe was RD after inflation, our calculations show in fact that this background is tiny, with an amplitude of h 2 Ω (o) GW (fp) 10 −29 , and peaked at high frequencies fp ∼ 300 MHz. The smallness of this background reflects simply the fact that the initial energy of the Higgs condensate represents only a tiny fraction of the inflationary energy. If the Universe was MD after inflation, although the background will be peaked at slightly smaller frequencies, its amplitude today can only be even smaller than in the RD case. The amplitude of the background is expected, however, to be enhanced significantly if the Universe underwent a KD regime after inflation. In that case, our calculations show that the background today could have an amplitude up to h 2 Ω (o) GW (fp) 10 −16 . This larger background is, however, peaked at very high frequencies, of the order of fp 10 11 Hz.
The generation of the GW background we have studied in this paper is, in some sense, universally expected, as long as the SM is not strongly coupled to the inflationary sector. However, given that the background is always peaked at very high frequencies, and its amplitude today is very small (in all cases), our prediction will remain, most likely, as a curiosity of the SM.