Postinflationary vacuum instability and Higgs-inflaton couplings

The Higgs-inflaton coupling plays an important role in the Higgs field dynamics in the early Universe. Even a tiny coupling generated at loop level can have a dramatic effect on the fate of the electroweak vacuum. Such Higgs-inflaton interaction is present both at the trilinear and quartic levels in realistic reheating models. In this work, we examine the Higgs dynamics during the preheating epoch, focusing on the effects of the parametric and tachyonic resonances. We use lattice simulations and other numerical tools in our studies. We find that the resonances can induce large fluctuations of the Higgs field which destabilize the electroweak vacuum. Our considerations thus provide an upper bound on quartic and trilinear interactions between the Higgs and the inflaton. We conclude that there exists a favourable range of the couplings within which the Higgs field is stabilized during both inflation and preheating epochs.


Introduction
The discovery of the Higgs boson at the LHC [1,2] in July 2012 furnished the final piece of the Standard Model (SM) of particle physics; however, it has also raised important new questions. One of these relates to the issue of the electroweak vacuum stability and the fate of the Higgs field in the early Universe, particularly during the inflationary and reheating eras [3].
For the currently preferred values of the top quark mass and the strong coupling, the self-coupling of the Higgs field turns negative at a high energy scale of order µ c ∼ 10 10 GeV [4][5][6] (see [7] for its gauge (in)dependence). This would suggest that there exists another, deeper vacuum state than the one we currently occupy. One finds then that the electroweak vacuum is metastable with the lifetime longer than the age of the Universe. Although this does not pose an immediate problem, the existence of the deeper vacuum raises cosmological questions. In particular, one must explain how the Universe ended up in an energetically disfavored state and why it stayed there during inflation [8]. Even if one fine-tunes the Higgs field initial conditions before inflation, light scalar fields experience large fluctuations of order the Hubble rate H during the exponential expansion epoch [9]. Unless H is sufficiently small, the Universe is overwhelmingly likely to end up in the catastrophic vacuum [10].
These problems can be solved by coupling the Higgs field to the scalar curvature [3] or by taking into account the Higgs-inflaton coupling [8]. We focus on the latter possibility in this paper and neglect the effect of the non-minimal coupling to gravity. 1 As shown in [11], Higgsinflaton interaction is inevitable in realistic models of reheating. Indeed, the inflaton energy must be transferred to the Standard Model fields which leads to a (perhaps indirect) coupling between the inflaton and the SM particles. The latter induces Higgs-inflaton interaction at loop level, where H is the Higgs doublet and φ is a (real) inflaton. Here λ hφ and σ hφ typically receive log-divergent loop contributions and thus require renormalization. In other words, these couplings are generated by the renormalization group (RG) evolution [11]. Their magnitude can be large enough to alter the Higgs evolution completely, in particular, by inducing a large effective Higgs mass which drives the Higgs field to zero. This mechanism is operative in the range 10 −10 < λ hφ < 10 −6 , (1. 2) with the upper bound coming from the requirement that the Higgs-inflaton interaction preserve flatness of the inflaton potential, and the lower limit dictated by the condition that the Higgs effective mass be greater than the Hubble rate during inflation. The trilinear interaction should be subdominant, λ hφ φ 2 σ hφ φ so that the effective mass term does not depend on the sign of the inflaton field. This is usually the case in explicit reheating models [11].
In this work, we study the effect of the above couplings after inflation. Although the Higgs-inflaton interaction can stabilize the Higgs potential during inflation, during preheating its effect can instead be destabilizing (see also [12]). The parametric resonance [13,14] due to the quartic interaction h 2 φ 2 and the tachyonic resonance [15,16] due to the h 2 φ term can lead to very efficient Higgs production. This causes large fluctuations and the Higgs variance h 2 that can exceed the critical value beyond which the system becomes unstable. We find that these considerations place important upper bounds on both λ hφ and σ hφ such that the range of favored couplings (1.2) reduces.
The field of Higgs dynamics in the early Universe has been very active in the recent years. Higgs field fluctuations during inflation in the metastable Universe have been studied in [17,18] and [19][20][21]. The Higgs condensate dynamics assuming stability of the Higgs vacuum were analyzed in detail in [22,23]. These considerations are affected by the presence of further Higgs interactions which are usually not included in the Standard Model. The effect of the non-minimal coupling Higgs to gravity on the Higgs dynamics was recently refined in [24]. In this framework, it was also noted that the resonances during preheating can destabilize the electroweak vacuum [12]. The effect of the quartic Higgs-inflaton coupling on the Higgs production during preheating was considered in detail in [25] (see also [26,27]). Our present work goes beyond these previous studies in that we consider a more realistic case of both quartic and trilinear interactions present, which brings in new and important qualitative features. We also refine the earlier analysis of the pure quartic case. Finally, we discuss implications of our findings for realistic reheating models. This paper is organized as follows. In the next section, we present our setup. In section 3, we consider the effect of the quartic Higgs-inflaton interaction on Higgs production during preheating. Section 4 is devoted to the more realistic case of both trilinear and quartic interactions present.

Framework
In this section, we present our inflationary setup. For concreteness, we study the Higgs production within the simple m 2 φ 2 chaotic inflation model with while our results easily generalize to other large field models. In the unitary gauge H = (0, h/ √ 2) T , the relevant Lagrangian is given by where the self-coupling λ h (h) is determined by the RG equations of the Standard Model. During inflation, φ undergoes a slow-roll evolution. On the other hand, for λ hφ > 10 −10 and a sufficiently large initial inflaton value φ 0 M pl , the Higgs mass is dominated by the inflaton interaction, m eff h λ hφ /2|φ| [8]. 2 Then Higgs field evolves exponentially quickly to zero.
Not long after the end of inflation, the inflaton field undergoes oscillations As long as the energy density of the Universe is dominated by the inflaton oscillations, the scale factor behaves as a = (t/t 0 ) 2/3 and the amplitude of oscillations decays as For concreteness, we assume that φ(t) in eq. (2.3) becomes a good approximation to the evolution of the inflaton at Φ 0 0.2M Pl . (2.5) Soon thereafter we can accurately approximate the time dependence of Φ(t) as The inflaton induced Higgs mass term also oscillates which can lead to efficient Higgs production. As one can see in eq. (1.1) the second term grows with respect to the first one as Φ(t) decreases due to the expansion of the Universe. Therefore the effect of the trilinear term becomes important at some stage even though it was negligible during inflation. Since the consequent effective Higgs mass term can have either sign ∝ σ hφ φ, the tachyonic resonance becomes effective. Both of the resonances play an important role and will be studied in the next sections.
Before we proceed, let us clarify our assumption about the running coupling λ h (µ), where µ is the renormalization scale. During the resonances, the Higgs quanta are produced coherently with the corresponding occupation numbers being very large. Thus we may treat h semi-classically. In this regime, we may take where h 2 plays the role of the relevant energy scale at which the coupling should be evaluated. Since we are only interested in the high energy regime, in our numerical analysis we use the step-function approximation where h SM c ∼ 10 10 GeV is the critical scale of the Standard Model at which λ h flips sign.

Pure Parametric Resonance
Let us first consider the case where the trilinear interaction is negligible, σ hφ ≈ 0. The Higgs-inflaton interaction is quartic so that we recover the well-known parametric resonance setting [14]. The equations of motion for the Higgs field are quadratic in h apart from the quartic self-interaction. During the parametric resonance regime, the effect of the latter can be approximated as h 4 → 6h 2 h 2 , which is known as the Hartree approximation. In that case, the equations of motion for different momentum modes decouple. In terms of the rescaled Higgs momentum modes X k ≡ a 3/2 h k , where a is the scale factor, one has [14] In the last term, w = p/ρ = − 1 + 2 3Ḣ H 2 is the equation of state parameter of the Universe, which vanishes in the matter-like background. We thus neglect this term.
If the Higgs-inflaton coupling λ hφ is substantial, the Higgs modes experience amplification due to broad parametric resonance. The parameter characterizing the strength of the resonance is such that q 1 corresponds to the broad resonance regime. In this case, the modes grow exponentially leading to a large Higgs field variance h 2 . The fluctuations can be so significant that they exceed the size of the barrier separating the electroweak vacuum from the catastrophic one at large field values. In this case, vacuum destabilization occurs. In what follows, we will estimate the corresponding critical size of λ hφ .
As was shown in [14], in the broad resonance regime the Higgs modes evolve adiabatically away from the inflaton zero-crossings and can be described by the WKB approximation where α k , β k are some constants. Adiabaticity is broken for certain modes near the inflaton zero-crossing, where the frequency ω k evolves very quickly. There the system can be treated in analogy to the Schrödinger equation as a scattering of plane wave solutions. The adiabatic constants α k and β k can be thought of as Bogolyubov coefficients. We assume a vacuum initial condition for the Higgs modes with α k = 1 and β k = 0. The occupation number of Higgs quanta after j mt/π zero crossings is then and can be written in terms of the corresponding Floquet index µ j k as [14] n j+1 k e 2πµ j k n j k .
(3.5) µ j k can be calculated via scattering of plane waves in a parabolic potential [14], Here a j is the scale factor after j zero crossings. The term sin θ tot is determined from the phase accrued by the modes and behaves in a stochastic manner for different momenta (see [14]). We take it to be zero for our estimates and use the consequent average value of the Floquet index. The occupation numbers at late times (a j 1) can then be approximated by Here the factor of 1/2 is due to the vacuum fluctuations in the initial state, although it plays no tangible role in our analysis.
Using the saddle point approximation, the Higgs field variance can be written as where κ 2 max =μ −1 j is the momentum in units of m which contributes most significantly. Here we have assumed that ω k is dominated by the inflaton-induced term. 3 Already after the first zero-crossing h 2 exceeds the critical scale ∼ 10 10 GeV of the Standard Model and therefore the Higgs self-coupling λ h = λ h h 2 can be taken to be negative from the beginning.
For our analysis, we take λ h = −10 −2 at large field values. Note that the fact that h 2 exceeds the SM critical scale does not necessarily lead to vacuum destabilization since the presence of the Higgs-inflaton coupling pushes the barrier separating the two vacua to larger values of order However, the position of the barrier is modulated by | cos mt| so it is not immediately clear what vacuum stability would require.
To derive the stability condition, one can use the following reasoning. Around each inflaton zero crossing, the effective Higgs mass squared is dominated by the Higgs self-interaction term λ h h 2 . Such a tachyonic term leads to exponential amplification of the Higgs field by a factor of order e m eff h ∆t , where m eff h is the modulus of the effective Higgs mass term and ∆t is the (short) period during which the Higgs self-interaction term dominates. ∆t is given explicitly by The tachyonic amplification is insignificant as long as m eff h ∆t does not exceed unity, that is, Clearly, this condition eventually gets violated since ∆t grows as the inflaton amplitude Φ decreases. However, if the resonance ends before this takes place, no destabilization occurs. Using λ hφ Φ 2 2m 2 at the end of the resonance [14], one finds that the stability condition can be written as Here we have neglected a smaller additive constant in the square brackets. If this condition is violated, the Higgs field grows explosively since the amplification factor e m eff h ∆t increases with h 2 itself. This leads to fast vacuum destabilization. On the other hand, if this condition is satisfied, it implies that the Higgs potential is dominated by the inflaton coupling term on the average and h does not fluctuate beyond the barrier (3.9). This result is consistent with the bound obtained in [25]. Figure 1 shows our numerical evolution of the Higgs fluctuations for different values of λ hφ . To produce this plot we have solved the mode equations in the Hartree approximation using Mathematica software. We see that for λ hφ greater than a few times 10 −8 , the Higgs field grows above the critical value and blows up at mt > 40. The destabilization time however should not be taken at face value since the Hartree approximation turns out to be rather crude for this purpose.
We have also performed a more sophisticated lattice simulation which takes into account the Higgs self-interaction without resorting to the Hartree approximation. We used the LATTICEEASY package [28] for this purpose. In our simulations, we choose the box size of 10/m (the L parameter of LATTICEEASY) with 64 grid points per edge (the N parameter). We have checked that a larger and finer grid does not change the results significantly. In Figure 2, we plot the destabilization time versus λ hφ . The green curve shows mt at which the system is destabilized, that is, the Higgs field variance blows up. The red line marks the end of the resonance such that the points below it correspond to vacuum destabilization during the resonance as studied in this section. 4 Our theoretical bound on λ hφ is marked by the vertical dashed line. We see that the latter describes the general situation reasonably well and λ hφ above 3 × 10 −8 typically leads to vacuum destabilization during the resonance. On the other hand, we also see the limitations of our approach. In particular, Figure 2 shows that the strength of the resonance does not behave monotonically with λ hφ . This is expected since we have taken the term sin θ tot to be zero, whereas in reality it either enhances or suppresses the resonance such that there can be certain values of λ hφ satisfying our bound yet leading to an unstable configuration. Formally, the area around λ hφ ∼ 2 × 10 −8 appears to be stable during the resonance and the destabilization occurs shortly after the resonance. However, one can classify this region as unstable since in reality the end of the resonance is not sharply defined due to various approximations we have made. Apart from these complications, we find that our simple considerations give a fairly good description of the system behavior during the parametric resonance.
Comparing Fig. 1 and Fig. 2, one finds that the commonly used Hartree approximation overestimates the destabilization time. This is to be expected since the quantity h 4 experiences greater fluctuations than h 2 h 2 does. Nevertheless certain questions such as the effect of perturbative Higgs are easier addressed using our Mathematica routine which employs the Hartree approximation. Hence we use both numerical approaches.  So far our discussion has ignored perturbative decay of the Higgs quanta, which reduces the efficiency of the resonance and can potentially invalidate our conclusions. The main decay channel is provided by the top quarks which are effectively massless for our purposes. The corresponding decay width is Taking y t (m eff h ) ∼ 1/2, m eff h λ hφ /2 |φ| and averaging | cos mt|, we find that the perturbative decay reduces the number of the Higgs quanta by a factor 2 or so in the region of interest (λ hφ ∼ 10 −8 ), see the left panel of Figure 3. Therefore it does not significantly affect our bound on λ hφ . On the other hand, for larger λ hφ ∼ 10 −7 , the Higgs decay can reduce h 2 by an order of magnitude thus delaying (but not avoiding) vacuum destabilization.
For completeness, in the right panel of Figure 3, we present a typical example of the occupation number evolution for different momenta. We find that at late times the Higgs field is typically dominated by the modes with momenta k ∼ m.

Effect of the Trilinear Interaction
The trilinear Higgs-inflaton interaction brings in an additional effective mass term whose sign oscillates in time. This results in the tachyonic resonance [15] which amplifies the Higgs fluctuations. We find that the effect is important and cannot be neglected.
The parametric and tachyonic resonances have been studied separately in detail. In realistic models, both of them are present at the same time, yet their combined effect is not well understood (see however [16], [29]) 5 . In particular, the Higgs field goes through a sequence of exponential amplification periods and plateaus. In what follows, we study some of the important aspects of the system and obtain the corresponding bound on σ hφ .

Equations of Motion
The trilinear interaction introduces an additional oscillating contribution to the effective Higgs mass. The frequency of this contribution is half the frequency of the quartic interaction. In particular, the Higgs dispersion relation in eq. (3.1) becomes where, as in eq. (3.1), we have neglected terms proportional toḢ ∼ H 2 . These terms become small, as compared to m 2 , soon after the end of inflation. Let us introduce and the q(t) parameter, which is defined in eq. (3.2). Then the equation of motion for the (rescaled) Higgs field can be written as where ρ 1 and ρ 2 are integration constants, y (z) are periodic functions of period π and µ is a characteristic exponent, or Floquet exponent, which in general is a complex number. When µ attains a real part, the solution grows exponentially. We discuss the most important properties of these solutions in Appendix A. In particular, the stability chart of the Whittaker-Hill equation is quite different from that of the Mathieu equation in the parameter range of interest (see fig. 9). In reality, the Universe expansion cannot be neglected and leads to the end of the resonance. Hence the Whittaker-Hill equation only provides a simple approximation to the equations of motion. The duration of the resonance is essential for our considerations since it determines the size of h 2 . Let us consider it in detail.

Duration of the Resonance
The essential difference between the solutions of eq. (4.4) in an expanding and static Universes is that in the former case the boundaries between the stability and instability regions are no longer clearly defined: they are smeared [14]. Despite this fact, we will use figure 9 as a helpful illustration. For that matter, the time dependence of A (k, z), p (z) and q (z) can be introduced adiabatically: as they evolve, one can think of them as tracing a trajectory in the three dimensional space, crossing through stable and unstable regions. Once these parameters decrease substantially and the trajectory converges to the lowest stable region, the resonance ceases.
The definition of the parameter A in eq. (4.6) contains two terms, which are time dependent. To compute the duration of the resonance, we first estimate the relative size of these two contributions at the end of the resonance. Using eq. (2.4) and the definition of q in eq. (3.2), we can write where the subscript f refers to values at the end of the resonance. Typically, the excited modes towards the end of the resonance are k/m ∼ 1. Since m/Φ 0 ∼ 5 × 10 −6 , we have In this work we are interested in the range of values of λ hφ given in equations (1.2) and (3.12).
Taking also q 1/3 f ∼ 1, the above ratio lies in the range 10 −2 ...1. That is, at the end of the resonance, the q-term dominates and it suffices for our purposes to consider the evolution of the k = 0 mode only. This restricts our parameter space to the plane A = 2q. The stability and instability regions for constant p,q can be obtained by the methods discussed in Appendix A. The result is shown in figure 4, where the labeled curves display the trajectories p(t), q(t) for different λ hφ and σ hφ . The vertical line p = 0 corresponds to the parametric resonance and one recovers the standard results of ref. [14].
The resonance stops when q (z) and p (z) reach the last stable region around p = 0 and q = 0 in figure 4. To estimate the time when this happens, we approximate the boundary of the lowest stable region by a linear relation q = 0.48 − 0.53 |p|. This approximation is shown by bold red lines in figure 4. One has to keep in mind however, that in the expanding universe the boundaries between stable and unstable regions are smeared. Hence, even when a trajectory in (q, p) parameter space reaches the last stable region, the resonance continues for some time, depending on the phase. Thus, the end of the resonance corresponds to where δ is a "fudge" factor to be determined from simulations. Our results show that δ varies from 0 to about 1/4. An analogous result for the parametric resonance was obtained in [14], in which case the resonance stops somewhere in the range 1 ≤ q f ≤ 4/3. 6 In our parameter range, we find that δ 0.1 gives a good approximation for most cases.  In figures 5 and 6, we plot numerical LATTICEEASY computations of occupation numbers n k with k ≈ 0.63m and h 2 for several models. 7 The values of t f from eq. (4.11) are shown by dashed vertical lines. We conclude that the agreement is quite good. Eq. (4.11) does not apply for very small values of σ hφ such that the tachyonic resonance is inefficient. In particular, the amount of time the system spends in the last instability region (just above the red line in fig. 4) is so small that no substantial amplification occurs. For such models, the dynamics of the resonance are close to those of the pure parametric case [14]. It is also worth recalling that λ hφ in eq. (4.11) is not allowed to be too small so that the ratio in eq. (4.9) is below unity.

Vacuum Destabilization by a Mixed Resonance
As in the parametric resonance case, the Higgs field fluctuations can grow large enough so that the system moves over to the catastrophic vacuum. This transition is facilitated by the presence of the trilinear term which results in very large Higgs occupation numbers. In what follows, we study the destabilization effect due to σ hφ . That is, we choose λ hφ for which the system is stable and analyze how large a σ hφ one can add without destabilizing the vacuum. As before, we focus on the destabilization during the resonance, i.e. before t f in eq. (4.11).

Simplified Bound on σ hφ
The analysis of the mixed trilinear-quartic case is substantially more complicated than the pure quartic case. As seen from the stability chart, the system goes through a series of stable and unstable regions with a varying exponent µ(t). We will thus content ourselves with only an order of magnitude estimate of the critical σ hφ .
Towards the end of the resonance, the Higgs-dependent potential is dominated by the trilinear term 1 2 σ hφ φh 2 since the quartic interaction decreases faster with time. The destabilization occurs when this term becomes overtaken by the Higgs self-interaction 1 4 Therefore, one can estimate the critical variance by where the Hartree approximation has been used and the oscillatory behavior of φ has been ignored.
On the other hand, the Higgs variance as a function of time can be calculated via the occupation numbers as in (3.8). The dominant contribution is given by modes around the comoving momentum k * which maximizes n k . For the parameter range of interest, we find that k * ∼ m towards the end of the resonance and the width of the k-distribution is of order k * /2. The corresponding n k * is a rather complicated function of time containing sections where it undergoes an exponential increase. For our purposes, we simply interpolate it by e µ * mt with some effective exponent µ * . We then obtain The destabilization occurs if h 2 reaches the critical value during the resonance. The latter stops around 2|σ hφ |Φ end m 2 . Taking this into account and dropping order one constants, one finds |λ h | ∼ 10 9 GeV , (4.14) with a end = (Φ 0 /Φ end ) 2/3 being the scale factor at the end of the resonance. Here we take a typical value 8 µ * ∼ O(10 −1 ) (cf. fig. 9). Note that the main σ hφ -dependence of the result comes from the duration of the resonance, mt end ∼ σ hφ × M pl /m 2 , while that of µ * and ln a end is milder. Although this estimate is very crude, we find that the bound is within a factor of a few from our numerical results. Here we have neglected both the λ hφ -dependence and the dependence on the sign of σ hφ .
Note that both the λ hφ and σ hφ bounds do not appear to depend explicitly on the critical scale of the Standard Model. This dependence is hidden in our assumption about λ h at the energy scales of interest. As long as λ h ∼ −10 −2 in that range, our bounds apply.
Finally, we have considered a chaotic φ 2 inflation model which fixes a large H infl ∼ 10 14 GeV. For models with small H infl < 10 10 GeV, the Higgs fluctuations during inflation are not dangerous and the Higgs-inflaton coupling can be set negligibly small. 9

Simulation Results
Our LATTICEEASY simulations show that the bound on σ hφ depends both on λ hφ and the sign of σ hφ . The latter is due to the fact that even though the Whittaker-Hill equation enjoys the symmetry z → z + π/2, p → −p, the time translation invariance is broken by the Universe expansion. Figs. 7 and 8 display the bounds on σ hφ as a function of λ hφ . We see that the upper bound varies between 10 8 GeV and 6 × 10 8 GeV in the region of interest.
We should note that these plots are somewhat simplified in that it is tacitly implied that |σ hφ | below the critical value leads to a stable system. In practice, this is not always the case and the destabilization time can be a non-monotonic function of σ hφ . However, these effects do not change our results drastically.

Comments on the Late Time Behavior
So far we have discussed vacuum destabilization during the resonance. The initial stage of preheating is dominated by a single process, that is, resonant Higgs production. At later stages, other processes such as rescattering, thermalization, etc. become important.
As seen in fig. 2, the Higgs vacuum can be destabilized much after the end of the resonance. The simple reason for it is that h 2 and the position of the barrier h c ∝ Φ scale differently in time. If only the quartic coupling is present, where α is between 1 and 3/4, depending on which k-modes dominate h 2 . This can be seen from the first equality in (3.8) and the fact that the comoving occupation numbers are constant after the end of the resonance, while the ω k scaling depends on the balance between k 2 /a 2 and the inflaton-induced mass term. In any case, h 2 decreases slower in time than h c does so that after a sufficiently long period the Higgs fluctuations go over the barrier. Fig. 2 shows that the relevant time scale is of order 100mt. Analogous considerations also apply to the mixed trilinear-quartic case.
However, the true dynamics of the system on a larger time scale are complicated. The Higgs interacts with other fields of the Standard Model which becomes important after the resonance. As noted in [25], thermalization effects can generate a thermal mass term for the Higgs thereby stabilizing the vacuum. Also, non-perturbative production of particles via the Higgs couplings can reduce h 2 [23]. These effects are subtle and require a careful investigation which is beyond the scope of our present work. On the other hand, the resonance regime is quite well understood and thus we believe our bounds on λ hφ and σ hφ are solid.

Implications for Reheating Models
In this section, we consider implications of our bounds for model parameters of representative reheating scenarios. We choose two examples considered in [11]: reheating via right-handed neutrinos and reheating via non-renormalizable operators.
In general one expects the Higgs-inflaton couplings to be present already at the tree level. However, if they are for some reason suppressed, λ hφ and σ hφ are generated by loop corrections. Therefore, the loop-induced couplings can be regarded as the corresponding lower bound. In what follows, we consider two conservative scenarios in which λ hφ and σ hφ are entirely due to loop effects.

Reheating via Right-handed Neutrinos
In this model, the inflaton decays into heavy right-handed neutrinos which subsequently decay into SM particles. This option is attractive since the inflaton-neutrino coupling is allowed already at the renormalizable level. The relevant interaction terms are where l L is the lepton doublet, the Majorana mass M is chosen to be real and we have assumed that a single ν R species dominates. The quartic and trilinear Higgs-inflaton couplings are generated at 1 loop and the result is divergent. In other words, such couplings are required by renormalizability of the model. As the renormalization condition, we take λ hφ (M Pl ) = 0, σ hφ (M Pl ) = 0 such that at the inflationary scale the couplings are generated by loop effects.
In the leading-log approximation, we find where µ is the relevant energy scale. In what follows, we assume real couplings and take µ ∼ m since this is the typical momentum of the Higgs quanta towards the end of the resonance. 10 In any case, the dependence on µ is only logarithmic. The value of λ ν is constrained by inflationary dynamics. In order not to spoil flatness of the inflaton potential, the coupling must satisfy λ ν < 10 −3 [11]. Taking λ hφ < 3 × 10 −8 and |σ hφ | < 10 8 GeV (see fig. 7), we find the following bounds on the neutrino Yukawa coupling and the Majorana mass, Although these constraints are not particularly strong, they are non-trivial. In particular, they imply that the neutrino Yukawa coupling cannot be order one.

Reheating via Non-renormalizable Operators
A common approach to reheating is to assume the presence of non-renormalizable operators that couple the inflaton to the SM fields. Let us consider a representative example of the following operators, where Λ 1,2 are some scales, G µν is the gluon field strength and q L , t R are the third generation quarks. These couplings allow for a direct decay of the inflaton into the SM particles. It is again clear that a Higgs-inflaton interaction is induced radiatively. In order to calculate the 1-loop couplings reliably, one needs to complete the model in the ultraviolet (UV). The simplest possibility to obtain an effective dim-5 operator is to integrate out a heavy fermion. Therefore, we introduce vector-like quarks Q L , Q R with the tree level interactions where the heavy quarks have the quantum numbers of the right-handed top t R , their mass M is taken to be above the inflaton mass scale and the couplings to the third generation are assumed to dominate. One then finds that O 1 appears at tree level with 1/Λ 1 = y Q λ Q /M, whereas O 2 appears only at 2 loops with 1/Λ 2 ∼ y Q λ Q y t α s /(64π 3 M) and can be neglected.
Using the renormalization condition that the relevant couplings vanish at the Planck scale and the fact that the heavy quarks contribute only at scales above M, we get in the leadinglog approximation where y t is the top Yukawa coupling and we assume M M Pl . As in the previous example, one of the couplings is constrained by the inflationary dynamics, |λ Q | < 2 × 10 −3 /(ln M Pl /M) 1/4 [11], since it generates a correction to the inflaton potential. The heavy quark mass must be well below the Planck scale, M M Pl , and the bound on λ Q depends 10 Choosing a higher µ would result in slightly looser bounds.
very weakly on M in the allowed range. Therefore, in practice one may take |λ Q | < 10 −3 .
Our results λ hφ < 3 × 10 −8 , |σ hφ | < 10 8 GeV lead to a stronger bound. For real couplings, we get where in the second inequality we took M ∼ m to obtain the most conservative bound and y t (M) ∼ 1/2. This implies, in particular, that the minimal value of the suppression scale Λ 1 = M/|λ Q y Q | is around the Planck scale and the maximal reheating temperature is of order 10 9 GeV (see [11] for details).

Conclusions
This work is devoted to an in-depth analysis of the Higgs-inflaton coupling effects in the reheating epoch. We have focussed in particular on the preheating stage when the parametric and tachyonic resonances are active. Our framework includes both the quartic and trilinear couplings since these are present simultaneously in realistic models. The resulting mixed parametric-tachyonic resonance is described by the Whittaker-Hill equation. While inheriting certain features of the two resonances, it brings in new effects which require a thorough investigation. Within this framework, we have analyzed the issue of electroweak vacuum stability during the preheating epoch assuming that the Higgs self-coupling turns negative at high energies. Even though the Higgs-inflaton couplings can stabilize the system during inflation, resonant Higgs production thereafter can lead to vacuum destabilization. The relevant quartic and trilinear Higgs-inflaton couplings are generated by the renormalization group equations in realistic models, and even their tiny values make a difference. Using both analytical methods and lattice simulations in a representative large field (φ 2 ) inflation model, we obtain upper bounds on the couplings from vacuum stability during preheating. These allow for a range of couplings, roughly 10 −10 < λ hφ < 10 −8 and |σ hφ | < 10 8 GeV, which ensure stability both during inflation and preheating.
Our analysis is limited to the timescale of the mixed resonance. This leaves out the issues of the late-time behavior of the Higgs fluctuations which can further limit the allowed range for the couplings. The required analysis is highly involved and we leave it for future work.

A.1 Computation of the Floquet Exponent
The Whittaker-Hill equation is given by where a constant A can be thought of as an eigenvalue of the differential operator on the LHS. The analysis of this equation can be found in [29][30][31][32].
According to the Floquet theorem, the solution can be written as a series of the form Plugging this Ansatz into the above equation, we obtain a recursive relation For given values of A, q and p we can find the Floquet characteristic exponent µ by solving for the roots of the determinant It is possible to prove (see, e.g., refs [29,30]) that this determinant can be written in a compact form as sin 2 iµ π 2 = ∆ (0) sin 2 √ A π 2 . (A.6) From this equation we can easily find µ = − i π arccos 1 + ∆ (0) cos The advantage of this representation of solutions is that it can be evaluated numerically very efficiently. Indeed, as one can see from the definitions of γ 2n and ξ 2n in eqs. (A.4), the off-diagonal elements of ∆ (iµ) decrease as ∝ n −2 as they depart from the center of the matrix.

A.2 Boundary Between Stability and Instability Regions
The stability of the solution in eq. (4.7) is determined by the characteristic exponent µ. In general, µ is a complex number µ = α + iβ. If the real part of µ is non-zero, that is α = 0, the given solution is unstable. For stable solutions α = 0 and their periodicity is determined by the value of the imaginary part β. If β is a rational fraction, the solution is periodic, while for irrational β the solution is non-periodic. Particularly interesting are the cases where β is an integer. If β = 2l, where l ∈ Z, then solutions are either even or odd periodic functions with a period π. For β = 2l + 1 those solutions are even or odd periodic functions with a period 2π. These solutions of period π and 2π lie on the boundary between the regions where the families of stable and unstable solutions reside, that is, the so called stability and instability regions in (A, q, p) space. To find the equations for these boundary surfaces, we can use the following Ansätze The function y 1 (z) describes even solutions of period π; y 2 (z) is an odd function of period 2π; y 3 (z) is an even function of period 2π and y 4 (z) is an odd function of period π. To find the coefficients, one plugs these functions back into the Whittaker-Hill equation. This gives us four recursive relations which can be written in a matrix form, 12) where no summation over J is implied; C J stands for C 1 = (C 0 , C 2 , C 4 , . . .) T , C 2 = (S 1, S 3 , . . .) T , etc. and A J represent eigenvalues A IJ (p, q) of an infinite square matrix M J , where I = 0, 1, . . . , ∞. These eigenvalues define the boundaries between the stability and instability regions.
To find A IJ (p, q), let us compute the four matrices M J explicitly. Plugging eq. (A.8) into eq. (A.1), we find where the first row corresponds to n = 0. Similarly, one obtains the other three matrices, (A.14) We compute the eigenvalues of these matrices numerically by truncating them at some high value of n. Some solutions of eq. (A.7) and eigenvalues of these matrices are shown in figure 9. The leftmost panel (q = 0) can be recognized as the familiar stability chart of the Mathieu equation.