Massless Preheating and Electroweak Vacuum Metastability

Current measurements of Standard Model parameters suggest that the electroweak vacuum is metastable. This metastability has important cosmological implications, because large fluctuations in the Higgs field could trigger vacuum decay in the early universe. For the false vacuum to survive, interactions which stabilize the Higgs during inflation -- e.g., inflaton-Higgs interactions or non-minimal couplings to gravity -- are typically necessary. However, the post-inflationary preheating dynamics of these same interactions could also trigger vacuum decay, thereby recreating the problem we sought to avoid. This dynamics is often assumed catastrophic for models exhibiting scale invariance since these generically allow for unimpeded growth of fluctuations. In this paper, we examine the dynamics of such"massless preheating"scenarios and show that the competing threats to metastability can nonetheless be balanced to ensure viability. We find that fully accounting for both the backreaction from particle production and the effects of perturbative decays reveals a large number of disjoint"islands of (meta)stability"over the parameter space of couplings. Ultimately, the interplay among Higgs-stabilizing interactions plays a significant role, leading to a sequence of dynamical phases that effectively extend the metastable regions to large Higgs-curvature couplings.


I. INTRODUCTION
A remarkable implication of the currently measured Standard Model (SM) parameters is that the electroweak vacuum is metastable. Specifically, given the measured Higgs boson and top quark masses [1], one finds that at energy scales exceeding µ ≈ 10 10 GeV the Higgs fourpoint coupling runs to negative values λ h (µ) < 0, signifying the existence of a lower-energy vacuum [2][3][4][5][6]. Although today the timescale for vacuum decay is much longer than the age of the universe [7,8], dynamics earlier in the cosmological history could have significantly threatened destabilization. That the false vacuum has persisted until the present day may thus provide a window into early-universe dynamics involving the Higgs [9].
In this respect, the evolution of the Higgs field during inflation is especially relevant. During inflation, light scalar fields develop fluctuations proportional to the Hubble scale H. Without some additional stabilizing interactions, the fluctuations of the Higgs would likewise grow and inevitably trigger decay of the false vacuum, unless the energy scale of inflation is sufficiently small [9][10][11].
Metastability thus provides a strong motivation for investigating non-SM interactions that could stabilize the Higgs during inflation-e.g., non-minimal gravitational couplings [12], direct Higgs-inflaton interactions [13], etc. However, the situation is actually more delicate, as one must also ensure metastability throughout the remaining cosmological history. While interactions such as those listed above may stabilize the vacuum during inflation, they often proceed to destabilize it during the post-inflationary preheating epoch, thereby recreating the problem we sought to avoid. Indeed, a balancing is typically necessary between the destabilizing effects of inflationary and post-inflationary dynamics.
To this end, a detailed understanding of the field dynamics after inflation is essential in determining which interactions and couplings are well motivated overall. A decisive component of this analysis is the form of the inflaton potential V (φ) during preheating. After inflation, the inflaton field φ oscillates coherently about the minimum of its potential. These oscillations determine the properties of the background cosmology, but they also furnish the quantum fluctuations of the Higgs field with time-dependent, oscillatory effective masses. These modulations are the underlying mechanism for Higgs particle production, as they ultimately give rise to non-perturbative processes such as tachyonic instabilities [14][15][16] and parametric resonances [17][18][19]. These processes can enormously amplify the field fluctuations over the course of merely a few inflaton oscillations.
For most inflaton potentials-such as those which are quadratic after inflation-particles are produced within bands of comoving momentum, and these bands evolve non-trivially as the universe expands. The rates of particle production for these bands also evolve and generally weaken as preheating unfolds. As a rule, the particle production terminates after some relatively short time, and one can classify a model that remains metastable over the full duration as phenomenologically viable [20][21][22][23][24].
However, there is a notable exception to this rule: models which exhibit scale invariance. Under a minimal set of assumptions, in which the inflaton is the only non-SM field, the scalar potential during preheating is restricted to the following interactions: that even if the inflaton-Higgs interaction does not appear as a direct coupling, it should be generated radiatively since inflaton-SM interactions are generally necessary to reheat the universe [25]. Additionally, we emphasize that we have made no assumptions regarding the potential in the inflationary regime, except that it smoothly interpolates to Eq. (1.1) at the end of inflation. Above all, the scale invariance implies that the system's dynamical properties are independent of the cosmological expansion, in stark contrast to other preheating scenarios. As a result, the salient features for vacuum metastability, in general, do not evolve: particles are produced within momentum bands that are fixed with time, and the production rates for the modes do not change. The field fluctuations grow steadily and unimpeded [19].
At first glance, the unimpeded growth in massless preheating appears catastrophic for electroweak vacuum metastability. That said, several considerations should be evaluated more carefully before reaching this conclusion. First, the Higgs particles produced during preheating have an effective mass proportional to the inflaton background value m h g|φ|. This dependence can significantly enhance the perturbative decay rate of the Higgs to SM particles. 1 For large enough coupling, the decay rate and the production rate could be comparable, thereby checking the growth of Higgs fluctuations and effectively stabilizing the vacuum. Second, as particles are produced, their backreaction modifies the effective masses of the field fluctuations. If the vacuum does not decay first, the energy density of the produced particles inevitably grows enough to disrupt the inflaton oscillations, terminating the parametric resonance and ushering in the non-linear phase of dynamics that follows.
In this paper, we assess the viability of massless preheating in the context of electroweak vacuum metastability and address each of the above questions. But our analysis includes an important generalization: we allow non-minimal gravitational couplings for both the inflaton and the Higgs. In some sense, these couplings are unavoidable since if we ignore them in the tree-level action, they are generated at loop level [33]; even so, we are inclined to include them for several other reasons. First, a non-minimal Higgs-curvature coupling provides an additional stabilizing interaction for the Higgs during inflation without introducing additional non-SM field content. Second, a non-minimal inflaton-curvature coupling allows us to present a complete and self-contained model. The effect of the curvature coupling is to flatten the inflaton potential in the large-field region, restoring viability to the quartic inflaton potential in the inflationary regime, which is otherwise excluded by observations of the cosmic microwave background (CMB) [34]. Furthermore, during inflation, the non-minimal gravitational interaction is effectively scale invariant, so such couplings fit naturally within the purview of scale-invariant theories of inflation [35][36][37][38][39][40][41][42]. In this way, our study provides an analysis of the preheating dynamics that emerges from this well-motivated class of inflationary models.
Along with the generalization of non-minimal gravitational couplings, we shall consider a generalization of the gravity formulation. That is, in addition to the conventional metric formulation, we consider the so-called Palatini formulation, in which the connection and metric are independent degrees of freedom [43,44]. In a minimally coupled theory, such a distinction is not pertinent, but with non-minimal curvature couplings, there is a physical distinction between these formulations [45]; we shall consider certain facets of both in this paper.
Undoubtedly, the presence of non-minimal curvature couplings can have a substantial impact on the preheating dynamics. Because the curvature interactions break scale invariance during preheating, particle production from these terms dissipates over time and terminates, unassisted by backreaction effects. Consequently, if the initial curvature contribution is at least comparable to the quartic contribution, the system experiences a sequence of dynamical phases: particle production is dominated by the former immediately after inflation and transitions to the latter after some relatively short duration. Moreover, even if the curvature coupling is small, it can play a significant role. The curvature interaction contributes either constructively or destructively to the effective mass of the Higgs modes, depending on the sign of the coupling. As a result, the generated Higgs fluctuations can have an orders-of-magnitude sensitivity to this sign, which is ultimately reflected in the range of curvature couplings that are most vacuum-stable.
By broadly considering the dynamics of such massless preheating models, we place constraints on the space of Higgs couplings for which the metastability of the electroweak vacuum survives. And contrary to constraints that arise in other scenarios, we do not find a simply connected region. Instead, we find a large number of disjoint "islands of (meta)stability" over the parameter space, which merge into a contiguous metastable region at large quartic coupling. Accordingly, unlike other preheating scenarios-which typically lead to an upper bound on the Higgs-inflaton coupling-we find that massless preheating requires a more complex constraint, and a lower bound ultimately describes the most favorable region of metastability. In this way, the constraints necessary to stabilize the Higgs during inflation and stabilize it during preheating work in concert rather than in opposition. This paper is organized as follows. We begin Sec. II with an overview of the class of models we consider for our study of massless preheating and our assumptions therein. We show how the inclusion of non-minimal gravitational couplings allows for a self-contained model with a viable inflationary regime in the large-field region. The post-inflationary evolution of the inflaton and other background quantities is also discussed. In Sec. III, we examine the production of Higgs and inflaton particles in this background without yet considering the backreaction from these processes. At the close of the section, we analyze the effect of perturbative decays on the growth of Higgs fluctuations. In the penultimate Sec. IV, we investigate backreaction and the destabilization of the electroweak vacuum. We delineate the resulting metastability constraints on models in which massless preheating emerges. Finally, in Sec. V, we summarize our main results and possible directions for future work.
This paper also includes three appendixes. In Appendix A, we provide an overview of vacuum metastability in the context of "massive preheating," where the inflaton potential is approximately quadratic after inflation and the scale invariance of the inflaton potential broken. These scenarios have been studied in the literature [20,21,23,24], and we reproduce their findings for the purpose of comparing and contrasting to our results. Meanwhile, in Appendix B, we include the analytical calculations for the tachyonic resonance which are not given in the main body of the paper, and in Appendix C, we provide the relevant details of our numerical methods.

II. THE MODEL
Let us consider a model in which the Higgs doublet H is coupled to a real scalar inflaton φ, and both of these fields have a non-minimal coupling to the Ricci scalar curvature R. The model is described most succinctly in the Jordan frame by the Lagrangian 2 where we have used that H = 0 h/ √ 2 T in the unitary gauge. 3 The potential is restricted to the scale-invariant interactions motivated in Sec. I: and now interpreted in the Jordan frame. The first term in Eq. (2.2) determines the evolution of the cosmological 2 We employ a system of units in which c ≡ ≡ 8πG ≡ 1 and adopt the (+, +, +) sign convention of Misner et al. [46]. Additionally, we follow a sign convention for the non-minimal couplings ξ X (X = φ, h) such that the conformal value is ξ X = 1/6. 3 The initial state is symmetric under the SU (2) L × U (1) Y gauge transformation. However, once the fluctuation of the Higgsdoublet field is amplified by the instability, it easily becomes classical with a finite Higgs expectation value. Thus, the radial mode and the phase directions are well defined at each local position, and this justifies our use of the unitary gauge. Additionally, note that the initial quantum fluctuations cause the orientations of the Higgs field to be randomly distributed, resulting in nontrivial topological configurations (Chern-Simons number). This leads to gauge-field production [47], the effects of which are an interesting topic for future work.
background, while the second term is meant to stabilize the electroweak vacuum during inflation. The last term is the source of electroweak instability, with λ h running to negative values at energy scales exceeding ∼ 10 10 GeV. Let us make some ancillary remarks on the gravitational formulation used in what follows. In principle, one can formulate general relativity in two ways: (i) using the metric g µν and taking the connection Γ α βµ to be given by the Christoffel symbols (the "metric formulation") or (ii) using the metric and the connection as independent degrees of freedom (the "Palatini formulation"). Of course, these formulations are physically equivalent for a minimally coupled theory, but this distinction carries weight given the non-minimal couplings in our theory. Not favoring one formulation over the other, we address both of these possibilities by introducing a parameter and when relevant, we shall state our results as functions of θ. Although the distinctions between these formulations can be realized in both the preheating and inflationary dynamics, they figure most prominently in the large-field region of the potential, including the inflationary regime; let us now focus our discussion on this regime.

A. The Inflationary Regime
In principle, the model we have stated in Eq. (2.1) is agnostic to the specific form of the inflaton potential in the inflationary regime. As long as basic criteria are satisfied-e.g., the vacuum remains stable during inflation, isocurvature perturbations are negligible, etc.-our study of the preheating epoch is insensitive to details of the inflation model. Even so, we are motivated to examine whether Eq. (2.1) could consist of a phenomenologically viable inflation potential in the large-field region, as this would furnish a complete and self-contained model.
On the one hand, a minimal coupling ξ φ = 0 would be the simplest possible scenario. Unfortunately, this yields the standard quartic inflation, which is known to predict too large a tensor-to-scalar ratio to satisfy observations of the CMB [34]. On the other hand, a finite ξ φ < 0 allows for a range of other possibilities [35][36][37][38][39]. In considering these, it is beneficial to transform to the Einstein frame, defined by the Weyl transformation g µν −→ Ωg µν in which 4 (2.4) While this restores the canonical graviton normalization, it also transforms the potential to and generates non-canonical kinetic terms for our scalar fields X = {φ, h} such that the Lagrangian becomes with a kinetic mixing matrix (2.6) In the large-field region φ 1/ −ξ φ , the inflaton potential then takes the form where φ is the canonically normalized inflaton field. In this region, the theory acquires an approximate scale invariance and effectively approaches so-called "attractor models" which are phenomenologically viable [49][50][51]. In particular, for an inflationary epoch of N e-folds, one finds a spectral index n s = 1 − 2/N and tensor-to-scalar ratio r = 12α/N 2 , where we have defined A plot which illustrates the inflaton potential V E ( φ) for a number of possible curvature couplings ξ φ and gravity formulations θ, all of which lead to massless preheating after inflation, is shown in Fig. 1. The inflationary trajectory is shown by the thick part of each curve, and inflation ends at the location of the point marker.
There are several constraints on our model parameters necessary to achieve viability. First, matching the observed amplitude of primordial curvature perturbations [34] fixes the relationship between λ φ and ξ φ : Second, we must ensure that the electroweak vacuum remains stable throughout inflation. This requirement bounds the effective mass squared of the Higgs as under the slow-roll approximation, where h is the canonically normalized Higgs field. The slow-roll-suppressed contributions to the Higgs mass neglected above can become more important toward the end of inflation, so the precise condition may be modified. In some inflation models, this causes inflation to end prematurely [52]; this may happen for ξ φ −1 in the Palatini formalism, which is outside the scope of this paper. The constraint in Eq. (2.9) is simplified if we consider that during inflation 1 −ξ φ φ 2 and thus Eq. (2.9) reduces to λ φ ξ h − g 2 ξ φ > 0. We shall assume that this constraint is satisfied to the extent that the Higgs is stabilized strongly at the origin during inflation. In other words, we assume that the Higgs mass is larger than the Hubble scale H during inflation, and imposing this assumption, we obtain (2.10) Third, we must also ensure that quantum corrections to the scalar potential from Higgs loops are controlled and do not ruin the above predictions for inflationary observables. To this end, if we require that g 2 λ φ , then the Higgs-loop contribution is subdominant. Taking a value λ φ = 10 −10 [based on Eq. (2.8)] gives the constraint explicitly as g 2 /λ φ 10 5 , which is well within the parameter space that we shall consider.
Finally, the accelerated expansion of the universe ends once the energy density of the inflaton background falls below ρ φ = 3V (φ)/2, i.e., when the field falls below . (2.11) Afterward, φ begins to oscillate about the minimum of the potential and we identify this point in time as the beginning of the preheating era; we focus our discussion throughout the rest of the paper on this epoch.

B. Cosmological Background After Inflation
Let us now discuss the evolution of the cosmological background after the end of inflation. Assuming the vacuum has been sufficiently stabilized, the Higgs field is negligible, 5 and the cosmological background is determined solely by the dynamics of the inflaton and its potential V E (φ) V (φ) in the region φ < φ end , where the notation V (φ) is introduced for the inflaton potential approximated in the small-field region.
The form of the inflaton potential in this region may still be sensitive to ξ φ . Notably, for |ξ φ | 1 an interesting distinction appears between metric and Palatini gravity. For metric gravity, the potential is quartic at sufficiently small φ 1/|ξ φ | but becomes quadratic in the intermediate region 1/|ξ φ | φ 1/ |ξ φ |. It follows that Higgs fluctuations are amplified by two different types of parametric resonance depending on the field region. 6 By contrast, in Palatini gravity the potential is purely quartic and has no quadratic region [55,56]. Our model therefore offers an explicit example of purely massless preheating that is consistent with inflation, even for a large non-minimal coupling ξ φ . Furthermore, we can achieve the same effect with ξ φ = O(1), as the intermediate quadratic region vanishes. Given that the study of massless preheating is our main interest and that this occurs in the small-field region, we shall henceforth assume that |ξ φ | 1 (with ξ φ < 0).
After the end of inflation, the evolution of the background inflaton field is then well approximated bÿ in which the dots correspond to time derivatives and H ≡ȧ/a is the Hubble parameter given in terms of the scale factor a = a(t). As discussed in Sec. I, the approximate scale invariance of the system makes its field dynamics and resonance structure rather unique. These features have been studied extensively [19] and here we briefly review them. The scale invariance is made transparent by writing Eq. (2.12) in terms of the conformal time η ≡ t dt /a(t ) and conformal inflaton field ϕ ≡ aφ: where the prime notation corresponds to η derivatives. Note that we have ignored a term proportional to φ 2 R, as this term negligibly impacts the inflaton evolution. The approximate scale invariance of the theory is manifested by the fact that the inflaton equation of motion in Eq. (2.13) is independent of the cosmological expansion. The solutions carry a constant amplitude ϕ of oscillations 5 The two-field evolution that one finds in breaking from this assumption is non-trivial and has been investigated in Ref. [53]. 6 In this context, we refer the reader to Ref. [54], in which the electroweak instability was studied (with g = 0 and ξ h = 0) in the large non-minimal inflaton coupling regime −ξ φ 1.
and are given in terms of a Jacobi elliptic function 7 in which x ≡ λ φ ϕη is the conformal time measured in units of the effective inflaton mass λ φ ϕ. The constant x 0 is used to match to the inflaton configuration at the end of inflation, and the period of oscillations in x is with K(X) ≡ π/2 0 dθ(1 − X 2 sin 2 θ) −1/2 defined as the complete elliptic integral of the first kind.
A cosmological energy density which is dominated by the coherent oscillations of a scalar field may behave in a variety of ways depending on the scalar potential. For example, the class of potentials V (φ) ∝ φ 2n (for integer n > 0) yield a cosmological background that behaves as a fluid with the equation-of-state parameter [57] w ≡ P ρ = n − 1 n + 1 , (2.17) in which P is the pressure and ρ is the energy density, averaged over several oscillations. For the quadratic (n = 1) and quartic (n = 2) potentials this demonstrates the well-known result that scalar-field oscillations in these potentials correspond to perfect fluids with matterlike (w = 0) and radiation-like (w = 1/3) equations of state, respectively. This behavior reflects our observation in Eq. (2.13) that ϕ evolves independently of the cosmological expansion. Namely, since the inflaton energy density scales like radiation ρ φ ∝ 1/a 4 , the amplitude of inflaton oscillations scale as φ ∝ 1/a. Therefore, the corresponding conformal amplitude ϕ is fixed. It also follows that in a radiation-like background the scale factor is proportional to x and grows according to Inflation ends once the kinetic energy grows sufficiently to have 3V (φ)/2 ≤ ρ φ , which corresponds to the time x end ≡ √ 12/φ end . Note that chronologically one always has x 0 < x end and these are related explicitly by for simplicity we have used ξ φ = 0 in this expression. 7 We define the Jacobi elliptic sine sn(X, Y ) = sin Z and cosine cn(X, Y ) = cos Z functions through the relation (2.14) Finally, in addition to the inflaton background, it is important that we examine the scalar curvature after inflation. In general, the curvature is a frame-dependent quantity with ΩR J = 4V J (φ) −φ 2 in the Jordan frame. Nevertheless, at sufficiently small φ we have Ω ≈ 1 and where we have employed the solution in Eq. (2.15). The scalar curvature thus oscillates about zero and over several oscillations averages to R = 3H(1 − 3w) = 0. However, as we shall find upon examining particle production, neither the curvature terms nor their time dependence can be neglected, as they can impart a significant contribution to the preheating dynamics.

III. PRODUCTION OF HIGGS PARTICLES
Having established the evolution of the classical inflaton field in Sec. II B, we can now discuss Higgs particle production in this background. We write the quantized Higgs field h in the Heisenberg picture as a function of the fluctuations h k (t) of comoving momenta k: where a † k and a k are creation and annihilation operators, respectively. For a given comoving momentum, these fluctuations follow equations of motion where Γ h k is a phenomenological term accounting for perturbative decays of the Higgs [18,29,31] and ω h k is the energy of the mode. These modes are coupled to both the oscillating inflaton background and the scalar curvature R, and therefore their effective masses carry an implicit time dependence. In the Jordan frame, However, in the Einstein frame, ξ φ = 0 generates a kinetic mixing between the inflaton field and the Higgs modes, thereby producing an inflaton-dependent friction term. We can absorb this friction into the effective mass term by the field redefinition H k ≡ aΩ −1/2 | h=0 h k [54], yielding the equations of motion in which the transformed modes are given by and we have defined the effective non-minimal coupling Note that we have neglected terms which are higher order in a −1 , such that the scalar curvature is given by In this way, we have absorbed most ξ φ effects and dependence on the gravity formulation into the single effective parameter ξ. Finally, we confirm that if we take ξ φ = 0, the scale invariance is restored for the conformal value ξ h = 1/6, as one would expect. Solving the equations of motion in Eq. (3.4), we can track the production of Higgs particles. In particular, the comoving phase-space density of particles associated with a mode of comoving momentum k is given by The physical mechanism that drives this production differs considerably between the ξa 2 R and g 2 ϕ 2 terms. For the former, when the inflaton field passes through the minimum of its potential, one may find that a range of Higgs modes become tachyonic ω 2 H k < 0. The tachyonic instability is strongest for the smaller-momentum modes, and these modes produce particles for longer durations. By contrast, oscillations in ϕ drive particle production from the latter. The resulting time-dependent modulations of ω H k give rise to parametric resonances for Higgs modes within certain momentum bands.
Another crucial distinction between these mechanisms is found by evaluating their overall scaling under the cosmological expansion. The main term in Eq. (3.5) responsible for tachyonic production redshifts as ξa 2 R ∝ 1/a 2 , while the term responsible for the parametric resonance g 2 ϕ 2 does not dissipate at all. Tachyonic production therefore always terminates after some finite time, even for the zero-momentum mode. On the other hand, production from the parametric instability continues unimpeded and ceases only once the evolution is disrupted by backreaction effects, as we cover in Sec. IV. Indeed, this distinction is ultimately traced to the Higgs-curvature interaction breaking the scale invariance that is preserved by all of the other relevant terms.
We first break our analysis into two limiting cases: one in which the quartic inflaton-Higgs interaction is dominant and the other in which the curvature interaction is dominant. Then, we explore the interplay between these interactions and finally begin to analyze the impact of perturbative Higgs decays on the field dynamics.

A. Production from Parametric Instability
Let us first consider the case that the curvature coupling ξ is negligible and thus ignore the tachyonic production. Then, Higgs particles are produced purely from the parametric resonance associated with the g 2 ϕ 2 term in Eq. (3.5). The dominant production in this case arises from the fact that as the inflaton passes through the origin, the effective masses may evolve non-adiabatically: which triggers a burst of particle production. 8 The growth of the number density for a given mode is exponential and follows log n k 2µ k x over several oscillations, where µ k is the characteristic exponent. Note the distinction between this regular exponential growth and the stochastic growth one finds for theories without an approximate scale invariance, e.g., those with a quadratic inflaton potential [18,19]. The stochastic nature of the resonance appears in these scenarios because the accumulated phase of each mode evolves with the cosmic expansion, destroying the phase coherence. In the scaleinvariant theory, no such time dependence may arise and phase coherence is maintained among the modes. A more in-depth comparison of the quadratic and quartic theories, with an emphasis on the results of this paper, is provided in Appendix A.
The size of a particular growth exponent µ k is determined by a combination of the comoving momentum k and the quotient g 2 /λ φ . In Fig. 2, we have numerically solved the mode equations and plotted contours of µ k in the space of these two quantities, rescaling the momentum as κ ≡ k/( λ φ ϕ). It is natural to separate the resonances into two different classes. The couplings with instability bands that include the zero-momentum mode, i.e., those between 2n 2 − n < g 2 /λ φ < 2n 2 + n, for n ∈ N, contain the broadest resonances and generally give the most copious particle production-we refer to these collectively as the "broad regime." Meanwhile, those couplings with only finite-momentum bands contain the most narrow resonances and give generally weaker production, and we refer to these as the "narrow regime." The most weak and narrow bands occur at the boundaries g 2 /λ φ = 2n 2 + n. The distinctions between these two coupling regimes have important dynamical implications and play a major role in this paper.
Although, in principle, each mode with µ k > 0 contributes to particle production, in practice the maximum exponent of the band µ max ≡ max k µ k dominates. In the broad regime, the maximum exponent is µ max ≈ 0.24, found at the central values g 2 /λ φ = 2n 2 , and this value is universal over the entire range of couplings. Conversely, the µ max values in the narrow regime are not universal. These features are most easily observed in the top panel of Fig. 2, where µ max is seen to oscillate back and forth between the broad and narrow regimes. Indeed, the growth rate of particle production is highly non-monotonic as one dials the inflaton-Higgs coupling. 8 A number of the results we present in this subsection are the subject of Ref. [19]; we summarize only the most relevant aspects.
The instability bands of the parametric resonance arising from the inflaton-Higgs interaction (bottom panel) and the corresponding maximum characteristic exponents µmax ≡ max k µ k for each coupling (top panel). The contours in the bottom panel show the value of the exponent µ k such that a given occupation number grows as n h k ∝ e 2µ k x . The couplings which lead to the strongest growth are found at the center of bands which have an unstable zero mode, i.e., for g 2 /λ φ = 2, 8, . . . , 2n 2 for n ∈ N, while the weakest are found at the edge of these bands g 2 /λ φ = 3, 10, . . . , 2n 2 + n-we refer to these as the "broad" and "narrow" regimes, respectively. There is a universal µmax ≈ 0.24 for the former, while for the latter µmax is a non-trivial function of g 2 /λ φ , given in Eq. (3.9) and plotted over an extended range in Fig. 3.
In fact, the narrow resonances grow to join the broad resonance at µ max ≈ 0.24 in the formal limit g 2 /λ φ → ∞.
The non-monotonicity of µ max thus diminishes as a function of g 2 /λ φ , albeit at an extremely slow pace. In what follows it proves useful to quantify this observation, so we have numerically computed µ max for the discrete minima of µ max in the narrow regime g 2 /λ φ = 2n 2 + n and found that they are well approximated by the function The constants A = 2.82 and B = 3.34 × 10 5 reproduce the numerical results to better than 1% accuracy over the range 3 g 2 /λ φ 10 4 , which spans the full range of narrow resonances relevant to this paper. We display the function in Eq. (3.9) and the numerical results together in Fig. 3 for comparison. Note that the weakest resonances globally are found in the limit g 2 /λ φ → 0, where the bands become increasingly narrow and follow µ max ≈ 0.15g 2 /λ φ [19]; we shall discuss this small-coupling limit further in Sec. III C 3.
We have plotted numerical solutions of the phase-space density n h k in the broad and narrow regimes, for adjacent coupling bands, shown by the black curve in the top and bottom panels of Fig. 4, respectively. The mode corresponding to the most rapid growth is shown in both cases: for the broad regime this is the κ = 0 mode, but for the narrow regime µ max corresponds to a finite momentum. We neglect the blue curves for the moment, as these first enter our discussion in Sec. III D.
For our purposes, the non-monotonic nature of µ max as a function of g 2 /λ φ has extensive implications. In contrast to many other preheating scenarios, the magnitude of our coupling has no bearing on the growth rate of a given fluctuation-only the associated band within the repeating resonance structure is important. That said, since the width of the momentum bands increases with g 2 /λ φ , the total number density n h does, in fact, depend on this coupling. We obtain the total comoving number density of the produced Higgs particles by using the saddlepoint approximation to integrate over each band: As this density partly determines the variance of the Higgs fluctuations, it plays a major role in assessing the metastability of the electroweak vacuum. Accordingly, we shall continue to calculate n h in all of the regimes.

B. Production from Tachyonic Instability
Let us now consider the opposite coupling limit g → 0, i.e., the limit in which particle production is driven not 3 1000 2000 3000 4000 The maximum exponents µmax ≡ max k µ k for the narrow-resonance couplings g 2 /λ φ = 2n 2 + n (for n ∈ N). These are computed numerically (point markers) and compared to the analytical fitting formula in Eq. (3.9) (solid curve). Note that the analytical function is evaluated only at the same discrete values g 2 /λ φ = 2n 2 + n in this figure, i.e., the minima of the top panel of Fig. 2. In the g 2 /λ φ → ∞ limit the broad/narrow regimes become degenerate with µmax ≈ 0.24, but this asymptote is approached slowly.
x end 10 20 30 40 50 60 70 x ≡ λ φ ϕη FIG. 4. The evolution of background quantities ϕ, a 2 R (top panel) and Higgs phase-space density n h k (center/bottom panels) during the early stage of preheating. The latter two panels show the broad/narrow regimes: taking n = 12, the center panel uses a coupling g 2 /λ φ = 2n 2 and has a broad range of resonant modes, while the bottom panel uses the adjacent g 2 /λ φ = 2n 2 + n and has only a narrow range of resonant modes; the mode with the largest growth exponent is shown in each case. Additionally, these two panels show the effect of allowing for perturbative Higgs decays (first discussed in Sec. III D, shown in blue) and non-minimal gravitational couplings ξ h , evenly spaced over the range 0 ≤ ξ h ≤ 50.
by parametric resonance but by a tachyonic instability. The effective masses ω H k are given by where we have defined a quantity r h ≡ ξϕ 2 /(2a 2 ) that indicates the strength of the curvature term at a given time and utilized Eq. (2.20). We recall from Sec. II A that vacuum stability during inflation requires λ φ ξ h > ξ φ g 2 , and, therefore, we consider only ξ h > 0 for the moment. Indeed, for sufficiently small momenta, one finds modes that cross the tachyonic threshold ω 2 H k < 0 when the inflaton is near the minimum of its potential, triggering an exponential growth in the corresponding Higgs fluctuations. Although this growth is typically short-lived due to the redshifting of the curvature term a 2 ξR ∝ 1/a 2 , these may still be a source of copious particle production soon after inflation and thus serve as a legitimate threat to destabilize the electroweak vacuum.
There have been a number of studies devoted to calculating the rate of tachyonic particle production in different settings [15,16], and we employ several of those techniques in this section. Similar to Sec. III A, near the turning points ω 2 H k = 0 the masses may change nonadiabatically and one must then compute the Bogoliubov coefficients to proceed. However, unlike in the previous section, the ω 2 H k may experience two distinct adiabatic segments of evolution. As long as r h 1, the effective masses change adiabatically away from the turning points in both the tachyonic and non-tachyonic segments. Under this assumption, we can apply the WKB approximation by calculating the phase accumulated by modes both during the tachyonic Applying these methods, one finds that after passing through a tachyonic region j times, the phase-space density for a given mode is written generally as [16] Hence, employing Eq. (3.11) we can compute the accumulated quantities X k and Θ k for g = 0. The details of this calculation appear in Appendix B and give (3.13) After successive bursts of tachyonic production, the growth exponent in Eq. (3.12) accumulates a value j 2jX k 4 dxX k /T , with the zero mode receiving the greatest share of the number density. The accumulated phases Θ k supply only oscillatory behavior or modify the distribution over momenta. Given that our primary concern is the overall growth of the Higgs number density, we shall neglect these quantities. Then, we can estimate that which holds for the duration of the tachyonic instability. Let us focus on the κ = 0 mode, which experiences the strongest growth in the tachyonic regime. At first glance, the growth may actually appear weak in comparison to the parametric resonance [in Eq. (3.10)] since it is not exponential in time-it merely obeys a power law. The difference is that the power scales with √ ξ and has no upper bound, in line with studies of the tachyonic instability in other settings [16,21]. This feature sharply contrasts with the parametric resonance, in which the µ max growth rate is bounded universally from above (as we observed in Fig. 2), regardless of the coupling. Of course, as the universe expands the tachyonic production soon terminates. In particular, as r h redshifts to values below unity, the tachyonic masses Ω h k are suppressed and the adiabatic assumption breaks down. Using |Ω h k | Ω 2 h k as the threshold for where this breakdown occurs, we find r 2 h (r h − κ 2 ) 3 . This condition implies that the span of modes exposed to the instability is bounded above by κ √ r h and that a given mode is active for x √ 6ξ/max(1, κ). The modes shut down successively, starting with the largest-momentum modes, such that tachyonic production ends at the time The total number density of Higgs particles is given by integrating Eq. (3.14) over the phase space, but a cutoff κ r 2 h should be imposed on each momentum band per the discussion above. We again use the saddlepoint approximation to perform the integration and obtain a comoving number density which is applicable for x x ξ . Although there is a brief transient phase of non-adiabatic production for x x ξ , particle production from the curvature term proceeds only through the substantially weaker narrow resonance, which shuts down entirely soon thereafter. We cover the details of this regime in Sec. III C 3 below.

C. Production in the Mixed Case
Let us now promote our discussion to the mixed case, in which both couplings g 2 /λ φ and ξ are non-zero. As such, the Higgs modes evolve with the effective masses (3.17) that follow directly from Eq. (3.5), and by analogy to There are several immediate implications; let us discuss these with a focus on the effect of the Higgscurvature coupling. First, since only the curvature term dissipates with the cosmological expansion, the dynamics may progress through several distinct phases, most noticeably if r h g 2 /λ φ at early times. Second, the reintroduction of the g 2 ϕ 2 term lifts the effective masses and thereby opens the ξ < 0 region to viability. Indeed, the ξ < 0 region was excluded in Sec. III B only because the inflationary constraints [in Eq. (2.9) and Eq. (2.10)] would have been violated, but for mixed couplings this half of parameter space can be reincorporated. Third, the small-coupling regime (g 2 /λ φ 1 and r h 1) stands apart from most of our discussion thus far. The particle production in this regime is driven entirely by narrow parametric resonances arising from two different types of modulations (ϕ/ϕ) 2 and (ϕ/ϕ) 4 . If these modulations can be approximated as sinusoidal, then a perturbative treatment is likely effective. We investigate all of these possible scenarios below.
The simplest possibility is an initially small curvature term |r h | g 2 /λ φ , since this does not allow a tachyonic instability to develop during preheating. (We confine our discussion to g 2 /λ φ 1 for the moment and cover the small-coupling regime separately.) Nonetheless, the curvature term in this case has an impact on the Higgs dynamics. Depending on the sign of ξ, these terms may add constructively or destructively, enhancing or suppressing the effective masses ω H k , respectively. For ξ > 0, the effect is constructive and the resonance bands are widened; e.g., the broad instability bands are extended to κ 2 ≤ 2g 2 /(π 2 λ φ ) + r h . We have seen that Fig. 4 falls within this category, and this extension explains the variation among different ξ. Conversely, for ξ < 0 the effect is destructive: not only do the bands narrow, but the non-adiabaticity that drives the parametric resonance is weakened. Explicitly, the maximum characteristic exponent is reduced to Surely, after a sufficient duration one has r h g 2 /λ φ and the exponent µ max ≈ 0.24 is restored. In other words, the full capacity of the parametric resonance is delayed until a time x g , which is given approximately by After this time the curvature term continues to dissipate and has an increasingly negligible influence. The growth of fluctuations then proceeds according to Sec. III A.
2. Dominant ξRh 2 (with g 2 /λ φ 1) An initially large r h g 2 /λ φ 1 ensures the Higgs dynamics are dominated by tachyonic production until a time x ξ . Unlike the converse situation above, where we could ignore the curvature interaction after some duration, there is no regime in which we can ignore the g 2 φ 2 interaction entirely. Aside from the fact that the 0 1 2 3 4 parametric resonance from this term will always play a role, its presence also alters the size of the tachyonic masses. We can incorporate this effect by performing an analysis along the lines of Sec. III B with g 2 /λ φ = 0. Then, for a given mode, we find that when g 2 /λ φ r h g 2 /λ φ + κ 2 is satisfied the tachyonic instability is active and proceeds adiabatically. In order to compute the integral analytically, we assume that Then, the accumulated quantities [from Eq. (3.12)] are found to be (3.20) As expected, Eq. (3.20) shows that tachyonic production is still concentrated in the small-momentum modes. However, the evolution of n h k in the presence of g 2 /λ φ = 0 is markedly different. Rather than a power law, the phase-space density rapidly asymptotes to a constant as again neglecting the oscillatory component cos Θ k . The termination of the tachyonic resonance is pushed to an earlier time depending on the quartic couplings: Afterward, the modes grow as n h k ∝ e 2µ k x , driven by the parametric resonance. The phase-space density can be found at these times by matching to Eq. (3.21) at x ξ . Finally, the number density of Higgs particles produced during this tachyonic phase is found by integrating Eq. (3.21). Using the saddlepoint approximation we find In Fig. 5 we have plotted the phase-space density n h k , calculated numerically over the momentum space of the Higgs modes. The quartic coupling chosen g 2 /λ φ = 21 sits on a narrow resonance and therefore produces particles only at finite momenta (highlighted by the gray region). Meanwhile, the different curves show various non-minimal couplings ξ. Because the parametric resonance yields unimpeded particle production, the finitemomentum peak in n h k continues to grow with time, while production for the other modes ceases once x x ξ . The large point markers along the vertical axis show agreement with our analytical result for n h0 where the tachyonic production is dominant, found by evaluating Eq. (3.21) at zero momentum.
The same logic as above may be applied to the ξ < 0 region of parameter space. As this region is dominated by tachyonic particle production, we may again employ Eq. (3.12) to calculate the generated Higgs spectrum. While this calculation is straightforward in principle, it requires some additional technical details that we include in Appendix B, but we summarize the results here. The phase-space density for the Higgs is given by in this limit. Compared to ξ > 0, the power indices are a similar magnitude if we neglect the effect of ξ φ . In Fig. 6, we demonstrate the growth of the zero-mode n h0 in the mixed-coupling case, showing both our numerical (thin curves) and analytical (thick curves) results. For comparison, we have included both the positive (blue curves) and negative (red curves) ξ regimes, 20 40 60 80 x ≡ λ φ ϕη The thin (thick) curves show our numerical (analytical) results, respectively. For ξ > 0, the tachyonic instability drives particle production early on [following Eq. (3.21)], but once x x ξ it is driven by the parametric resonance n k ∝ e 2µ k x . The quartic and curvature contributions add constructively and enhance the instability. Meanwhile, for ξ < 0, these contributions add destructively, and the growth exponent µ k is effectively shut off [following Eq. as well as the narrow and broad resonance regimes. The destructive effect of ξ < 0 is evident, as changing the sign of ξ induces a gap of many orders of magnitude in the zero-mode density. As a result, we should expect that electroweak vacuum metastability shows a preference for the ξ < 0 half of parameter space. We shall return to this topic when we examine Higgs destabilization in Sec. IV.

Small-Coupling Regime
Let us now focus on the small-coupling regime, i.e., where both g 2 /λ φ 1 and |r h | 1. The latter inequality is always imminent, so we are free to also interpret this regime as the late-time behavior of a model with g 2 /λ φ 1 and an arbitrary curvature coupling. Indeed, as soon as |r h | 1 the tachyonic instability transitions into a narrow resonance. In fact, both of the instabilities in this regime experience a narrow parametric resonance, and our analysis should follow a different approach. Namely, the oscillatory terms in Eq. (3.17) act as small modulations to ω 2 H k , and we can treat these terms perturbatively [19]. We expand the elliptic functions as  with the leading five coefficients given in Table I. Noting that the series converges quickly after the first three terms, we truncate at order = 2 and the expansion serves as a good approximation. Moreover, given that r h are slowly varying compared to the timescale of inflaton oscillations, the Higgs mode equations take the form of the Whittaker-Hill equation: where we have defined z ≡ 2π(x − x 0 )/T . The coefficients of the frequency terms are given by in which we have defined the quantities Using Floquet theory it is relatively straightforward to calculate the characteristic exponents µ k for the unstable modes [58,59]. We can therefore find the maximum exponents µ max ≡ max k µ k over the small-coupling parameter space. These results are displayed in Fig. 7.
An interesting feature of Fig. 7 is that, given the scalefactor dependence of r h , we can interpret the figure as showing the flow of µ max as a function of time. For some initial r h = 0 value, µ max flows along contours of constant g 2 /λ φ (e.g., the gray arrows) until it reaches the r h = 0 axis. Notably, for ξ < 0 the evolution of µ max over this time is not necessarily monotonic. For larger quartic couplings, one may enter the µ max = 0 region (enclosed by red curves) and then exit it while approaching r h = 0.
Note that in the regime of negligible curvature couplings, Eq. (3.27) reduces to the form of the Mathieu equation. The Mathieu equation is familiar from models of massive preheating, but unlike such models the parametric instability from the g 2 φ 2 interaction does not terminate due to the approximate scale invariance of the theory. The fact that the curvature contributions are neither long-lived nor tachyonic means that the quartic contribution becomes the most salient in the small-coupling regime. Even for couplings as large as r h = O(g 2 /λ φ ), the resonance bands are widened, but their influence only induces a logarithmic correction to the growth rates. For these reasons, the r h → 0 limit is usually the most pertinent in the context of vacuum metastability. Here, the modes grow as n h k ≈ e 2µ k x /2 for momenta around the narrow bands. The primary resonance is centered at κ 2 = 2π/T , and one finds the maximum exponent µ max = F 1 T g 2 /(8πλ φ ) ≈ 0.15g 2 /λ φ [19]. Using the saddlepoint approximation, we integrate to find (3.30)

D. Perturbative Higgs Decays
Thus far, for simplicity we have neglected the perturbative decay of Higgs particles after their production. These decays may not appear important when considering the substantial rates of steady non-perturbative particle production. After all, perturbative decays are known to have only a minor effect in preheating with massive inflaton potentials (see Appendix A). However, given the unique properties that appear for massless preheating, it is worth examining the Higgs decays in closer detail.
Before continuing on this path, let us briefly address the other ways in which the Higgs may transfer its energy density, most notably the non-perturbative production of gauge bosons. Even though the Higgs mass term oscillates coherently due to the inflaton-field motion, the gauge-field mass terms do not oscillate uniformly. There-fore, in the early stages of preheating, we expect that resonant production of W and Z gauge bosons from the Higgs evolution is less significant than resonant production of the Higgs field itself. Still, the averaged value of the gauge-boson mass terms are modulated by the time dependence of h 2 , so a substantial amount of gauge bosons could be produced. As long as the Higgs field is much smaller than the critical value in Eq. (3.31), a simple analysis shows that the gauge-field production from fast oscillations of h 2 is in the narrow regime. More precise estimation of gauge-field production requires a dedicated study (see also footnote 3). In the following, we shall therefore assume that non-perturbative gaugefield production induced by the Higgs field can be ignored and that perturbative decays are the dominant mechanism through which energy is transferred out of the Higgs field.
In the early stages of preheating, we expect that perturbative decay channels for the Higgs are kinematically accessible since the homogeneous background value for h is negligible. Even if a background value is generated, the masses of SM particles are lighter than the Higgs during preheating as long as h h crit , in which is the Higgs value at the barrier separating the false and true vacua, i.e., the local maximum of the potential. Therefore, the dominant decay channel of the Higgs is into top quarks 9 at a rate Γ h 3y 2 t m h 16π (3.32) in the rest frame, where y t is the top Yukawa coupling, evaluated at the Higgs mass scale, and we have denoted as the effective Higgs mass. In general, we should include a Lorentz factor γ = ω h k /m h which suppresses the decay rate for relativistic particles. The tachyonic instability produces particles with physical momenta k/a φ 2 ξλ φ /2, which translates to γ O(1) and does not appreciably suppress the rate. Likewise, the parametric instability with g 2 /λ φ 1 produces particles with k/a [λ φ φ 2 g 2 /λ φ ] 1/2 , which is non-relativistic. The exception is small couplings g 2 /λ φ 1, for which Higgs production is dominated by particles with a Lorentz factor γ k/(am h ) λ φ /g 2 . While we include the full effect of this time-dilation on the decay rate in our numerical computations, our analytical calculations are performed assuming that γ max(1, λ φ /g 2 ).
Naturally, as perturbative decays exponentially suppress the number density, they work against the nonperturbative production processes above. The field dependence of the effective Higgs mass m h implies an interesting chronology of events. If the parametric resonance is the dominant production mechanism, then a burst of particles is produced as φ passes through the origin. Meanwhile, decays are strongest when φ is maximally displaced. The result is that rather than n h k maintaining a constant value during its adiabatic evolution, it is exponentially suppressed at the rate Γ h . However, due to its dependence on evolving background quantities, the decay rate is a non-trivial function of time. The effect on the number density is to dissipate it as where the integration is performed over times when decays are kinematically possible. An important limit is that in which the non-minimal couplings ξ h , ξ φ are negligible, as this also represents the system for times x x ξ . The pivotal observation is that the effective decay rate scales as aΓ h ∝ |ϕ| and thus has no overall scale-factor dependence, allowing a direct competition between the production and decay rates to determine the fate of the electroweak vacuum; this contrasts with massive preheating, in which the decay exponent carries a logarithmic time dependence. (We refer to Appendix A for a more thorough comparison.) Explicitly, the time-averaged conformal decay rate aΓ h is given by where we have assumed the regime g 2 /λ φ 1. Indeed, we find that the number density of Higgs particles does not grow for couplings exceeding Remarkably, this result suggests that electroweak vacuum metastability can be achieved in massless preheating. Moreover, metastability is given not by an upper bound but by a lower bound on the quartic coupling (see also Refs. [60,61]).
In the opposite limit, where g is negligible and the curvature term dominates, neither the production or perturbative decays evolve exponentially. In particular, we have aΓ h ∝ log a and a growth rate that is dominant at early times. The result is that perturbative decays do not quell the early phase of tachyonic production, but if the vacuum survives, decays increasingly dissipate the fluctuations as the system evolves.
Looking back to Fig. 4, we have demonstrated the effect of the perturbative decays on particle production by plotting n h k both with (blue curves) and without (grayscale curves) the decays incorporated. Clearly, in the absence of decays, the phase-space density remains constant when the inflaton is away from φ ≈ 0, as the n h k correspond to an adiabatic invariant. However, when the Higgs decays are turned on, the curves show a dissipation during the adiabatic evolution. The dissipation does not always appear to be at a constant rate, reflecting that in our numerical computations, the decay rates are not averaged as in Eq. (3.35) but are given instantaneously as a function of the effective Higgs mass m h .
Before concluding this section, we put aside the Higgs for a moment and make some remarks regarding the quantum fluctuations of the inflaton field. An extensive part of the discussion in this section can be applied by analogy to the production of inflaton particles. These particles are generated through the inflaton self-coupling, so that the effective masses [analogous to Eq. (3.3)] are in the Jordan frame, or also in the Einstein frame after neglecting higher-order terms in O(a −2 ). In fact, production via this inflaton "self-resonance" is understood by the simple replacement g 2 /λ φ = 3 in the analysis of Sec. III A [19]. For example, from our result for the Higgs number density n h in Eq. (3.10), we can estimate the inflaton number density n φ with appropriate replacements. Interestingly, the characteristic growth exponent µ max ≈ 0.036 for inflaton production corresponds to the weakest of the g 2 /λ φ ≥ 1 regime. In Sec. IV below, we shall find that this self-resonance and its properties are extremely material to the discussion of vacuum destabilization and the end of preheating.

IV. BACKREACTION AND VACUUM DESTABILIZATION
While Sec. III provides the groundwork for our study of electroweak vacuum metastability, an essential ingredient is still absent from our analysis. In particular, recalling our assumption that the Higgs field has a negligible background value h ≈ 0, the term in the potential which threatens to destabilize the vacuum (the self-interaction λ h h 4 /4) has not yet entered our analysis. In order to incorporate the unstable vacuum, we must consider not only the growth of Higgs fluctuations but also the backreaction that these fluctuations have on the system. As discussed in Sec. I, although the fluctuations in massless preheating appear to grow unimpeded, they will inevitably grow large enough to either trigger vacuum decay or disrupt the background inflaton evolution through backreaction.
A complete treatment of the backreaction and nonlinearities that appear, especially in the later stages of preheating, requires numerical lattice simulations. However, we find that, in order to probe the Higgs stability, or estimate the onset of the non-linear stage, the oneloop Hartree approximation [62] is sufficient. That is, we assume the factorization h 4 → 6 h 2 h 2 − 3 h 2 2 for the quartic term in the Lagrangian (and an analogous expression for the inflaton), where the expectation values are given by In turn, the effective masses for the Higgs and inflaton modes-from Eq. (3.3) and Eq. (3.37), respectively-are modified as the variance of the fluctuations grows: The fluctuations also couple to the inflaton background, modifying the equation of motion as Given that the motion of φ drives the particle-production processes, once these oscillations are disrupted the dynamics of the type discussed in Sec. III is shut down. The termination of these processes generally coincides with the energy density of the fluctuations growing comparable to that of the inflaton background and thus corresponds to the end of the linear stage of preheating.
We shall include all of the Hartree terms in our numerical calculations (as detailed in Appendix C), but the corrections from the self-interactions 3λ h h 2 and 3λ φ φ 2 will tend to carry the most weight for our analytical calculations. We find that from simple analytic estimates we can reproduce the salient lattice simulation results in the literature. For further validation, our numerical analyses are also compared to the literature on models for which the inflaton potential is quadratic in Appendix A. For further discussion on the non-linear dynamics and backreaction in preheating, we refer the reader to Refs. [18,19].
The main results of Sec. III consist of the produced number density n h , so it is necessary to translate between these quantities and the variance h 2 . While one can fully express h 2 in terms of Bogoliubov coefficients, the result is a sum of two terms, one of which is rapidly oscillating and not important for our purpose of producing order-of-magnitude estimates [18]. As long as the produced Higgs particles are non-relativistic one can write h 2 n h /(a 3 ω h ), with a similar expression applying to the inflaton or any other relevant fields. Note that the number density is well defined only when evolving adiabatically, so the points at which we evaluate ω h k must remain consistent with this assumption.

A. Onset of Non-Linear Stage
In the context of vacuum metastability, a natural concern is that particle-production processes arising from scale-invariant terms do not terminate in the linear stage of preheating. That is, these processes terminate only once field fluctuations become so large that they significantly backreact on the system at some time x NL . The longer the linear stage lasts, the more concern for vacuum destabilization since the Higgs fluctuations must remain sufficiently controlled over this entire duration.
In our model, either the Higgs or inflaton fluctuations may bring about an end to the linear stage, but a spectator field could potentially play this role as well (see Sec. V for further discussion along these lines). Let us first examine the necessary conditions for the Higgs to play this role. There are two competing effects that influence the background inflaton field [19]. On the one hand, as the inflaton loses energy to particle production, ϕ falls and thereby decreases the effective frequency of oscillations λ φ φ. On the other hand, as the Higgs variance grows the inflaton obtains an effective mass g 2 h 2 , increasing the oscillation frequency. The inflaton oscillations are disrupted once these two quantities are similar in magnitude, i.e., once g 2 h 2 /(λ φ φ 2 ) approaches unity. Of course, throughout this process the size of h 2 must remain controlled h 2 h 2 crit so that we do not destabilize the Higgs. According to Eq. (3.31), this requirement leads to a rather tight constraint on the couplings: which is approximately g 2 /λ φ 10 4 for our benchmark parameter choice λ φ = 10 −10 . However, we have already found that in this regime the rate of perturbative decays dominates over that of particle production [see previous discussion regarding Eq. (3.36)]. We can thus conclude that, in the absence of backreaction sourced by any other fields, the linear stage of preheating must be terminated by the inflaton fluctuations and not the Higgs. Unfortunately, the inflaton fluctuations grow slowly with the weak characteristic exponent µ max ≈ 0.036 (refer to the end of Sec. III D), so this can take a significant amount of time. The analysis above applies just as well to the inflaton variance φ 2 with the replacements g 2 → 3λ φ and Γ h → Γ φ , where Γ φ is the modeldependent inflaton decay rate. Then, the onset of the non-linear (NL) stage occurs at where W −1 denotes the negative branch of the Lambert W -function and we have neglected Γ φ . Using the parameters above we find x NL ≈ 413. Remarkably, this figure is consistent to less than 2% error with lattice-simulation results [63] that give x NL = 76 − 14.3 log λ φ ≈ 405.
vacuum barrier x end 10 20 30 x ≡ λ φ ϕη If the variance grows to the extent that ω 2 H k is dominated by the (negative) 3λ h H 2 term, the Higgs experiences a runaway tachyonic growth that rapidly destabilizes the electroweak vacuum, causing the sharp divergence in the curves. A point marker at the end of a curve indicates destabilization, and the gray curves show H 2 in the absence of decays. Identically to Fig. 4, the top (bottom) panels correspond to broad (narrow) parametric resonance, respectively. In the latter, ϕ 2 grows sufficiently to end the linear preheating stage at xNL ≈ 400, consistent with lattice results in the literature [63].

B. Vacuum Destabilization
In the above, we have established a clear criterion for model viability: for a given choice of couplings, if the electroweak vacuum survives for a duration longer than x NL , then it survives the linear stage of preheating. The transient non-linear stage that follows tends to asymptote to thermal equilibrium, during which resonant particle production stops and the energy density is redistributed among the modes to approach a thermal distribution. We expect that if metastability is maintained until x NL then it also survives the non-linear stage; we briefly discuss these issues and how massless preheating contrasts with other scenarios in the conclusions in Sec. V.
Our task then shifts to calculating the vacuum decay times x dec as a function of the model parameters. A robust method to conservatively estimate x dec is as follows [21,23]. As the Higgs fluctuations grow, the effective mass term associated with the Higgs self-interaction 3λ h h 2 becomes negative and large enough to make some Higgs modes tachyonic over a time interval ∆t. These modes grow exponentially in this interval at a rate 3|λ h | h 2 . In turn, this amplifies h 2 , which then makes the tachyonic masses even larger-i.e., the system enters a positive feedback loop which rapidly destabilizes the vacuum. We can estimate x dec as the time at which the growth exponent becomes larger than O(1): and we have used that ∆t = ∆x/( λ φ φ). Alternatively, in situations where the Higgs modes oscillate slowly relative to the background inflaton field, the simpler condition 3λ h h 2 m 2 h is appropriate, in which we evaluate the effective Higgs mass m h [of Eq. (3.33)] at the amplitude of the inflaton oscillations φ = φ.
The process of vacuum destabilization as it unfolds is illustrated in Fig. 8 by including backreaction in our numerical calculations. The panels reflect the same parameter choices made for Fig. 4 in the previous section, where we focus on the broad g 2 /λ φ = 2n 2 and narrow g 2 /λ φ = 2n 2 + n parametric resonances within the same coupling band (with n = 12). The curvature coupling ξ = 30 is subdominant but non-negligible. Both the Higgs and inflaton variance are shown, and the former shows the results both with (blue) and without (gray) Higgs decays. In the top panel, the perturbative decays delay destabilization but it ultimately occurs at x dec ≈ 35. The runaway growth of the Higgs fluctuations is observed toward the endpoint of each curve, as they finally cross the vacuum barrier H 2 ≥ (ah crit ) 2 . By contrast, in the lower panel we see that the narrow resonance would lead to destabilization in the absence of decays, after not much longer x dec ≈ 55. However, we see that accounting for perturbative decays prevents vacuum destabilization. Moreover, we observe the onset of the non-linear stage at x NL induced by the growth of inflaton fluctuations, as ϕ 2 grows sufficiently large. As expected, this stage occurs at x NL ≈ 400.
Let us now survey x dec over the parameter space of couplings and produce analytical estimates in each region. We shall work largely in parallel to Sec. III and employ the results from that section. That is, we shall first cover the ξ = 0 and g = 0 limits separately and then progress to the mixed case, providing a general picture of vacuum metastability at the end of the section.

From Parametric Instability
Let us first consider destabilization of the vacuum from the parametric resonance. The effective mass of a Higgs mode ω 2 h k ≈ k 2 /a 2 + g 2 φ 2 + 3λ h h 2 may first become tachyonic once the inflaton field passes through φ = 0. We can estimate x dec by following the discussion surrounding Eq. (4.8), but we must first calculate the Higgs variance. Using that h 2 n h /(a 3 ω h ), with Eq. (3.10) and Eq. (3.30) for the different coupling regimes, we find , (4.9) where we have defined the effective rate in terms of the conformal decay rate of Eq. (3.35). The Lorentz factor γ max(1, λ φ /g 2 ) accounts for the dilation of relativistic Higgs decays, which is negligible in the g 2 /λ φ 1 regime but important for g 2 /λ φ 1. Finally, the vacuum decay time is found in both regimes: (4.11) The g 2 /λ φ 1 result is based on the destabilization condition in Eq. (4.8). Meanwhile, for g 2 /λ φ 1 the Higgs modes oscillate slowly relative to the inflaton background, so that 3λ h h 2 g 2 φ 2 is an appropriate condition. Additionally, note that µ max depends on the regime of the resonance, as given explicitly in Sec. III A.
Manifestly, the effective rate µ h shows the competition between particle production and decays. For large enough coupling, the decays can eventually dominate, because µ max is globally bounded from above; the decay rate of produced Higgs particles scales with the coupling, but the rate of particle production does not.
Let us compute the critical value of the coupling for which x dec ≥ x NL , i.e., the smallest coupling for which the metastable vacuum survives preheating. We shall separately consider the couplings corresponding to the broad and narrow resonances. In the broad regime, we find g 2 /λ φ 2.4 × 10 3 . Likewise, using Eq. (3.9) for µ max , for the narrow resonances we find g 2 /λ φ 213. The dramatic gap between these thresholds reflects the difference in the strength of the resonance between these two regimes. Indeed, confirming our earlier observations In particular, the two green curves for g 2 /λ φ 1 show the broad and narrow regimes, evaluated using µmax ≈ 0.24 and Eq. (3.9), respectively. As these µmax correspond to the minimum and maximum growth exponents, they estimate the envelope of the highly non-monotonic x dec . In the right panel, the analytical estimates are also shown by green curves and correspond to Eq. (4.13) and Eq. (4.14), plotted over their regions of validity.
in Eq. (3.36), we find that perturbative decays of the Higgs stabilize the electroweak vacuum for sufficiently large quartic couplings, and these couplings significantly differ based on the regime of the parametric resonance.
To paint a more complete picture and confirm our analytical results, we solve the mode equations numerically within the one-loop Hartree approximation and scan continuously over the range of couplings g 2 /λ φ ; this is shown in the left-hand panel of Fig. 9. The black and gray curves show the numerical results with and without perturbative decays, respectively. In line with our calculations in Eq. (4.11), we find that x dec is highly nonmonotonic with respect to the quartic coupling. Indeed, this behavior originates from the structure of the resonance bands (in Fig. 2), in which the minima (maxima) of the oscillations correspond to the broad (narrow) resonance regimes, respectively. Therefore, evaluating our analytical result in Eq. (4.11) at the broad and narrow µ max values estimates the envelope of the x dec oscillations. This envelope is shown by the thick green curves at g 2 /λ φ ≥ 1 in Fig. 9. We also show our analytical estimate for the g 2 /λ φ ≤ 1 region. Meanwhile, the pivotal role of the perturbative decays is highlighted by the gray curve, which shows the numerical x dec with perturbative decays turned off. Indeed, in the absence of decays, the vacuum is stabilized only for very small coupling g 2 /λ φ 0.25.

From Tachyonic Instability
Likewise, let us consider the pure curvature coupling limit, in which g = 0. We once again confine our discussion to ξ > 0 since otherwise the vacuum decays. From Sec. III, we know that the growth of fluctuations from this tachyonic source is very different from the paramet-ric instability. Most notably, the tachyonic production ceases once x x ξ , even for the zero mode. We estimate the variance again using h 2 n h /(a 3 ω h ), and since the effective mass is curvature dominated we have ω 2 h k ξR ξλ φ φ 4 . Using Eq. (3.16) then yields for x < x ξ . Afterward, there is no particle production, so n h is fixed and ω h ∝ φ 2 ∝ 1/a 2 . The variance then redshifts as h 2 ≈ n h /(a 3 ω h ) ∝ 1/a, which is slower than the scaling one would find in the quartic-dominated case. The fact that h 2 redshifts in this way is a crucial observation. If the ratio h 2 /h 2 crit grows at times x x ξ due to its redshifting behavior, the vacuum may decay even after particle production shuts down. For the quartic interaction, this ratio would be fixed. However, for the curvature interaction, the variance redshifts more slowly h 2 ∝ 1/a, and the barrier redshifts more rapidly h 2 crit ξ h R/|λ h | ∝ 1/a 4 . Consequently, the ratio redshifts as h 2 /h 2 crit ∝ a 3 and grows. Even if the Higgs remains stable throughout the tachyonic phase, it destabilizes once the fluctuations grow beyond the extent of the barrier h 2 h 2 crit ; this leads to a vacuum decay at which is valid for x dec > x ξ , or equivalently ξ 16. Note that the presence of even a small quartic coupling may shield the vacuum from this effect, since the quartic term will eventually dominate g 2 φ 2 ξR and halt the growth of h 2 /h 2 crit , which happens once x 12ξλ φ /g 2 .
Otherwise, destabilization occurs during the tachyonic phase and the condition in Eq. (4.8) can be used again to produce an estimate. This task is now more complicated since backreaction is not the only source of tachyonicity. The vacuum could decay in one of two ways: (i) by the 3λ h h 2 term or (ii) directly by the tachyonic mass of the curvature term. A coupling of order ξ O(10) is sufficient to produce the latter, and in this case the vacuum decays rapidly. For smaller ξ, however, the false vacuum can survive much longer. Using Eq. (4.8) with the approximation that ∆x T /2 leads to (4.14) We consider this result consistent with our assumptions only if Eq. (4.14) gives an earlier decay time than Eq. (4.13); in terms of the curvature coupling, this implies a region of validity ξ 4 for Eq. (4.14).
Our numerical computations of the vacuum decay time x dec are shown in the right panel in Fig. 9. The pure curvature coupling case corresponds to the black curve and the analytical estimates constitute the thick green curve. We have also included the numerical results for several small quartic couplings in the range 10 −3 g 2 /λ φ 10 −1 . These results are in line with our rough estimates and confirm that the presence even of a small quartic coupling can prevent decay of the false vacuum. We have neglected ξ O(10) in our analytical estimates since these correspond to rapid destabilization, for which the precise values of x dec are not essential.

The Mixed Case
Let us finally consider the general mixed case in which both couplings ξ and g 2 /λ φ are non-zero. As we have observed in Sec. III C, if the strength of the curvature coupling is at least comparable to g 2 /λ φ the system undergoes a sequence of distinct phases of particle production. In the context of vacuum stability, this means that the electroweak vacuum must initially survive a phase of tachyonic production until x ξ [6ξ/ g 2 /λ φ ] 1/2 and then survive the parametric instability until non-linear dynamics set in at x NL . In what follows below, we shall assume such a scenario for our analytical calculations.
Rather than concern ourselves directly with x dec , we produce an estimate of the constraint for vacuum metastability by examining where h 2 h 2 crit . In the ξ > 0 region we can write h 2 using Eq. (3.23) together with h 2 n h /(a 3 ω h ). We arrive at the expression Owing to the exponential growth of h 2 , the constraint on our parameter space has only a logarithmic sensitivity to h 2 crit and the coefficient in Eq. (4.15). That said, the constraint is sensitive to the rate of perturbative decays, and we should reintroduce this rate as we calculate our result. Along these lines, we find where all of the logarithmic terms have been collected into the quantity log C + . The dependence on x dec and aΓ h are important. To manage x dec , a natural assumption is that x dec x ξ in the region we are evaluating, since otherwise the tachyonic effect is not strong enough to destabilize the vacuum. We therefore have at most x dec (6ξ/ g 2 /λ φ ) 1/2 . We can solve the more complicated inequality that results for ξ (and we shall do this numerically below), but based on our observations of the tachyonic effect in this section x dec = O(10) serves as an appropriate simplifying approximation. As discussed in Sec. III D, the perturbative decays arising from the quartic interaction are generally more relevant, and these scale parametrically as Γ h ∝ g 2 /λ φ . These considerations lead to a constraint of roughly ξ O(0.1)g 2 /λ φ to ensure metastability. An interesting feature of the mixed-coupling scenario is that it opens the possibility of metastability in the ξ < 0 region, which is ruled out if g 2 /λ φ = 0. A similar procedure to the above is followed to produce a bound on the ξ < 0 region, the details of which are included in Appendix B. The constraint is expressed as and we have again taken x dec = O(10) and packaged the logarithmic terms in a quantity C − . The function giving H 1 (3, 0) ≈ 0.76 is defined in Appendix B. Note that for ξ 1 the term proportional to g 2 /λ φ on the left-hand side is subdominant and the remaining terms balance for |ξ| ∼ g 2 /λ φ . Therefore, neglecting numerical coefficients, the constraint becomes ξ O(g 2 /λ φ ). The negative curvature coupling thus yields a weaker bound, which is expected based on our observations from Sec. III C. Note that the inflation constraint ξ h − ξ φ g 2 /λ φ 1/12 in Eq. (2.10) yields a bound that is roughly coincident.
A more complete picture of the constraints is achieved by numerical calculations, and we display these in Fig. 10. For this figure, we have calculated the time of vacuum decay x dec over the {g 2 /λ φ , ξ} parameter space, making the simplifying assumption that ξ φ = 0. Note that the gray region excludes the subset of couplings that destabilizes the vacuum during inflation. Additionally, note that the curvature coupling is shown on a log scale over the negative and positive axes, with the exception of the range −1 ≤ ξ ≤ +1 where the scale is linear. The numerics run until a time x NL so that the white regions effectively show where x dec /x NL ≥ 1 and thus where the vacuum remains metastable throughout preheating.  21  22  23  24  25  26  28  29  30  31  32  33  34 x dec x NL FIG. 10. The decay time x dec of the electroweak vacuum computed numerically and normalized by xNL. The gray region is excluded due to violating Eq. (2.10). The white regions show "islands of (meta)stability" in which the false vacuum survives preheating. Along the g 2 /λ φ direction, the pattern reflects the band structure in Sec. III A, with the least stable regions centered around g 2 /λ φ = 2n 2 (for n ∈ N). For larger g 2 /λ φ , perturbative Higgs decays enlarge the metastable regions until they become contiguous at g 2 /λ φ 2 × 10 3 [in agreement with Eq. (4.11)]. The effect of the curvature interaction is to form an envelope over the metastable regions, and the yellow curve shows our estimate [from Eq. (4.16)] for ξ > 0.
Echoing the behavior observed in Fig. 9, we find that over the vast majority of parameter space, the metastable regions are not connected. Instead, we find a large number of disjoint "islands of (meta)stability" scattered throughout. Indeed, the constraint for metastability cannot be fully expressed as a simple bound on the couplings-the fate of the electroweak vacuum depends on g 2 /λ φ in a highly monotonic way. Naturally, the most unstable regions appear where the characteristic exponent µ max is maximized-i.e., around the broad resonances g 2 /λ φ = 2n 2 (for n ∈ N). Conversely, the most stable regions appear around the narrow resonances, where µ max is minimized [see Eq. (3.9)]. As g 2 /λ φ is taken to larger values, the metastable regions grow and start to merge, as shown by the magnified inset panel. The integers written over the metastable regions correspond to the couplings of the narrow resonances written in the form g 2 /λ φ = 2n 2 + n. The couplings eventually become large enough for perturbative decays to stabilize the Higgs, re-gardless of the resonance regime, forming a contiguous metastable region for g 2 /λ φ 2 × 10 3 . This observation agrees with our analytical estimates in Eq. (3.36) and Eq. (4.11) for the broad regime. In principle, the only limitation in taking larger couplings is that we do not ruin the flatness of the inflaton potential, which (as briefly discussed in Sec. II A) roughly requires g 2 /λ φ 10 5 . Beyond these observations, we notice that the interplay of the couplings is also clarified in Fig. 10. In particular, the figure shows that the Higgs-curvature interaction has the effect of cutting off the metastable regions and imposing an envelope that scales with the couplings. This envelope is what we have analytically estimated in Eq. (4.16) and Eq. (4.17). We have plotted the numerical solution to the inequality in Eq. (4.16) using a yellowdashed curve. This curve is indeed consistent with the approximate bound that we estimated ξ O(0.1)g 2 /λ φ . We have not included the analogous curve for the ξ < 0 region since this closely coincides with the constraint (the gray region) on the stability of the vacuum during inflation. On the whole, the interplay between the two interactions effectively extends the range of metastability for the Higgs-curvature coupling, shielded by the presence of a similarly large Higgs-inflaton coupling.

V. CONCLUSIONS AND DISCUSSION
Some amount of degeneracy often exists in inflationary models in the sense that observational constraints may be satisfied over a degenerate subspace of the model parameters. Undoubtedly, one reason that studies of postinflationary dynamics are essential is that this dynamics may break such degeneracies by leading to qualitatively different preheating histories, subject to a qualitatively different set of possible constraints. For example, as we have encountered in this paper, the stability of the Higgs during inflation can be ensured by either a direct coupling to the inflaton or a non-minimal coupling to gravity. The dynamical roles played by these interactions are similar during inflation but differ dramatically during preheating, in which the non-trivial interplay between these interactions implies a rich structure of metastable regions and associated constraints on the Higgs couplings. Indeed, the broader motivation for this work is to explore dynamics that may reveal independent probes for models of early-universe cosmology.
In this paper, we have examined the massless preheating dynamics that arises in models composed of scaleinvariant interactions, in which the inflaton potential is effectively quartic after inflation. Most notably, we have focused on the implications these models have for electroweak vacuum metastability. We have provided constraints on the couplings for which the Higgs remains stabilized during and after inflation. Among these couplings, we have included the possibility that the Higgs and the inflaton have non-minimal couplings to gravity. While our study is motivated by addressing the metasta-bility of the vacuum, our results are readily generalized to other approximately scale-invariant models.
Several comments are in order. First, while we have considered non-minimal gravitational couplings for both the Higgs and inflaton fields, the Higgs coupling ξ h has received the bulk of our attention. And while the effects of ξ φ have been included in our analysis through the effective curvature coupling ξ [in Eq. (3.6)], our treatment can be extended in several ways. For example, we have not taken into account the higher-order effects on the field fluctuations that appear as a result of ξ φ = 0 corrections to the background-field solutions. For this reason, our results are applicable only for |ξ φ | 1, which is the scope assumed for this paper. These higher-order effects are challenging to incorporate analytically in massless preheating. Notably, even the zeroth-order background solution [the elliptic function in Eq. (2.15)] is considerably more complicated than the sinusoidal form found in massive preheating. It would be worthwhile to study the effects of ξ φ ≈ −O(1) without resorting to such approximations for the background equations of motion.
Second, our study is performed under the assumption that there is no new physics beyond the SM that significantly alters the renormalization-group evolution of the Higgs quartic coupling. This assumption could be reasonable, as there are no hints of new physics in collider or weakly interacting massive particle dark-matter searches. Either way, the cosmological history following the preheating stage-i.e., the reheating epoch and, more precisely, the inflaton decay channels-should be consistent with this assumption. The reheating epoch could then be realized in a number of ways. For instance, the inflaton could couple to the right-handed neutrino, which is responsible for neutrino masses and leptogenesis [64]. In this case, reheating could be completed by the decay of inflaton to right-handed neutrinos and their subsequent decay to SM particles. Alternatively, if the inflaton is stable, it also opens the possibility that the inflaton is a candidate for dark matter [61,65] or dark radiation [66], provided that preheating converts most of the inflaton energy density into SM degrees of freedom. Indeed, the details and variations on these possibilities ultimately depend on the sign and size of the inflaton mass-squared term.
Along another direction, some comments are in order regarding the stage of preheating that occurs after x NL in our study. After entering the non-linear stage of preheating, the energy density held in fluctuations is redistributed among the modes by rescattering processes. In general, the Higgs can destabilize during this time, and this possibility has been addressed for massive preheating in Ref. [21], based on the results of Refs. [67,68]. In massive preheating, this destabilization arises because the effective Higgs mass induced by the inflaton-Higgs coupling redshifts as g 2 φ 2 ∝ 1/a 3 , while that induced by the Higgs self-coupling redshifts as 3λ h h 2 ∝ 1/a 2 . In other words, the tachyonic contribution continues to grow relative to the stabilizing mass term and can eventually trigger de-cay of the vacuum, even though it was stabilized at the end of the linear preheating stage. The details of the thermal Higgs mass after preheating and the evolution of the background temperature thus become important to ensure metastability as thermalization begins. That said, in massless preheating, these two mass contributions redshift at the same rate g 2 φ 2 ∼ 3λ h h 2 ∼ 1/a 2 , so these same concerns over destabilization are not as relevant. Nevertheless, the stages of evolution leading to thermalization are, of course, interesting in the context of massless preheating for a host of other reasons. Furthermore, other considerations such as the thermal corrections to the vacuum tunneling probability become important during the reheating epoch, such as those studied in Ref. [69].
Thermal effects could be important even before x NL , through direct decay channels (if any) of the inflaton into SM particles or through the Higgs decay products. Whether the Higgs decay products become quickly thermalized is a nontrivial issue, since the particle-production rate can greatly exceed the thermalization rate. In this case, one may not necessarily assume that the Higgs decay products form a thermal plasma, even if the thermal relaxation time is much shorter than the age of the Universe at that moment. Nevertheless, if instant thermalization of the decay products is assumed, we can discuss thermal backreaction on the Higgs dynamics. If the thermal mass of the Higgs becomes comparable to the zerotemperature effective mass, our analysis would be modified, and we would expect broader regions in Fig. 10 to become (meta)stable. Our preliminary study shows that the effect of the Higgs thermal mass is insignificant over a wide range of the parameter space. However, it is important to correctly account for the full evolution of the decay products without the thermalization assumption, which we leave for future work.
As we recall from Sec. IV A, the time x NL at which the non-linear stage of preheating begins is central to our results. We concluded, in our minimal model realization, that the inflaton self-resonance is responsible for the onset of the non-linear stage, as its fluctuations grow from the quartic self-resonance. However, a natural extension is to consider the presence of a spectator field that does not couple to the Higgs. If such a spectator field is likewise subject to particle production, then x NL could be triggered by the growth of the spectator fluctuations instead; this is not difficult to arrange, considering that the characteristic exponent for the inflaton self-resonance is rather small at µ max ≈ 0.036. A reduction in x NL would imply that our metastable regions, such as those that appear in Fig. 10, grow in extent. This enhancement could be substantial. For example, taking a quartic coupling between the spectator and the inflaton, with the maximal exponent µ max ≈ 0.24, we find that x NL ≈ 63, roughly a factor of 6 smaller than the inflaton self-resonance.
Moreover, the possibility of spectator fields opens other avenues of exploration. For instance, such fields could generate important cosmological observables such as non-Gaussianities in the density perturbations [70]. Additionally, spectator fields could naturally serve as dark-matter candidates. The corresponding relic abundance generated from preheating would depend non-trivially on the couplings and spins of the spectator fields. It is worth noting that such dark-matter production from the thermal bath with a generic equation of state-including the radiation-like equation of state relevant to our studyhas been a topic of recent interest [71,72].
There are also a number of possible extensions that would be interesting to explore in the context of multifield inflation. For instance, in models which allow for a sizeable angular inflaton velocity at the end of inflation [73,74], the inflaton-dependent modulation of the effective Higgs masses could be suppressed. As a result, the rate of particle production could be suppressed as well. On the other hand, the onset of the non-linear stage of preheating would also occur at a later time x NL . The end result of the competition between these two effects, and the more complicated preheating dynamics altogether, requires a dedicated study [75]. An alternative multi-field extension could involve hybrid inflation [76]. In this case, vacuum stability during inflation is ensured by the inflaton-Higgs coupling as usual, but this coupling will not oscillate much during the preheating phase since most of the energy is transferred to the so-called waterfall fields [77]. These features imply that the electroweak vacuum would be relatively stable during both inflation and preheating.
It would also be interesting to consider the presence of additional terms that break the scale invariance of the theory, as these can have rich dynamical implications. While the breaking of the scale invariance due to the nonminimal coupling terms are relevant only around the end of inflation, effects of lower-dimensional terms become more and more important at a later time. For instance, a non-vanishing inflaton mass m φ might alter the parametric resonance once the inflaton oscillations reach a sufficiently small amplitude [19]. In particular, for couplings g 2 /λ φ λ φ /m 2 φ the resonance becomes stochastic once φ m φ / λ φ . The parametric resonance then increasingly behaves as in massive preheating (see Appendix A), with the stochastic resonance giving way to the narrow resonance onceφ m φ /g. Not only does this evolution change the spectrum of produced particles and the production rates, the introduction of m φ = 0 also allows the resonance to terminate at some time before x end . A small inflaton mass can thus modify our constraints for electroweak metastability in a non-trivial way.
Yet another possible extension of our work concerns the gravity formulation. In this paper, given the nonminimal couplings of our scalar fields to gravity, we have observed some of the physical distinctions that may appear between the metric and Palatini formulations. Still, a more general extension could be explored along the lines of the Einstein-Cartan formulation in the presence of Holst and Nieh-Yan terms [78][79][80]. The Einstein-Cartan theory generalizes our treatment, reproducing the metric and Palatini formulations in certain limits and allowing for a continuous interpolation between them.
All in all, if the metastability of the electroweak vacuum may be ensured in models from which massless preheating emerges, we may adopt a renewed interest in observables that are sensitive to the details of preheating. A well-suited example consists of gravitational waves (GWs) produced during this epoch [81]. Massless preheating is distinctive in this sense because the amplitude of the produced gravitational radiation does not dissipate and the frequency of the waves does not redshift. While these properties may present challenges for the observational prospects of GWs from such scenarios, several methods show promise [82]. For example, radio telescopes may be used to examine the distortion in the CMB from by-products of high-frequency GWs interacting with background magnetic fields [83]. Indeed, such observations could even be used to probe broader early-universe properties, such as the energy scale for inflation [84].
The inflaton amplitude then redshifts as φ ∝ 1/a 3/2 and [in line with Eq. (2.17)] the cosmological equation of state is matter-like. The inflaton field evolution is sinusoidal: where in analogy to Sec. II B we have defined ϕ = a 3/2 φ and the dimensionless time x ≡ mt. The constant x 0 is fixed by the conditions at the end of inflation, and we shall assume φ end = O(1). The scale factor evolves as and thus for consistency x end = 2 2/3φ −3/2 end . To compute the growth of fluctuations we write the equations of motion for the Higgs modes in the Einstein frame. By defining H k ≡ a 3/2 Ω −1/2 | h=0 h k , we can remove the damping terms, and up to O(a −3 ) we have in which the effective masses are We have defined the rescaled momenta κ ≡ k/m, and in the second line we have employed Eq. (A3). Additionally, in analogy to Eq. (3.6) we have defined the effective curvature coupling The rate of cosmological expansion is much smaller than that of the inflaton oscillations, so that the mode equations in Eq. (A5) approximately take the Mathieu form H k + {A k − 2q cos[2(x − x 0 )]}H k = 0, with the slowly varying parameters given by An important observation is that the contributions to the effective Higgs mass from the quartic interaction g 2 φ 2 ∝ 1/a 3 and non-minimal coupling ξ h R ∝ 1/a 3 redshift identically in the massive preheating scenario. Because much of the complexity we found in our study of massless preheating arose from the mismatch between these two terms, the massive preheating scenario is simpler in this respect. For example, whether particle production for a given mode occurs through the tachyonic or parametric instability is determined by the sign of A k − 2|q|. Namely, A k − 2|q| > 0 yields the parametric instability, and A k − 2|q| < 0 yields the tachyonic instability. We can write these explicitly as Therefore, unlike in the massless preheating scenario, the type of particle production that drives the zero mode is a fixed time-independent property based on the couplings. Notably, in the g, ξ φ → 0 limit we recover the well-known result that 3/16 < ξ h < 3/8 is necessary for the complete absence of tachyonic production [21]. Another crucial distinction between the massless and massive preheating dynamics is found by reexamining the perturbative Higgs decays. In particular, using the restframe decay rate Γ h in Eq. (3.32), the decays introduce a dissipative exponent − dx Γ h /γ that suppresses the growth of the Higgs number density. In terms of the Mathieu parameters the effective Higgs mass is given by where the integration region is over the preheating times x ≥ x end for which m 2 h ≥ 0 and we have taken γ = 1. Because this integral has an identical form to the accumulated phase of the zero mode Θ 0 [refer to the general quantity defined in Eq. (3.12)], we can employ previous calculations for Θ k available in the literature [16]. The decay exponent amounts to a sum over each nontachyonic oscillation period − dxΓ h ∝ dx Θ 0 , and we can extract the time dependence by noting that the integrand scales approximately as Θ 0 ∝ |q| ∝ 1/x so that The logarithmic time-dependence of the decay exponent in Eq. (A11) demonstrates a critical distinction between the massive and massless preheating scenarios. For the latter, the time-dependence of the decay exponent was found [in Eq. (3.34) and Eq. (3.35)] to be linear, making it possible for the perturbative decays to efficiently counter the growth of Higgs fluctuations, even shutting them down for sufficiently large quartic coupling. By contrast, the logarithmic time-dependence in massive preheating shows that the capacity for perturbative decays to stabilize the electroweak vacuum is relatively negligible. This distinction makes the features of the massless preheating dynamics somewhat unique.
In order to further illustrate this comparison, we have plotted the evolution of the Higgs variance H 2 in Fig. 11 for both the quadratic and quartic scenarios, using the quartic coupling g 2 /m 2 φ = g 2 /λ φ = 800 and curvature couplings ξ h = ±200. The inflaton mass m φ = 10 −5 and coupling λ φ = 10 −10 are taken in the respective cases.  11. A comparison of the Higgs variance H 2 between massless and massive preheating, i.e., between the quartic and quadratic inflaton potentials, respectively, taking m φ = 10 −5 and λ φ = 10 −10 . The backreaction effect is not included in this figure. The top and bottom panels show a difference in sign for the curvature coupling ξ h . We have also included curves that show the result if perturbative decays are turned off. Notably, this effect is minimal in massive preheating but significant in the massless case. The contrast between the preheating scenarios is also evident in the short-lived growth of fluctuations for the massive case which ceases for |q| 1 4 . Note that although x and H have the same qualitative meaning in both preheating scenarios their explicit definitions differ due to the difference in background cosmology.
In contrast to the unimpeded growth that arises from the quartic potential, we observe the end of particle production that occurs in the quadratic case once |q| 1/4. And comparing the curves, which show the results with perturbative decays both included and not included, we see manifestly the relative importance of perturbative decays in these scenarios. Indeed, although the couplings are the same for each of the comparisons in Fig. 11, the results for vacuum stability can greatly differ. For instance, upon including backreaction one finds that, for the couplings in the bottom panel, the false vacuum survives in the quadratic case, but it decays in the quartic case. Notwithstanding, owing to the complex pattern of metastability regions for the quartic case, we can change this outcome with only a small change in the coupling.
An analytical study of the preheating dynamics and vacuum decay times is found in the literature for the un-  Fig. 10).
The results are given relative to a fiducial time x * ≈ 2 × 10 3 for which the resonance is deep in the narrow regime over the parameter space shown. Our most immediate observation is that the metastable region is connected, which contrasts with the more complex disjoint regions found in the massless preheating scenario. The features appear consistent with the bounds found in previous studies in the literature [21,24].
mixed case [21], while the mixed case is studied mostly numerically in Ref. [24]. For the purpose of our comparison we shall focus on numerical results for x dec . Along these lines, we provide an analog to Fig. 10 for the massive preheating scenario in Fig. 12. Note that in order to make a clear comparison between this result and the massless preheating result in Fig. 10, we have constructed Fig. 12 using the same axes and axis scales. The gray exclusion regions correspond either to the Higgs destabilizing during inflation ξ h −g 2 /(2m 2 φ ) or the flatness of the inflaton potential being ruined by quantum corrections ξ h 1 2 (10 −6 − g 2 )/m 2 φ . The line showing q = 0 is significant since along it the effective masses of the Higgs modes do not oscillate in time and the resonant particle production vanishes.
The most immediate observation we make in Fig. 12 is that the structure of the metastable region is relatively simple in comparison to Fig. 10. For g 2 /m 2 φ 10 3 , the upper bound of the metastable region is roughly linear but flattens for smaller values of g 2 /m 2 φ . These results appear consistent with both Ref. [21] and Ref. [24]. In particular, in Ref. [21], bounds on the unmixed couplings were given approximately as g 2 /m 2 10 3 and ξ 10, as reflected in our figure. While the first approximation gives a more accurate fit, the second is linear in both α and β and useful for our estimates. We have plotted the latter approximation using solid curves in the right panel of Fig. 13, with point markers indicating the numerically integrated result. Using this approximation, we find that the growth exponent is given up to O(κ 2 ) by where j is the number of times the effective Higgs mass has passed through the tachyonic region. Note that Eq. (3.25) is recovered in the κ → 0 limit. The growth exponent readily gives the phase-space density n h k via Eq. (B1) if we neglect the oscillatory component Θ k .
Using the saddlepoint method to integrate over the momenta, we obtain the Higgs comoving number density which is valid for times x x ξ . Using that the Higgs variance is h 2 n h /(a 3 ω h ), we are finally in a position where we can produce the constraint in Eq. (4.17).

Appendix C: Numerical Methods
In this appendix, we summarize the numerical methods used to compute the evolution of dynamical quantities in this paper. Note that analogous methods are applied for the comparison to massive preheating in Appendix A.
We perform our computations in the cosmological time t and establish our initial conditions at the end of inflation. According to the analysis in Sec. II A, the initial conditions at t end are given by where φ end = 1 is assumed. Furthermore, we assume for our numerics that the inflaton potential is given by a purely quartic potential V (φ) = λ φ φ 4 /4, and we have neglected the difference between the canonical and noncanonical inflaton field. While this affects the initial velocity by an O(1) factor, we have checked that such changes have a negligible impact on our results. The inflationary constraint in Eq. (2.10) is assumed so that the homogeneous value h for the Higgs field is stabilized strongly at the origin. Even after the end of inflation, accounting for h has a negligible influence on our results, so we ignore its equation of motion overall. We solve the equations of motion for the background inflaton field φ, inflaton fluctuations φ k , and Higgs fluctuations h k . As discussed in Sec. IV, the backreaction is accounted for using the Hartree approximation h 4 → 6 h 2 h 2 − 3 h 2 2 for the Higgs self-interaction and using an analogous substitution for the quartic inflaton potential φ 4 → 6 φ 2 φ 2 − 3 φ 2 2 [62].
We discretize the momenta linearly over the range k ∈ [0, Λ] with N = 300 lattice points. The momentum cutoff Λ is chosen based on the structure of the parametric and tachyonic instabilities; i.e., we use the maximal momentum of each unstable band. For example, with ξ h = 0 and g 2 /λ φ 1 our numerics use a cutoff near the maximum of the resonance band Λ = λ φ ϕ(g 2 /λ φ ) 1/4 . The equations of motion are then given bÿ h k + (3H + Γ h k )ḣ k + ω 2 h k h k = 0 , (C4) where the mode energies are The Hubble parameter is computed using the first Friedmann equation H 2 = ρ tot /3 for the total energy density ρ tot . Note that ρ tot includes both the energy density of the inflaton background ρ φ =φ 2 /2 + V (φ) and that of the field fluctuations, although the latter do not provide a significant contribution until the system approaches the non-linear stage of preheating. We evolve these equations of motion until the Higgs destabilizes or the onset of non-linear dynamics is reached at x NL ≈ 400.
The decay rates Γ h k for the Higgs modes are calculated in the same way as Eq. (3.32), except in our numerical calculations we generally write Γ h k = Γ h /γ, where the Lorentz factor γ carries some momentum dependence.
The initial conditions for the field fluctuations in Eq. (C3) and Eq. (C4) are taken as Gaussian such that φ k = 1/ 2ω φ k and h k = 1/ √ 2ω h k at the initial time.
However, this initial condition is UV divergent and sensitive to the cutoff scale Λ introduced in our discretization. Employing adiabatic renormalization [85], we calculate the variance so that it is initially vanishing: The scaling of the counterterm with a is chosen so that it scales with the physical momentum k/a. Finally, we incorporate the running of the Higgs fourpoint coupling λ h (µ) by evaluating the renormalization scale at µ = h 2 . The effect of the backreaction is most relevant at energies near h SM ≈ 10 10 GeV where, based on SM measurements, the coupling changes sign. Therefore, in our numerics we take the approximate which is similar to that used in other numerical studies [21,23] of electroweak vacuum metastability.