Higgs-curvature coupling and post-inflationary vacuum instability

We study the post-inflationary dynamics of the Standard Model (SM) Higgs field in the presence of a non-minimal coupling $\xi|\Phi|^2R$ to gravity, both with and without the electroweak gauge fields coupled to the Higgs. We assume a minimal scenario in which inflation and reheating are caused by chaotic inflation with a quadratic potential, and no additional new physics is relevant below the Planck scale. By using classical real-time lattice simulations with a renormalisation group improved effective Higgs potential and by demanding the stability of the Higgs vacuum after inflation, we obtain upper bounds for $\xi$, taking into account the experimental uncertainty of the top-Yukawa coupling. We compare the bounds in the absence and presence of the electroweak gauge bosons, and conclude that the addition of gauge interactions has a rather minimal impact. In the unstable cases, we parametrize the time when such instability develops. For a top-quark mass $m_t \approx173.3 {\rm GeV}$, the Higgs vacuum instability is triggered for $\xi \gtrsim 4 -5$, although a slightly lower mass of $m_t \approx 172.1 {\rm GeV}$ pushes up this limit to $\xi \gtrsim 11 - 12$. This, together with the estimation $\xi \gtrsim 0.06$ for stability during inflation, provides tight constraints to the Higgs-curvature coupling within the SM.


I. INTRODUCTION
The Standard Model (SM) potential may become negative at very high energies [1,2]. This has prompted an important effort to determine whether the electroweak vacuum is, in the present, stable or unstable. Current measurements of the top quark and Higgs masses suggest that we live in a meta-stable Universe: the probability of the Higgs field to decay into a higher-scale negativeenergy vacuum is non-zero, but the estimated decay time is much larger than the present age of the Universe [3].
However, the situation is quite different in the early Universe. In this case, high energies and high spacetime curvature can make the vacuum more unstable. In particular, this may happen during inflation [3][4][5][6][7][8][9][10][11], or during the successive period of (p)reheating [12][13][14][15][16]. The dynamics of the Higgs field Φ during and after inflation, as well as the potential instability of the Higgs vacuum, depend very sensitively on the strength of its non-minimal coupling to the scalar curvature, defined as ξ|Φ| 2 R, with R the Ricci scalar. This interaction is necessary to renormalise the theory in curved space [17,18], and given that ξ runs with energy, it cannot be set to zero at all energy scales. Gravitation is very weak in comparison with the other interactions, so current particle-physics experiments provide only very weak constraints to this coupling [19]. The coupling ξ can be considered, therefore, as the last unknown parameter of the SM.
If ξ 0.1 the Higgs is effectively light during inflation, and behaves as a spectator field, forming a condensate with a large vacuum expectation value (VEV) [3,[20][21][22]. If it exceeds the position of the potential barrier, the Higgs reaches its true negative-energy vacuum and generates patches of anti-de Sitter space, resulting in a catastrophic outcome for our Universe [3][4][5][6][7][8][9][10][11]. One simple way to prevent this from happening is to consider a suf-ficiently low inflationary scale, so that even if the Higgs is excited during inflation, its amplitude never reaches the potential barrier. Another way of ensuring vacuum stability is to consider values of the top quark mass 2-3 sigma below its central value m t = 172.44 ±0. 13(stat) ±0.47(syst) GeV [23], so that the second minimum in the Higgs potential is either shifted to sufficiently high energies, or it is simply not present. In any case, if the Higgs field remains stable during inflation, it starts oscillating around the minimum of its potential short after inflation ends, rapidly decaying into the SM gauge bosons and fermions via nonperturbative parametric effects [22,[24][25][26]. This may have relevant cosmological consequences, like the realization of successful reheating into the SM without additional mediator fields [27], the realization of baryogenesis via leptogenesis in certain extensions of the SM [28][29][30], or the production of a primordial background of gravitational waves (GW) peaked at high frequencies [24,31].
On the other hand, if ξ 0.1, the height of the potential barrier increases, and the Higgs is no longer a light degree of freedom during inflation [3,7]. In this case, the Higgs field adquires an effective mass of the order m 2 Φ ξR ∼ 12ξH 2 * H 2 * during inflation, with H * the inflationary Hubble rate. This prevents the Higgs from developing large amplitude fluctuations during inflation. However, the situation is quite the opposite after inflation ends. The post-inflationary oscillations of the inflaton φ around the minimum of its potential induce rapid changes in the spacetime curvature R, which becomes negative during a significant fraction of time in each oscillation. The effective mass of the Higgs field becomes tachyonic during those moments, m 2 Φ ∝ R < 0. If ξ is sufficiently large, the Higgs field may be significantly excited during the tachyonic periods, potentially triggering the vacuum instability [12]. This issue has been studied lately, using both analytical and numerical techniques, as well as classical real-time lattice simulations [12][13][14]. The results of all of these works agree qualitatively, finding ξ O(1) − O (10) as an upper bound for achieving stability after inflation. A similar lattice analysis of the values of the Higgs-inflaton coupling inducing the instability of the Higgs vacuum has also been carried out in [15], while an analysis of the combined effects of both Higgs-curvature and Higgs-inflaton couplings has been done in [16].
In this work, we use classical field theory lattice simulations to constrain the range of allowed ξ values which ensure the stability of the Higgs vacuum after inflation. We do a systematic parameter analysis of the Higgs postinflationary dynamics. We use in the simulations the renormalization group improved Higgs effective potential, and study the impact of the initial conditions and number of Higgs components in the results. We include also an analysis of how the time scale at which the Higgs field develops the instability depends on ξ and the topquark mass.
Furthermore, we consider the more realistic situation (so far not considered in the literature) where the Higgs field is coupled to the electroweak gauge bosons. We mimic the SM gauge interactions an Abelian-Higgs analogue model, but argue that this is enough to demonstrate the effect of restoring the gauge interactions. Our Abelian set-up captures well the gauge boson field effects onto the Higgs post-inflationary dynamics, as we expect the non-Abelian terms of the Lagrangian to be subdominant, especially at the earliest times.
We have assumed throughout this work a chaotic inflationary model 1 2 m 2 φ φ 2 , with m φ 10 13 GeV fixed by the observed amplitude of the CMB anisotropies [32]. We note that, although the simple chaotic inflation model with a quadratic potential is in tension with observational data [32,33], it is possible to modify the large filed behavior to make it consistent with observations [33,34], without changing its post-inflationary dynamics.
The structure of the paper is as follows. In Section II we present a brief review of the inflaton and Higgs dynamics after inflation in the presence of a Higgs-curvature non-minimal coupling. In Section III we present the equations of motion and the initial conditions of the different fields, as well as some qualitative aspects of our lattice simulations. The following three sections present the results from our lattice simulations, with increasing degree of complexity. In Section IV we consider a free scalar field with no potential. This is useful to understand better the results in Section V, where we introduce the renormalisation group improved Higgs potential. We determine the values of the coupling ξ that give rise to an unstable Universe, and parametrize the time scale at which the instability takes place, as a function of ξ and m t . In Section VI we repeat the same analysis, but including also the gauge bosons in the lattice. In Section VII we discuss our results and conclude. Finally, in the Appendix we discuss the details of the lattice formulation of the system. In this paper we use signature (−, +, +, +), and hence consider a flat background with Friedman-Lemâitre-Robertson-Walker (FLRW) metric like ds 2 = −dt 2 + a 2 (t)dx i dx i , where a(t) is the scale factor, and t is the cosmic time. We refer to the reduced Planck mass as m p = (8πG) −1/2 2.44 · 10 18 GeV.

II. INFLATON OSCILLATIONS AND HIGGS PARTICLE CREATION
We consider throughout the paper the inflationary chaotic model V (φ) = 1 2 m 2 φ φ 2 , where φ is the inflaton, and its mass is fixed to approximately m φ = 1.5 × 10 13 GeV in order to explain the CMB anisotropies. If φ O(10)m p , the field is in a slow-roll regime, causing the inflationary expansion of the Universe. However, when H(t) ≈ m φ with H(t) the Hubble parameter, the inflaton field starts oscillating around the minimum of its potential, ending the inflationary stage. Let us define t * as the time when H(t * ) = m φ holds exactly. This moment signals the onset of the inflaton oscillations. The coupled equations of motion of the inflaton and scale factor areφ To obtain the initial conditions for the homogeneous inflaton, we have solved numerically the coupled inflaton and Friedmann equations, Eq. (1) and Eq. (2), imposing the slow-roll conditionsφ −m 2 φ φ 2 /3H 2 ,φ m 2 φ φ 2 well before the end of inflation. From the numerical solution, we obtain the time t * at which H = m φ holds exactly. At this moment we find Using Eqs. (1)-(2), the Ricci scalar can be expressed in terms of φ andφ, The inflaton field after inflation behaves, approximately, as a damped oscillator with decaying amplitude [35] Each time the inflaton field crosses around zero, φ ≈ 0, we have R(t) < 0 from Eq. (4). This can be clearly seen in Fig. 1, where we plot both the inflaton and the Ricci scalar as a function of time.
Let us focus now on the post-inflationary dynamics of the Higgs field. The relevant piece of the Standard Model for us is , ξ is the strength of the coupling of the Higgs to the Ricci scalar, F µν the gauge field strength (assumed Abelian for simplicity), and D µ = ∂ µ − iA µ is a gauge covariant derivative, which describes the interaction of the Higgs field with the electroweak gauge bosons. Due to the coupling to the scalar curvature, the Higgs gets an effective mass m 2 Φ (t) = ξR(t). Therefore, the Higgs becomes effectively tachyonic with m 2 Φ < 0, during the intervals when the Ricci scalar becomes negative. Due to this, there is a strong excitation of the Higgs field, a phenomenon known as tachyonic resonance [36].
We can estimate both the period of time that the Ricci scalar becomes negative, as well as the maximum momenta excited by the resonance. The inflaton crosses zero periodically at m φ t n = (n−1/2)π, n = 1, 2, 3, .... We can determine a typical envelope amplitude between the nth and the (n + 1)-th crossings, as φ a /m p = 8/3/πn. When the inflaton crosses around zero, the Ricci scalar On the other hand, the greater the coupling ξ, the larger the range of Higgs tachyonic modes excited while the curvature is negative. We estimate this as an infrared (IR) band from k = 0 up to a cutoff Λ, where a n is the scale factor at t n (initially a 1 1). Let us consider the unitary gauge, so that the SM Higgs doublet can be written as a real degree of freedom, Φ = ϕ/ √ 2. Let us redefine the Higgs amplitude as h ≡ ϕ/a 3/2 so that in cosmic time, this re-scaling eliminates the friction term in the Higgs EOM. If we ignore the presence of the gauge bosons and of the Higgs self-interacting potential, the equation of motion of its Fourier modes is where ∆ ≡ − 3 4ȧ 2 a 2 − 3 2ä a , so that ∆ k 2 /a 2 for subhorizon scales. We can then set ∆ → 0, and using Eqs. (4)- (5), write the previous EOM as a Mathieu equation where z ≡ m φ (t − t * ) and These equations of motion have been extensively studied in the context of parametric resonance in m 2 φ φ 2 preheating (see for example [35,37]). The main difference with respect to standard parametric resonance is that now, we are not constrained to the case A k > 2q, and hence we have great resonance bands which induce a stronger particle creation effect in the broad resonance regime q 1. However, note that due to the expansion of the Universe, φ a (z) decreases, and hence this pushes the Higgs into a narrow resonance regime, where this effect is much weaker. The dynamics of this theory was studied in [7] with the properties of tachyonic resonance of [36], and after that numerically in [13] and in the lattice in [14][15][16].

A. Higgs potential
Let us now consider the effect of plugging back the Higgs potential. In particular we consider the renormalisation group improved Higgs potential valid for large field amplitudes above the electroweak scale, ϕ v ≡ 246 GeV. Here λ(ϕ) is the renormalised Higgs self-coupling at the renormalisation scale µ = ϕ, where the running behavior has been computed up to three loops [1,2]. The running is very sensitive to the strong coupling constant α s , Higgs mass m h , and the Yukawa top coupling y t , the latter being currently the strongest source of uncertainty. We show in Fig  the public package of [38]. We observe that the Higgs potential possesses a maximum (a barrier) at a given scale ϕ + , and crosses zero at ϕ o , becoming negative for higher energies. We show these values for different top quark masses in Table I If λ > 0, the Higgs tachyonic resonance effect weakens, whereas if λ < 0, the tachyonic effect is enhanced. We need to take into account both cases in the lattice simulations due to the running of λ with energy.

III. LATTICE FORMALISM
We describe now our lattice formulation, which will be applied in the following three sections.
The equations of motion in the continuum can be derived from the minimization of action (6). We take the Higgs as a complex doublet, with four real components As we will argue later, we have neglected the purely nonabelian terms of the gauge field self-interactions. The covariant derivative is then simply written as The equations of motion, in the temporal gauge A 0 = 0, areΦ The third of these equations is the Gauss constraint, a relation between fields that must be obeyed at all times. In the lattice we solve a discrete version of equations (14)- (16), obtained from a discrete gauge-invariant action, see Eqs. (38), (39). Details of the lattice formalism are given in the Appendix. Note also that this is not, strictly speaking, an Abelian-Higgs model, as we are introducing two Higgs complex fields instead of simply one. The form of a(t) in these equations, as well as the Ricci scalar R(t) = R[a,ȧ,ä], is obtained from the selfconsistent solution of the inflaton and Friedmann equations (1) and (2). As we shall see, for the values of ξ we consider in this work, the energy of the Higgs field is always several orders of magnitude subdominant with respect the one of the oscillating inflaton. Hence, we just ignore the contribution of the Higgs field to the Friedmann equation. Note that the inflaton is taken as a homogeneous field, and we do not introduce it explicitly in the lattice, it simply dictates the form of a(t) and R(t) as a function of time.
In Section IV we study tachyonic resonance in the lattice, taking the Higgs as a free field without selfinteraction. The Higgs will then be excited only due to the rapidly changing spacetime background. In Section V we re-introduce back the Higgs potential, but ignore yet its interaction with the gauge bosons. We determine under those circumstances, what values of ξ lead the Higgs field to become unstable, so that it rolls rapidly into the true vacuum. In Section VI we finally incorporate a gauge structure into the simulations, and study their effect on the post-inflationary Higgs dynamics, reevaluating again the critical values of ξ.
Choosing µ = ϕ, we obtain the renormalisation group improved effective potential 1 . Mimicking (11) with µ = ϕ, we introduce in the lattice the Higgs potential where we assume |ϕ| v. As seen before, the Higgs selfcoupling λ(ϕ) runs with the value of ϕ. We introduce this running in our simulations as a local function of the lattice point n, i.e. λ(ϕ(n)): as the value of |ϕ| changes from lattice point to lattice point, so does too the value of λ. More specifically, we introduce a quartic logarithmic polynomial λ in (|ϕ|) = 4 n=0 c n (log |ϕ|) n , interpolating the 3-loop calculation of the running obtained in [38] for the relevant range of Higgs amplitudes |ϕ| (see Fig. 2). As we have mentioned, the running of the potential depends strongly on the value of the top quark mass, the current world average being m t = 172.44 ±0.13(stat) ±0.47(syst) GeV [23]. We take this uncertainty into account by providing different sets of {c n } constants, for different interpolations of the running corresponding to each value of m t .
The interpolation can only describe appropriately the running of λ for certain values of |ϕ|, failing at very low and very high field amplitudes. This mismatch is not a problem, because we have checked that those values are never achieved anywhere in the lattice, before the instability of the Higgs field is developed. On the other hand, when the Higgs has become unstable and decays towards the negative-energy vacuum, the amplitude of the Higgs field starts increasing very fast, achieving the region in which the interpolation fails. However, our aim in this work is to determine the specific time at which the instability develops, not to characterize the dynamics of the Higgs field once the instability has started, so this is neither a problem. In fact, in order to ensure numerical stability during the Higgs field transition, it is convenient to modify the high-energy running of λ, so that it generates a second vacuum at lower energies than those at which it should really be according to the running predicted in the Standard Model. This is achieved for c 4 > 0. In particular, we have chosen the constants so that the negative-energy vacuum is generated at approximately ϕ = ϕ v ≈ 10 16 GeV. If the Higgs field goes to the vacuum at ϕ = 0, the Higgs field is stable, while if it goes to the negative-energy vacuum, we say it has become unstable. We have explicitly checked that our characterization of the times of instability is independent on the particular choice of constants c n , as long as they fit the Higgs effective potential within the range ∼ 10 9 − 10 14 GeV.

A. Initial conditions
We start the lattice simulations at time t = t * , where we impose for all four components of the Higgs that their initial homogeneous amplitude vanishes, h n (t * ) = 0 (n = 1, 2, 3, 4). We then add on top a spectrum of fluctuations, which mimic the spectra of quantum vacuum fluctuations 2 , where a * ≡ a(t * ), and R * ≡ R(t * ) ≈ 10H 2 * from Eq. (4). The spectra of quantum fluctuations (19) is set in the lattice in a similar way as Latticeeasy [39]: We impose in momentum space the following spectra for the Higgs field and field derivatives where ω k,n ≡ (k/m φ ) 2 + ξR * , θ n1 and θ n2 are real phases which vary randomly from lattice point to lattice point in the interval θ n1 , θ n2 ∈ [0, 2π], and |h n | varies according to the probability distribution The ultraviolet cutoff k c is introduced in order to prevent the excitation of UV modes which are not expected to be excited by the tachyonic resonance, i.e. k c ≈ Λ where Λ is given by Eq. (7). Hence, the variance of (a component of) the Higgs field initially is with ν = 9/4 − (ξR * /H * ) 2 . However, for the couplings ξ > 4 we are considering, this spectra is almost identical to the FLRW case described by Eq. (19).
where we have taken a * = 1 in the second equality. Typical numbers chosen in our simulations are ξ ∼ 10 and k c ∼ 10H * , which give for the initial Higgs amplitude Typically ϕ 2 * ϕ + , and hence most or the entire Higgs field is already in the right hand side of the barrier when initial conditions are set. This, however, does not mean that the Higgs field will immediately become unstable, as the mainly positive sign of R may impede it. We shall discuss this issue in more detail in Sections V and VI. Let us also remark that this way of fixing the initial conditions is only appropriate if the tachyonic resonance regime of the system enhances the Higgs amplitude significantly over the value given in Eq. (23). If it does not, we cannot trust the lattice approach. Finally, let us also note that there is a contribution to the Higgs effective mass from its self-interactions, i.e. the effective Higgs mass should be better m 2 eff ≈ ξR * + λ ϕ * 2 . However, taking λ ≈ −0.01, ξ ≈ 10, and H * = m φ ≈ 6 × 10 −6 m p , we can check the second term is negligible with respect the first one.

IV. SIMULATIONS WITH A FREE SCALAR FIELD
In this Section, we study the case of a free, noninteracting scalar field, i.e. we solve only the first equation in Eqs. (14)-(16), setting λ = g s = 0. Although this is obviously not a physical case, it will be helpful to understand our later results better when we include the Higgs self-interactions potential. Therefore, we consider now a 4-component Higgs field, coupled to the spacetime curvature through the term ξRϕ 2 , with R[φ,φ] evolving due to the oscillating inflaton. We have done several lattice simulations of this system, varying the coupling ξ within the range ξ ∈ [4,70].
We show in Fig. 3 the spectra of the Higgs field for the particular cases ξ = 5 and ξ = 30. In both panels, red color corresponds to early times, while dark blue correspond to late times. In these spectra, a cutoff has been put in the distribution of initial fluctuations at a certain scale k c , as indicated in Eq. (20). The value of k c has been estimated from a previous set of lattice simulations without cutoff, in which we see that for k > k c , the field excitation due to the tachyonic resonance is negligible. Both spectra grow very fast, saturating eventually at a time t ≈ t res , defined below. Naturally, the spectra grows several orders of magnitude more in the ξ ≈ 30 case than in the ξ ≈ 5, as the tachyonic effect is stronger in the first case.
We show in Fig. 4, for the couplings ξ = 3, 6, 10, 15, 30, the Higgs conformal and physical amplitudes as a function of time, averaged over the whole volume of the lattice. We remind that this plot is for a four-component Higgs field, while for a single component we have ϕ 2 n ≈ ϕ 2 /4 for each n = 1, 2, 3, 4.
We expect the Higgs excitation to end when q 1, see Eq. (10). Using Eq. (5) and taking q ≈ 0.2 as the condition for the end of the tachyonic resonance regime, we find this time t res to be given by In the figure we indicate this with vertical dashed lines. We see that for t t res particle creation is exponential, and the greater the ξ, the stronger the growth of the conformal Higgs amplitude aϕ. However, as we approach t ≈ t res , the Higgs excitation stops. From then on, the dynamics of the Higgs field is dominated by the expansion of the universe. More specifically, we have found that the late-time behaviour of the Higgs amplitude is where the particular numerical value of the exponent depends on the particular value of ξ considered. We indicate this in the left panel of Fig. 4 with dashed lines. As expected, Eq. (25) indicates that ϕ ∝ a −1 (t). We have found that a rough estimate for the Higgs amplitude for late times is where the first factor accounts for the initial excitation of the Higgs modes, and the second accounts for the later energy dilution, given in Eq. (25). Before we move on, it is important to note that, as we decrease ξ, the amplitude of the excited IR modes decreases significantly, being comparable to the amplitude of the (non-excited) UV modes for very low couplings. This signals that the lattice simulations cannot be trusted for these low couplings, because there is no significant excitation of the Higgs field over the initial vacuum fluctuations. Correspondingly, for these low couplings, the contribution of the UV modes to the Higgs amplitude becomes increasingly important, and hence its value can depend strongly on where we put the cutoff k c of the initial fluctuations. Therefore, there is a minimum value ξ for which we can trust the lattice simulations. In this paper we have determined this condition as ϕ(tres) ϕ * a(tres) a * > 2, which means basically that the contribution to the Higgs amplitude from the Higgs excitation, is greater than the one from the Higgs initial vacuum fluctuations. With this, we find that we cannot trust simulations with ξ 4.

V. SIMULATIONS WITH UNSTABLE POTENTIAL
Let us now move to simulations with the full Higgs potential (11), including the four components of the Higgs field but yet without including gauge bosons. All the results of this Section have been obtained with lattice cubes of N 3 = 256 3 points, and minimum momentum p min = 0.18m φ .
To get a qualitative understanding of the dynamics, let us recall the linearised equation of motion (8) for the Higgs field modes h ≡ ϕa 3/2 . For high Higgs field values, ϕ > ϕ 0 , the self-coupling is negative λ(ϕ) < 0, and therefore the interaction term tends to increase the Higgs field value, and induce a transition to the negativeenergy vacuum. The more the Higgs field has been amplified by the tachyonic resonance, the faster the instability is. On the other hand, because the Ricci scalar remains larger time positive than negative each inflaton oscillating period, the non-minimal coupling term ξR(t) effectively creates a potential barrier than resists this increase. The amplitude of the curvature term decays as ξR(t) ∝ a −3 (t) ∝ t −2 , so it becomes gradually less important. However, if it delays the instability until the Higgs field amplitude has decreased below the barrier scale ϕ < ϕ 0 , then the Higgs field remains stable throughout the entire evolution. Because the amplification by the tachyonic resonance depends exponentially on the nonminimal coupling ξ [see Eq. (26)], whereas the effective barrier due to ξ depends on it only linearly, one expects that for high ξ, the instability takes place faster, and for low enough ξ it is prevented completely. Fig. 5 shows the volume-averaged amplitude of the Higgs field |ϕ| as a function of time, for different choices of the Higgs-curvature coupling ξ, obtained directly from lattice simulations. In this Figure, we have used the running of the potential corresponding to the top quark mass m t = 172.12 GeV, see Fig. 2. This potential has the barrier at ϕ + ≈ 7.8 × 10 11 GeV. We can see that, for initial times z 10, the amplitude grows (in a oscillating way) due to the Higgs tachyonic resonance regime, as described in Section II.
We see that for high values of the non-minimal coupling, ξ ≥ 16, the Higgs field becomes unstable during the tachyonic resonance, triggering a transition to the high-energy vacuum ϕ = ϕ v . For lower values of the non-minimal coupling, the tachyonic resonance ends before the Higgs has become unstable. After this the be- FIG. 6: Higgs field spectra κ 3 a 2 |ϕ k | 2 as a function of κ = k/m φ in the presence of a unstable potential (Section V), for ξ = 12.2 and mt = 172.12 GeV. The spectra is depicted at times m φ (t − t * ) = 0, 10, 20..., going from early times (red) to late times (dark blue).
haviour is initially similar to the free field case discussed in Section IV: the system settles in a quasi-equilibrium state in which the field amplitude gradually decreases due to the expansion of space. In the intermediate range of couplings, 12.2 ≤ ξ ≤ 14, the instability eventually takes place, at a time that we denote by t i . We indicate this with a vertical dashed line in Fig. 5.
For ξ ≤ 12, the field amplitude eventually decreases below the potential barrier, ϕ < ϕ + . By this time, the barrier stabilises the state, and therefore the instability does not take place at all. This demonstrates that physically the instability is due to the tachyonic resonance. Even though the amplitude of the initial vacuum fluc- we plot an oscillation-average to compare both terms more easily. These are the terms that appear in the Eq. (12) for the field modes. We also plot, with vertical lines, the corresponding time m φ ti at which the Higgs becomes unstable.
tuations is higher than the barrier scale, it is not high enough to lead to an instability before it is damped to safe values by the expansion of space. From the spectra shown in Fig. 6 we can see that the infrared modes have to be amplified by roughly three orders of magnitude by the tachyonic resonance in order for the instability to happen. In particular, this means that the use of classical field theory simulations is well justified in this case.
As explained above, we can explain the triggering of the Higgs instability in terms of the balance between the terms ξR and −λ(ϕ) ϕ 2 that appear in the EOM of the field modes, Eq. (12). We have plotted in Fig. 7 the timeevolution of these two terms for different ξ and top mass m t = 172.12GeV. Although ξR is periodically oscillating between positive and negative values, the resulting oscillation average is always positive. We observe that initially, the first term dominates over the second, but as commented, when (if) the absolute value of the second term becomes of the same order of magnitude that the first one, the Higgs field becomes unstable. This can happen during the initial regime of tachyonic resonance, or later on due to the different dilution rates of both terms.
In conclusion, as expected, we can define a critical coupling, ξ c ≈ 12 for m t = 172.12 GeV, so that for ξ ξ c the Higgs field is always stable, while for ξ ξ c the field becomes unstable at a certain time m φ t i , whose numerical value decreases as ξ gets greater. This general picture also applies for other values of the top mass. If we take the top quark mass a bit higher, ϕ + is lower, and hence the Higgs field takes a much longer time to settle on the safe side of the potential barrier. Due to this, for these masses the numerical value of the critical coupling ξ c is lower.
The order-of-magnitude fit of the time-dependence of the amplitude obtained for the free case in Eq. (26) also holds quite well in the self-interacting scenario, before the instability takes place. This signals that the effect on the Higgs dynamics of λ is not very important before the transition to the high-energy vacuum takes place. Inverting this fit, we can find an order-of-magnitude estimate of the time t o at which we recover λ(ϕ) > 0, with where ϕ o is given in Table I Fig. 8 the instability time as a function of ξ obtained from our lattice simulations. We have observed that the specific value of m φ t i depends on the particular random realization of the Higgs field initial conditions in Eq. (20), so for each point, we have done several simulations for different realizations of the initial conditions (this is discussed in more detail in Section V A). Points indicate the average value of m φ t i , while the shadow region surrounding each of the curves indicate the standard deviation, assuming a Gaussian distribution.
The behaviour of the five curves with ξ is quite similar. In all curves we can identify two critical values, ξ  Table  II. The meaning of these values is as follows: • ξ > ξ (2) c : For these values, we observe that the Higgs field always develops an instability, at a time m φ t i O(10), quite independently on the value ξ (at least for the cases we have simulated). This is seen as a plateau in the right part of the numerical curves shown in Fig. 8. Qualitatively, for this range of values, the Higgs field becomes unstable when it is still in the tachyonic resonance regime. One can see an example of this in Fig. 5 for m t = 172.12 GeV: for the cases ξ = 16, 18, 20, which verify ξ > ξ (2) c ≈ 14, the Higgs becomes unstable in the oscillatory regime, while for ξ = 12, 14, with ξ ξ (2) c , the instability is developed when the resonance has already finished.
• ξ (1) c < ξ < ξ (2) c : For these values, the Higgs field also develops an instability, which makes it end in the negative-energy vacuum, but this happens only after the tachyonic resonance has ended. However, for these values, the instability time m φ t i depends very strongly on the value ξ. A change of few units in ξ changes m φ t i in several orders of magnitude.
• ξ < ξ (1) c : Finally, for these values, we observe that the Higgs field is always stable, coming back to the safe side of the potential without having become unstable.
We indicate the values of both ξ   and ξ (2) c , obtained from lattice simulations, for different values of the top quark mass. The error in ξc signals the uncertainty with respect initial conditions. The meaning of this parameters is explained in the bulk text. red curves) in Table II. Note that, as expected, as we increase the value of the top quark mass, the position of the barrier in the Higgs potential moves to smaller field values, and hence the initial distribution of the Higgs field is much deeper in the negative-energy region. Due to this, ξ (1),(2) c are lower, and the Higgs field takes much longer to enter into the safe side of the potential. Let us note that the identification of these critical values is not unambiguous, and in particular, for couplings near the critical one ξ ≈ ξ (1) c , we observe that depending on the specific realization of the initial conditions, the Higgs may or may not become unstable. This source of uncertainty is indicated with a ± sign in Table II. Finally, let us note that our technical definition of the second critical coupling is such that for ξ > ξ (2) c , we have m φ t i < 20. The curved, dashed lines in Figure 8 indicate the approximated time at which the Higgs will enter into the safe side of the potential, using Eq. (27). The idea is that at the critical coupling ξ = ξ (1) c , the curve for m φ t i obtained from the numerical simulations (bands in colors in Fig. 8) will meet approximately the corresponding dashed ones. We can see in Fig. 8 that this works relatively well, taking into account that Eq. (27) is only an approximation.
In Fig. 8 it can also be seen that for m t = 173.95 GeV and m t = 174.56 GeV, the instability curves do not meet their corresponding curved-dashed lines for ξ 4, which are the cases that we cannot study in the lattice as discussed at the end of Section IV. Hence, for these masses we can only provide the upper bound ξ c 4.
As a final comment, let us note that in all our simulations we made the inflaton to oscillate indefinitely, even though this is clearly not realistic. The inflaton is expected to be coupled to other species which will eventually induce its decay due to parametric resonance effects at a time that we denote by t br . After this time, the energy density is no longer dominated by a coherently oscillating scalar field, and therefore Eq. (4) is no longer valid. This puts an end to the tachyonic resonance regime of the Higgs field. Therefore, the estimates for ξ c provided here will not be valid if t br t res . For example, as seen in [40], if the inflaton is coupled to a single scalar field χ with coupling g 2 φ 2 χ 2 , one finds m φ t br 40 for c , see Table II. For each data point we have done several lattice simulations corresponding to different realizations of the initial Higgs field conditions, see bulk text. The points indicate the average value m φ ti, while the envelope of each of the lines indicate the standard deviation σ ≡ N −1/2 i (xi −xi) 2 . For data points with ξ ≈ ξc, only some of the ten simulations do not become unstable, and hence we do not show the deviation in these cases. For ξ ξ (1) c all simulations are always stable (i.e. m φ ti = ∞), and hence data points are not drawn. g 2 6.9 · 10 −3 , so that t br t res for the values of ξ considered here. Therefore, in this case, our bounds can be applied.

A. Dependence of lattice simulations on the Higgs number of components and initial conditions
We address now how our results depend on the position of the momenta cutoff in the spectra of initial conditions, as well as on the number of Higgs components we put in our simulations.

Dependence on Higgs initial conditions
We have explained previously how the initial conditions of the Higgs field are set throughout the lattice. Basically, we impose at initial time t = t * vanishing homogeneous modes h n = 0 (n = 1, 2, 3, 4), and then we add quantum fluctuations to each of the components. These fluctuations are imposed only up to a certain cutoff momentum k c , so that for k > k c the fluctuations are set to zero. Also, the random nature of the initial conditions is implemented in the code through a pseudo-random number generator, so that different seeds produce different realizations for the initial conditions.
It is essential to fix the initial cutoff appropriately, so that the non-excited UV quantum modes, which cannot be treated in the lattice, are not excited as classical modes. In the results presented in Fig. 8, we have done several simulations with different initializations for each point. More specifically, for values ξ < 5, we have done ten simulations, five of them with κ c = 10 (κ c ≡ k c /m φ ), and the other five with κ c = 12. We have also varied the seed in each of the ten simulations. This matches quite well the analytical estimation for the classical estimation of modes during tachyonic resonance given in Eq. (7). For values 5 < ξ < ξ (2) c , the second set of five simulations has been done instead with cutoff κ c = 15. Finally, for points ξ > ξ (2) c , we have only done four simulations (two with κ c = 10 and two with κ c = 15), because for these points the dependence of our results on the initial conditions is negligible.
The top panel of Fig. 9 shows how the instability curves change for different choices of the initial cutoff κ c for the particular choice m t = 173.34 GeV. The inclusion of UV modes in the lattice beyond the physical cut-off, makes larger the Higgs amplitude h 2 , so that the negative λ h 2 /a 3 term in Eq. (12) is enhanced, and hence reduces the instability time m φ t i . For ξ 10, this effect is negligible, because as we saw in Section IV, the amplitude of the excited IR modes dominates over the UV ones, but it becomes increasingly important as ξ diminishes. As we decrease the coupling, the UV modes become more relevant, and if they are not appropriately eliminated, their contribution can make the instability time wrongly smaller. At very low couplings, this is related to the invalidity of the lattice approach, as explained in the last paragraph of Section IV.

Dependence on Higgs number of components
We now compare our results, in which we have taken the Higgs as a 4-component field (N c = 4), with a similar set of lattice simulations with a 1-component field N c = 1.
We expect differences between the two scenarios for several reasons. The first one is that, if we include a 4-component field, the tachyonic mass is exciting 4 scalar fields instead of one. If we neglect at first the Higgs self-interaction term, this means ϕ 2 n ≈ ϕ 2 /4 (n = 1, 2, 3, 4). Due to this, if we consider only a 1component Higgs, the magnitude of the negative selfinteraction term is being underestimated, and increases artificially the instability time m φ t i for a given coupling ξ, as well as the critical value ξ c . Also, in the 1-component case the Higgs potential possesses a discrete symmetry instead of a continuous one. This may change the dynamics of the transition to the negativeenergy vacuum.
To check this, we show in the bottom panel of Fig. 9 the dependence of the instability curve on the number of components, for the top quark masses m t = 172.12 GeV and m t = 173.34 GeV. We compare the cases N c = 4 (i.e. the case we have presented above), and N c = 1. As expected, for the 1-component case the critical coupling ξ c increases slightly. For the m t = 173.34 GeV case, we have ξ c ≈ 6 instead of ξ c ≈ 4, while for the m t = 172.12 GeV we have ξ c ≈ 13 instead of ξ c ≈ 12. Apart from that, we see that the particular shape of the instability curve is significantly changed, meaning that the effect of the interaction between the different Higgs components is relevant for the dynamics of the system.

VI. SIMULATIONS WITH GAUGE FIELDS
Until now, we have ignored the coupling of the Higgs field to the gauge bosons of the Standard Model. We now evaluate if the effects of this interaction modify significantly the results presented in the last Section.
Let us consider the following action in the continuum, which imitates the interactions of the Higgs with the four SU (2)×U (1) weak and hypercharge gauge bosons (W 1,2,3 and the Y ), in the Abelian approximation, Here . We take the SU(2) and U(1) couplings as constants with g 2 ≈ 0.3, and g 2 ≈ 0.3 (corresponding to their value at very high energies, according to the SM renormalization group). This action describes correctly the Higgs-gauge fields interactions of the SM, as long as the non-linear interactions of the gauge fields among themselves (due to the truly non-Abelian nature of the SM symmetries) can be ignored. This is typically a good approximation as long as the gauge fields are not largely excited.
As observed in [25], a system with N > 1 Abelian gauge-bosons can be equivalently described by a single effective gauge boson A µ . If we define the gauge boson field amplitudes in terms of the effective gauge boson as follows We plot the spectra of the conformal Higgs field for different times (Section VI). Continuous lines correspond to a Higgs coupled to gauge bosons (g 2 s = 1.2), while the dashed lines indicate the equivalent when such coupling is set to zero (g 2 s = 0). Here, we have chosen mt = 173.34GeV and ξ = 8. From early (red) to late times (purple), we have m φ (t − t * ) = 0, 2, 4, 8, 18, 59, 100, 161, 403. and substitute this in action (28), we recover action (6) for an effective gauge boson with gauge coupling g 2 s ≈ g 2 + 3g 2 . The effective gauge boson is then simply the sum of the original ones, i.e. A µ = Y µ + 3 a=1 W a µ . Naturally, what our lattice simulations do is to solve a discrete version of Eqs. (16), which we provide in Eq. (38) of the Appendix. Details of how we derive this equations and the assumptions we made are provided in more detail there. The results we present in this section are based on lattice simulations in cubes of N 3 = 128 3 points, with a minimum infrared momenta p min = 0.5m φ . This captures quite well the relevant range of momenta excited during the tachyonic resonance regime, for both the Higgs and the gauge fields.
Let us try to quantify the energy transferred from the  Higgs into the electroweak gauge bosons. Action (6) can be written as S = S m + S R , with S R ≡ d 4 x √ −gξR|Φ| 2 containing the Ricci-Higgs interaction term, and S m containing the other terms. We define the matter stressenergy tensor as T We show in Fig. 10 the evolution of the different contributions to the energy density (30) as a function of time, for the case ξ = 8. These energies have been divided by the inflaton energy ∼ 3 . We see that the Higgs and gauge fields energy is several orders of magnitude lower than the inflaton energy, which justifies neglecting their contribution to the Friedmann equation, as commented in Section III. At late times, the Higgs kinetic and gradient energies evolve as E ϕ K , E ϕ G ∼ a −4 , and thus eventually become sub-dominant with respect the magnetic energy.
We show in Fig. 11 the time-evolution of the Higgs spectra in the presence of a gauge interaction with g 2 s = 1.2, and compare it when such interaction is not present (g 2 s = 0). We clearly see that the gauge bosons have a very important backreaction effect on the Higgs field, propagating its spectra to the UV.
Finally, Fig. 12 shows the instability time m φ t i as a function of ξ obtained from lattice simulations, when we do include the coupling of the Higgs with the gauge bosons. We have simulated the cases m t = 172.12, 172.73, 173, 34GeV, and compared with the results obtained in Section V, when we ignored such coupling. Although the instability curves are slightly different with respect the case without gauge bosons, the values for the critical coupling ξ  Table III. In conclusion, our simulations demonstrate that the addition of gauge fields does not impact significantly in the post-inflationary dynamics of the system. The interaction of the Higgs with the electroweak gauge fields only changes marginally the results on the critical couplings ξ also indicates that the addition of the truly non-Abelian gauge bosons will not change the above conclusion, as the non-linear nature of the non-Abelian gauge field interactions cannot stimulate further the gauge bosons. Quite on the contrary, the non-linear structure of non-Abelian interactions typically prevents the stimulation of the gauge fields up to the level of excitation that (linear) Abelian interactions allow for.

VII. SUMMARY AND DISCUSSION
In this work, we have studied the post-inflationary dynamics of the Standard Model Higgs with lattice simulations, in the case where it possesses a non-minimal coupling ξ to gravity. This term is necessary for the renormalization of the theory in curved spacetime. We have assumed a chaotic inflation model with m 2 φ φ 2 potential. We include the running of λ(ϕ) in our simulations as a function of the value of the Higgs field at the lattice point. We have considered different runnings, corresponding to different experimental values of the topquark mass. The running is such that it generates two vacua to the Higgs potential: one at ϕ ≈ 0, and one at high-energies. With our lattice simulations, we have been able to obtain the critical coupling ξ c such that for ξ ξ c the Higgs field becomes unstable and decays into the negative-energy Planck-scale vacuum. Our lattice simulations also take into account the 4 components of the Higgs field and the cutoff of the spectra of initial fluctuations, which are necessary to correctly quantify the value of ξ c . We have done two sets of lattice simulations; one with only the Higgs field, including the effective expansion caused by the post-inflationary dynamics of the inflaton; and another in which we also include the coupling of the Higgs to gauge bosons (mimicked with an Abelian-Higgs-like model). We have observed that the effect of the gauge bosons can be very important in the Higgs post-inflationary dynamics. Tables II and III, together with the estimation ξ 0.06 from the stability of the Higgs field during inflation [7], provide tight constraints to the values of this coupling compatible with observations. However, we have assumed a chaotic inflationary model with potential m 2 φ φ 2 . It is expected that inflationary models with lower inflaton amplitudes during preheating will widen this range of values, as the value of the Ricci scalar |R(t)| decreases, and hence the excitation of the Higgs field due to the tachyonic resonance is less strong. If the Standard Model potential does not have a second negative-energy vacuum at high energies, we cannot find upper bounds for ξ c in this way.

The upper bounds in
In this work, we have neglected the terms coming purely from the non-Abelian structure of the SM Lagrangian, considering instead that the {W a , Y } bosons can be regarded as Abelian gauge fields. We have argued that considering linear Abelian interactions leads to a larger excitation of the gauge fields, so that the nonabelian terms can be safely ignored. We reach the important conclusion that the inclusion of gauge bosons in the system (even in the Abelian approach) does not change significantly the upper bound for ξ. The critical values  We present the discrete equations of motion we solve in our lattice simulations with gauge bosons, which as explained, include the four components of the Higgs field, and the three components of the effective gauge boson A µ , Eq. (29) (remember we set A 0 = 0). Let us define new spacetime variables as (z, z) ≡ m φ (t, x). Continuous action (6) can be written in terms of these variables as where we have defined R as Let us consider a lattice cube with N 3 points, with time step dt and lattice spacing dx. For a given field f , we define the following discrete forward and backward derivatives at a position n in the lattice as [dx 0 ≡ dt, dx i ≡ dx, (i = 1, 2, 3)] withσ µ a set of orthonormal vectors. Let us also define discrete lattice covariant derivatives as where V µ ≡ e −iAµdxµ is the link. We can write an action in the discrete equivalent to Eq. (32) as where the different pieces of the Lagrangian are Here, we take the scale factor as evaluated at semi-integer times, and we define a ≡ a +0/2 + a −0/2 2 .
We can minimize this Lagrangian with respect ϕ † and A i respectively, to obtain the discrete equations of motion. Setting A 0 = 0, they are On the other hand, minimizing the action with respect A 0 , we obtain the Gauss conservation law in the discrete, (39) In order to trust the results from our lattice simulations, the fields must preserve then the condition We have checked in our simulations that this is indeed the case, for all running time. Typically, ∆ G 10 −12 .
The mininum and maximum momenta captured by the lattice are p min = 2π L and p max = √ 3N 2 p min , with L ≡ dxN . We must choose N and L so that we capture the relevant range of momenta of this system. Fig. 13 shows the instability curve for the top quark mass m t = 173.34GeV, for different lattice parameters p min and N . We see that the results are consistent, independently on the particular features of the lattice.