Large Density Perturbations from Reheating to Standard Model particles due to the Dynamics of the Higgs Boson during Inflation

Cosmic Microwave Background (CMB) observations are used to constrain reheating to Standard Model (SM) particles after a period of inflation. As a light spectator field, the SM Higgs boson acquires large field values from its quantum fluctuations during inflation, gives masses to SM particles that vary from one Hubble patch to another, and thereby produces large density fluctuations. We consider both perturbative and resonant decay of the inflaton to SM particles. For the case of perturbative decay from coherent oscillations of the inflaton after high scale inflation, we find strong constraints on the reheat temperature for the inflaton decay into heavy SM particles. For the case of resonant particle production (preheating) to (Higgsed) SM gauge bosons, we find temperature fluctuations larger than observed in the CMB for a range of gauge coupling that includes those found in the SM and conclude that such preheating cannot be the main source of reheating the Universe after inflation.

Cosmic Microwave Background (CMB) observations are used to constrain reheating to Standard Model (SM) particles after a period of inflation. As a light spectator field, the SM Higgs boson acquires large field values from its quantum fluctuations during inflation, gives masses to SM particles that vary from one Hubble patch to another, and thereby produces large density fluctuations. We consider both perturbative and resonant decay of the inflaton to SM particles. For the case of perturbative decay from coherent oscillations of the inflaton after high scale inflation, we find strong constraints on the reheat temperature for the inflaton decay into heavy SM particles. For the case of resonant particle production (preheating) to (Higgsed) SM gauge bosons, we find temperature fluctuations larger than observed in the CMB for a range of gauge coupling that includes those found in the SM and conclude that such preheating cannot be the main source of reheating the Universe after inflation.

I. INTRODUCTION
Inflation [1][2][3] is a period of accelerated expansion that occurred in the very early epoch of our Universe. It was first proposed to explain the homogeneity, isotropy, and flatness observed in the cosmic microwave background (CMB) radiation [4,5], as well as the lack of relic monopoles. A mechanism for driving the dynamics of inflation comes in the form of a rolling scalar field, the inflaton [6,7]. In this framework, the density perturbations that are observed in the CMB are explained by the quantum fluctuations of the inflaton field. It is these perturbations that later develop into the large scale structure observed in the Universe [8,9]. The most stringent constraints to the theory of inflation come from the observations of the CMB by the Planck satellite, including the power spectrum [10] and bispectrum [11] of temperature anisotropies. * Electronic address: aliki.litsa@fysik.su.se † Electronic address: ktfreese@utexas.edu ‡ Electronic address: e.sfakianakis@nikhef.nl § Electronic address: pstengel@sissa.it ¶ Electronic address: luca.visinelli@lnf.infn.it The inflationary period must end by successfully reheating the Universe, which marks the transition into the radiation dominated cosmological era before Big Bang Nucleosynthesis (BBN) occurs. If inflation is driven by a rolling scalar field, reheating can occur through the decay of the inflaton into light degrees of freedom in the Standard Model (SM) or an intermediate sector. Typical mechanisms for reheating are perturbative decay of the inflaton [12,13] and resonant particle production [14]. In particular, the inflaton field φ might have decayed due to the presence of interaction terms in the Lagrangian such asψψφ and φF µνF µν . The former term is a Yukawa-type coupling to a fermion ψ and the latter is a Chern-Simons coupling to gauge bosons with a field strength F µν , as found in models of natural inflation [15]. Other scenarios for reheating include the decay of the inflaton condensate into its own quanta, which must ultimately decay into SM particles, or through gravitational interactions [16].
In this paper we consider reheating via inflaton decay to SM particles coupled to the Higgs boson. We note that the inflaton in this paper is not the Higgs boson; instead the Higgs is a light spectator that plays an important role in the reheating process. The scenario we consider is minimal in the sense that no new particles beyond the SM are introduced other than the inflaton itself. Along with the inflaton flat direction which is the main compo-nent in driving the expansion rate of the Universe during the inflationary stage, the Higgs boson and other light fields that are present at this epoch would act as spectators since they would not directly affect the evolution of the background geometry. However, light spectator fields would acquire large quantum fluctuations and thereby effective masses that vary from one Hubble patch to another. If the light spectators are also associated with the decay of the inflaton field in each Hubble patch, their stochastic dynamics can cause spatial fluctuations in the reheat temperature and large density perturbations.
Inhomogeneous reheating due to the stochastic behavior of a light spectator field is known as modulated reheating [17][18][19][20][21]. Examples of light spectators can include the SM Higgs boson, with mass ∼ O (125) GeV [22,23] and the hypothetical axion, with a mass typically well below the MeV range (see Ref. [24] for a recent review). In this paper we focus on modulated reheating caused by coupling of the inflaton decay products to the SM Higgs boson [25][26][27][28], which is taken to be a light spectator during inflation. We assume that the inflaton couples primarily to SM particles that develop masses when the Higgs field acquires a Vacuum Expectation Value (VEV). If the inflaton were instead to decay to a massless gauge mode in the broken phase (the analog of the photon at lower temperatures) or to un-Higgsed neutrinos, then the effective Higgs mass during inflation would be irrelevant to the reheating process, since the SM Higgs does not couple to these particles 1 .
Enqvist et al. [29] showed that during inflation, due to the quantum fluctuations of the Higgs field, the Higgs can develop a mass so that electroweak (EW) symmetry can be treated as effectively broken [30,31]. The expectation value of the Higgs amplitude over the entire inflating patch is vanishing h I = 0 due to the symmetric potential (where the subscript I is used to indicate the initial value at the onset of inflaton oscillations, i.e. at the end of inflation). However, due to the quantum fluctuations, the variance is non-zero. Thus, the typical Higgs amplitude, i.e. the effective VEV, in a typical Hubble patch at the end of inflation is given by a root mean square value h I = h 2 ∝ H I where H I is the Hubble scale at the onset of inflaton oscillations. Even assuming that the energy density of the Higgs field is always subdominant to that of the inflaton, its effective VEV would be driven by the stochastic dynamics to relatively large values. The effective nonzero Higgs VEV during inflation then gives mass to all SM particles that couple to the Higgs. In turn, this would affect the decay of the inflaton to SM particles. Although usually considered to be much lighter than the inflaton, SM particles with masses due to the Higgs boson condensate which develops during inflation [29,32] can have several interesting effects.
In our previous paper (Ref. [33] hereafter referred to as FSSV), we studied the effects of Higgs blocking, i.e. the delay in the reheating process due to the large particle masses acquired during inflation due to the effective Higgs VEV. As long as the particle masses exceed the inflaton mass, reheating cannot occur. Only once the Higgs condensate decays do the particle masses vanish and reheating can proceed. We studied Higgs blocking for the cases of both perturbative decay and resonant particle production. We also briefly discussed the potential for generating large temperature fluctuations due to the stochastic nature of the Higgs blocking, with variation from one Hubble patch to another. Subsequently, Ref. [34] calculated signatures for the effects of the Higgs on reheating in the Cosmological Collider framework and Ref. [35] showed the potentially large temperature fluctuations which can arise due to Higgs blocking when the inflaton decays to fermions with relatively large Yukawa couplings to the SM Higgs. In our present work we extend our previous results of FSSV and treat both Higgs blocking and Higgs modulation to derive the corresponding density fluctuations that arise during reheating. For the case of perturbative inflaton decays, our previous work in FSSV found Higgs blocking to be negligible for SM fermions with Yukawa couplings y < 1; here we will show that, even for this case, Higgs-modulation alone can indeed generate large temperature anisotropies.
In this work we consider a simple reheating scenario where the inflaton (in the post-inflationary epoch) proceeds along a single field direction with a potential approximated by that of a massive scalar field. We do not further specify any particular inflationary potential. We must again stress that the inflaton field in this study is not the Higgs [36] field, which is taken to behave as a spectator field. We treat the decay rate of the inflaton Γ φ as a time-and space-dependent quantity because of its relation to the value of the Higgs field. Perturbations in the Higgs field lead to variations in Γ φ , similar to those discussed for example in Ref. [17]. We use the perturbed Einstein equations introduced by Ref. [17] to derive the temperature fluctuations induced by the spatial dependence of the Higgs field. While Ref. [35] has a similar parameterization of Higgs effects on reheating, the amplitudes of temperature fluctuations we calculate in this work are notably larger for equivalent choices of parameters and we are able to show that large temperature fluctuations are even possible in the case of perturbative reheating without the full effects of Higgs blocking. 2 Furthermore, in our framework it is straightforward to extend our analysis to the case of non-perturbative preheating, for which we use a similar set of equations to calculate the density fluctuations induced by the spatial dependence of (Higgsed) SM gauge boson masses.
In Sec. II A we present equations which describe the evolution of energy densities in different Hubble patches which are characterized by different inflaton decay rates in each patch. In Sec. II B we show the corresponding perturbed equations describing the evolution of adiabatic matter and metric perturbations after the end of inflation. After elaborating on the dynamics of the Higgs field in Sec. II C, we explain our methods for calculating the comoving curvature perturbation in the cases of perturbative and resonant inflaton decays, as well as its connection to the temperature anisotropies observed in the CMB in Secs. II D, II E and II F, respectively. Our numerical results are presented in Sec. III for perturbative inflaton decays and Sec. IV for the case of reheating via resonant decays of the inflaton into gauge bosons. We summarize our findings and offer our prospects for future work in Sec. V.
To guide the reader to our main results: Our primary results for the case of (spatially-dependent) perturbative inflaton decay to SM particles can be found in Fig. 9. This figure shows constraints on the parameters (inflaton decay rate to SM fermions when masses are negligible, Yukawa coupling of SM decay products with the SM Higgs boson, and self-coupling of the Higgs during inflation) obtained by requiring that the amplitude of temperature fluctuations do not exceed CMB observations. For the case of resonant preheating to SM Higgsed gauge bosons, our main results can be seen in Fig. 10; one can see that for all reasonable parameter choices the density fluctuations (shown in terms of the Bardeen potential) are too large.

A. Unperturbed equations
We work using a flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric, described by the line element where a(t) is the scale factor, x are spatial coordinates, and t is cosmic time. We assume that towards the end of the inflationary period, the inflaton field begins to oscillate about the minimum of its potential, behaving as a massive scalar field with energy density ρ φ . We consider a collection of n Hubble patches, each of which is characterized by a different inflaton decay rate Γ j φ , for j ∈ {1, 2, . . . , n}. The perturbative decay of the inflaton field into radiation at the rate Γ j φ in the j-th Hubble patch is described by the coupled first-order equations where ρ j φ , ρ j r are the respective energy densities of the inflaton and of radiation, both in the same patch. The independent variable N j in Eqs. (2)-(3) is the number of e-folds since the beginning of the reheating stage in patch j at time t i , defined as where the Hubble rate is H j =ȧ j /a j with a j the scalefactor. The Hubble rate is related to the total energy density at the patch we are considering through the Friedmann equation where G is Newton's constant. At this point it is also useful to define the background energy densities, with respect to which perturbations are calculated in Sec. II B. These background quantities are not evaluated at a particular Hubble patch like the ones mentioned above, but are in fact averaged over all Hubble patches (spatially independent) and only evolve with time. Their definitions are similar to those of Eqs. (2)-(3) but include a decay rateΓ φ , also averaged over all Hubble patches, as follows The Hubble parameter in this case is given by and the number of e-folds is Hdt .

B. Perturbed equations
In this Section we consider the evolution of the adiabatic matter and metric perturbations over the background quantities defined above.
Working in the Newtonian gauge, we perturb the metric of Eq. (1) in the absence of anisotropic pressure perturbations where we introduced the gravitational potential perturbation Φ. We also write the perturbations in the energy density of the inflaton and radiation components at a Hubble patch j as where δ j φ , δ j r 1 are small perturbations over the averaged background quantitiesρ φ andρ r respectively, at the particular patch. Despite the fact that, in principle, the perturbations themselves are dependent on the patch at which they are calculated, in the following we are going to refrain from using the j superscript in our Equations. That is because, as we will explain in Sec. II D, in our calculation we will only trace a characteristic value of the perturbations and not their actual distributions. This way we are later on able to constrain parameter space in a much more computationally efficient way. The actual perturbation distributions are calculated in our companion paper [37] with the purpose of deriving the corresponding non-gaussianity of the temperature anisotropy spectrum.
The gravitational potential perturbation Φ couples to the energy density perturbations through We now turn to the first order perturbations to the Boltzmann Eqs.
(2)-(3), allowing also for fluctuations δΓ φ in the inflaton decay rate. We follow the method used in Ref. [17] to assess the effects of such a perturbation. On super-horizon scales k aH, the expressions describing perturbations in the matter and radiation energy densities read where δ Γ ≡ δΓ φ /Γ φ is the perturbation in the inflaton decay rate. The background Eqs. (2)-(3) together with Eqs. (12)-(14) form a set of five coupled differential equations.
In the standard scenario where δ Γ = 0, super-horizon perturbations remain frozen until they re-enter the horizon (dΦ/dN = 0). In this work, where there is modulated reheating, on the other hand, δ Γ = 0 leading to dΦ/dN = 0; thus superhorizon perturbations evolve with time. In fact, the decay rate becomes time and spacedependent once the effects of Higgs modulation and Higgs blocking (as defined in Sec. I and discussed in FSSV) are taken into account during the reheating process. An alternative way of understanding the super-horizon evolution of the potential perturbation is to consider the isocurvature perturbations temporarily produced by the inhomogeneous transfer of energy between the inflaton and the radiation bath across different Hubble patches. We elaborate more on this topic in Sec. II F.

C. Higgs field dynamics
Having defined a framework for inhomogeneous reheating in previous sections, we now consider the dynamics of the Higgs boson during and after inflation, which determine the evolution of density perturbations produced during Higgs-modulated reheating. In this section, we summarize the aspects of the Higgs dynamics most relevant for Higgs-modulated reheating, while a more comprehensive discussion can be found in FSSV. We assume the SM Higgs is minimally coupled to gravity and is a light spectator field during inflation.
We take the background value of the Higgs doublet and its potential to be where ν = 246 GeV and λ is the quartic self-coupling, which is taken to be positive during inflation. Here h is a real scalar field. During inflation, the Higgs field initially rolls classically down its potential and soon reaches a regime dominated by quantum fluctuations. The result is that the super-horizon modes of the Higgs follow a random walk during the final stages of inflation. After a sufficient number of e-folds, the Probability Density Function (PDF) 3 describing the Higgs field at the end of inflation is [38] f eq (h) = 32π 2 λ I 3H 4 In the previous equation, λ I and H I are the Higgs quartic self-coupling and the Hubble scale at the end of inflation, respectively. Furthermore, the Gamma function has the value Γ(1/4) 3.625. Due to the stochastic dynamics mentioned above, each Hubble patch at the end of inflation has a different effective Higgs VEV. The probability to find a particular Higgs VEV in a given Hubble patch at the end of inflation is given by the equilibrium PDF of Eq. (17). The reheating dynamics within each Hubble patch (after inflation has ended) are completely deterministic, once we specify the initial Higgs VEV value sampled from Eq. (17). The dynamics of the Higgs field's VEV after inflation has ended in a Hubble patch j reads where the derivative of the Hubble rate with respect to the number of e-folds is For the moment, we have neglected the backreaction of the SM gauge bosons on the Higgs field dynamics, which is considered in Sec. III B. Eqs. (18)- (19) describe the damped oscillations which the Higgs experiences during the reheating period [33].
Using the distribution of Higgs values at any point in time, we can define a characteristic value of the Higgs. The mean value of the Higgs VEV h across the observable Universe is zero, since for every patch with a positive value h j > 0 there will be another patch where the Higgs VEV has the value −h j . We therefore take the characteristic value of the Higgs within a typical patch to be the standard deviatioñ Hereh is not the second moment of the equilibrium Higgs distribution in Eq. (17). Instead, it is a time dependent quantity and must be evaluated from the actual Higgs distribution at every point in time. Specifically we will computeh(N ) as a function of the number of e-folds N after inflation. Henceforth we will use the language "characteristic Higgs VEV" for this quantityh (although not actually a VEV itself, but rather a typical value within a given patch). At this point we should comment on the validity of Eq. (17) as the Higgs distribution at the end of inflation.f eq represents the equilibrium PDF of any light spectator field present during inflation with a quartic self-interaction term dominating its potential at large field values (in our case the Higgs). The equilibrium PDF has been derived under two main assumptions. First, the calculation assumes a de-Sitter spacetime during inflation even though the small scale variation of the observed CMB power spectrum suggests a slight deviation from the pure de-Sitter limit. We also assume there is a sufficient number of e-folds during inflation N equil ∼ O(λ −1/2 ) for the Higgs PDF to evolve towards equilibrium.
While Eq. (17) is sufficient to both demonstrate the important effects of Higgs dynamics on reheating and derive interesting constraints, a more thorough stochastic analysis of the Higgs dynamics could provide further insight and even better (possibly model-dependent) constraints. Recent analysis of stochastic dynamics of light spectator fields has uncovered interesting results; for example Ref. [39] showed that light spectator fields can acquire larger field displacements during inflation when accounting for deviations from the de-Sitter approximation. If applied to our calculations, the parameter space of reheating models could be more tightly constrained than what is shown in Figs. 9 and 10. Furthermore, the existence of four degrees of freedom in the Higgs field (ignoring gauge bosons for a moment), means that the Higgs random walk will be four-dimensional, leading to larger VEV's than its onedimensional counterpart. The existence of gauge fields complicate the actual calculation, but as shown in Ref. [40], the end result for the Higgs VEV is closer to that of a four-dimensional random walk than the onedimensional system leading to Eq. (17). Taking a conservative viewpoint, we use Eq. (17) as the basis of our calculations, keeping in mind that incorporating a more realistic PDF for the Higgs field will result in even tighter constraints.
One other concern we should address is the stability of the SM Higgs potential at large field values. The value of the Higgs quartic self-coupling λ ≡ λ(µ) depends on the renormalization parameter µ and can become negative at a high scale µ inst ∼ 10 11 GeV due to its renormalization group (RG) evolution, possibly leading to an instability in the potential. The random walk of the Higgs field during inflation could thus send the Higgs into an antide Sitter minimum [41][42][43][44]. This instability is sensitive to the value of the top quark mass m t which, for its best fit values leads to a negative Higgs potential above a (gauge-dependent) instability scale Λ inst ≈ 10 11 GeV. For simplicity, we assume that any possible instability in the Higgs potential is cured either by new physics decoupled from the inflation scale or by displacing the value of m t below its best fit value, within its significant experimental and theoretical uncertainties [45].

D. Perturbative decay: patch-by-patch Higgs evolution
In the previous sections we defined both the background and the perturbation equations which govern the growth of inhomogeneities in the case of an inflaton decay rate which varies between different Hubble patches. Furthermore, we showed how the VEV of the Higgs field can differ from patch to patch due to its random walk during inflation. We now move on to explain the dependence of the inflaton decay rate on the Higgs VEV, which causes variation in reheating between Hubble patches. In this subsection, we consider the case of perturbative decay of the inflaton. In Sec. II E we will turn to non-perturbative decay.
For simplicity, in computing the decay rate, we assume a Yukawa-type coupling between the inflaton and the fermion to which it decays. In this case, the decay rate of the inflaton at patch j was shown in our previous work in FSSV to be where m φ is the inflaton mass and Γ 0 is the inflaton decay rate in the absence of the Higgs modulation/blocking One can see that the standard deviation value of the Higgs PDF decreases with increasing N, while the decay rates approach the unblocked value Γ0 as blocking is lifted in more and more Hubble patches.
(the massless fermion limit). The Heaviside function, defined as Θ(x) = 1 for x ≥ 0 and Θ(x) = 0 for x < 0, accounts for the phase-space blocking due to large effective fermion masses. Although our calculations specifically assume a Yukawa-type coupling between the inflaton and SM fermions, the results we present in the case of perturbative inflaton decay should remain basically the same for any final state particles that become massive due to the SM Higgs. The details of the phase space of the decay products would vary for different inflaton-SM interactions or choice of final states but would not substantially change our results. The mass that the fermion acquires due to the Higgs mechanism in a patch with Higgs VEV h j is where y is the Yukawa coupling for a given SM fermion.
Although SM Yukawa couplings are technically scale dependent parameters, the RG evolution is typically insignificant between the electroweak and inflation scales for the minimal field content we consider in our analysis (i.e. the inflaton is the only new particle in addition to the SM). While we do consider the non-trivial RG evolution of other relevant SM parameters, in particular that of the Higgs quartic coupling λ (as discussed in Secs. II C and III), we assume the SM Yukawa couplings at the inflation scale are equivalent to the corresponding values at the electroweak scale. We thus present our results in terms of a generalized Yukawa coupling but make specific interpretations based on the hierarchy of Yukawa couplings relevant for the pattern of electroweak symmetry breaking in the SM. Fig. 1 presents the PDFs of the binned Higgs and inflaton decay rate values across all Hubble patches at N = 0, 1, 2.5 and 6.5 e-folds after the end of inflation (from top to bottom), for H I = m φ , y = 1, Γ 0 /m φ = 0.1 and λ I = 10 −3 . Note that the scales of the (left panel's) horizontal axes have been chosen to vary in going from the top to the bottom panels of the figure, to accommodate the decreasing width of the Higgs PDF as a function of increasing N 4 . The decrease in the width of the Higgs PDF is due to the decreasing amplitude of the damped Higgs oscillations; according to Eq. (18) the Higgs VEV becomes negligible within a few e-folds after the Higgs condensate begins to oscillate as a massive field.
At the same time, the damping of Higgs oscillations causes the decay rates Γ j φ in the right panel to approach their unblocked value of Γ 0 = Γ j φ (h j = 0). There are two effects that take place as the Higgs VEV h j in a Hubble patch decreases. First, the argument of the Heaviside function in Eq. (21) becomes positive and Higgs blocking is lifted. After the lifting of Higgs blocking the decay rates are within the range 0 < Γ j φ /Γ 0 < 1. Also, 4 Also see Fig. 7, where the blue solid line corresponds to δh = h 2 .
Higgs modulation due to the second factor of Eq. (21) (which we subsequently refer to as the phase space factor) approaches unity when the Higgs VEV has become much smaller than the mass of the inflaton. After both Higgs blocking and modulation are extinguished the decay rates are given by Γ j φ /Γ 0 = 1. Because Higgs modulation/blocking are lifted at different times in different Hubble patches, in Fig. 1 we can see some patches with 0 < Γ j φ /Γ 0 < 1, as well as others with Γ j φ /Γ 0 = 1 at N = 6.5 e-folds. In the following we will show that, even in the absence of Higgs blocking, Higgs modulation from the phase-space factor in the decay rate can lead to the production of large density perturbations.
Furthermore, it is worth noting that the Higgs PDFs in the left panel of Fig. 1 lose their original shape and become irregular after a sufficient amount of e-folds has passed. The reason is that Higgs oscillations proceed at different frequencies for different Higgs initial values |h j I |. However, since the oscillation frequency is the same for initial conditions with the same absolute value |h j I |, the Higgs PDFs remain symmetrical with respect to the vertical axis.
The "average" value of the decay rateΓ φ is obtained by using the characteristic value of the Higgs VEVh from Eq. (20) and plugging into Eqs. (21) and (22). We find In principle, we could also defineΓ φ in terms of h = 0, which would give the constant value Γ 0 . However, we are interested in calculating the density perturbations arising from the relative effects of Higgs modulation/blocking between different Hubble patches. Associated perturbations to the decay rate should thus be calculated with respect to the decay rate given by Eq. (23) rather than the unblocked decay rate Γ 0 , which would fail to capture any of the relevant modifications to the reheating dynamics.
The general expression for the perturbation of the decay rate relative to the averageΓ φ can be calculated as where δh is a variation of the Higgs VEV. We take δh = h 2 (corresponding to half of the distribution's width) and plug it into Eq. (24). The resulting characteristic perturbation is Our calculations to determine the decay rate perturbation δ Γ (N ) at every time-slice (labeled by the number of e-folds N ) after the end of inflation proceed as follows. We will treat a large number n ∼ 10 4 of Hubble patches, solving the equations independently for each patch. To obtain initial values for the Higgs field in multiple patches, we begin by drawing a sample of n values of the Higgs field from the Higgs PDF in Eq. (17), obtaining a collection of initial conditions h j I with j ∈ {1, ..., n}. Each value of the Higgs field corresponds to a different Hubble patch. Each value of h j I is used as the initial condition for solving the system of four coupled differential Eqs. (18)-(19), (2)-(3) describing the Higgs and energy density evolution respectively. This results in a collection of n solutions for the Higgs field h j = h j (N ), from which we calculate the value of h 2 (N ) that characterizes the Higgs distribution at every time slice.
To reiterate, in this section we introduced a method in which different Higgs VEVs in different Hubble patches evolve separately and are averaged at every time-slice in order for us to define the averaged background quantities and the corresponding perturbations. The benefit of this method is that it allows us to directly use Eqs. (21)- (25) in conjunction with Eqs. (13)- (14), to compute the density perturbations in matter (inflaton) and radiation during reheating. In other words, we are able to determine the time evolution of the inflaton decay rate and its 5 An alternative way of determining these two functions would be to separately calculate the time-evolution of the decay rates Γ j φ (j ∈ {1, ..., n}) in n Hubble patches. We could, then, define the average background decay rate asΓ φ (N ) = Γ φ and the decay rate perturbation as δ Γ = Γ 2 φ /Γ φ . A similar method is used in the case of resonant inflaton decay, as discussed in Sec. II E. In the present section we choose to make the dependence on the distribution of Higgs VEVs manifest in order to highlight the role of the Higgs field in the perturbation dynamics. We have, however, verified that the two methods give results which agree within O(10%). perturbation, both of which are necessary components of the perturbed Einstein equations introduced in Ref. [17]. Results of applying this method to the case of perturbative inflaton decay will be presented in Sec. III below.

E. Resonant decay: patch-by-patch energy densities
The approach used in Sec. II D to calculate the gravitational perturbation Φ is not applicable in cases where the inflaton decays non-perturbatively through parametric or tachyonic resonance. The reason is that in these cases we cannot define the decay rate Γ φ as we did in the case of perturbative decay. Several models of parametric resonance have been studied since early works on the subject [46,47]. Recent models include tachyonic preheating in α-attractors [48][49][50], Higgs inflation and related models [51][52][53][54][55], as well as the formation of structures (such as oscillons) during preheating and their gravitational wave (GW) signatures [56]. To demonstrate the effects of Higgs modulation/blocking on preheating, we will consider the example of a model exhibiting tachyonic resonance [57,58], specifically the case of an inflaton coupled to an abelian gauge field through a Chern-Simons term [59]. This paradigm is inspired by natural inflation [15], where the inflaton is an axion, possessing a shift-symmetry. Recent attention has focused on gauge field couplings to the inflaton potentially generating large scale magnetic fields [60] and a significant amount of GWs [61][62][63].
The effective Lagrangian for the system we consider is where φ is the inflaton and f is proportional to the breaking scale of the associated U (1) symmetry, expressed in units of the Planck mass m Pl . For the electromagnetic four-vector A µ , we have F µν = ∂ µ A ν − ∂ ν A µ and F µν = µνβγ F βγ with the totally antisymmetric tensor µνβγ . We did not explicitly introduce the Higgs field in the above Lagrangian, but its effects are included through the gauge field mass, which is determined by M = g|h|/2, where g is the gauge coupling. We use abelian gauge fields as a proxy for the full electroweak SM sector, as explained in detail in FSSV, and leave a full analysis of non-abelian effects (see e.g. Ref. [64]) for future work. Note that we refer to 1/f as the Chern-Simons coupling which is typically also proportional to the U (1) charge of the inflaton and the square of the gauge coupling.
For the calculation of the gravitational perturbations in the case of resonant decay we follow a slightly different method, focusing on the the transfer of energy from the inflaton to the gauge fields. We numerically solve the linearized equation of motion for the gauge field modes Figure 2. For the case of resonant inflaton decay, the energy density budget as a function of e-folds N after inflation for f = 0.1 m Pl and gauge field mass M/m φ = 0, 1 (blue and red respectively). We work in the linear fluctuation approximation, in which we neglect back-reaction effects. The solid curves correspond to the energy density in the produced gauge fields (radiation) ρr. The black-dotted curve shows how the energy density of the inflaton would evolve without any transfer of energy to radiation from resonant particle production, ρ bg . The colored dashed lines show our approximation of the true inflaton energy density ρ φ when we take into account the energy loss due to gauge field production. We approximate this as ρ φ ≡ ρ bg − ρr. The blue and red dots show our estimation for the time of complete preheating.
where the superscript ± denotes the two helicities. Here where χ k is the Fourier transform of the transverse component of the gauge field. As discussed in our previous work FSSV, we consider the point where the energy density of the linearized fluctuations equals the energy density of the unperturbed inflaton background to be indicative of complete preheating. Although this approximation does not account for the back-reaction on gauge boson production, lattice simulations have shown that tachyonic resonance of this form can efficiently preheat the Universe [60].
Since tachyonic resonance is strong enough for the parameters chosen to completely preheat the Universe within O(1) e-folds, we neglect the evolution of the Higgs condensate, taking it to be fixed at the value it has at the end of inflation in each Hubble patch. By starting with a distribution of Higgs values among different Hubble patches we define a similar distribution of gauge field masses M . By computing the energy density in radiation (gauge fields) and matter (inflaton condensate) in each patch (see Fig. 2), we define the averaged values ρ φ = ρ φ andρ r = ρ r over all patches and the corresponding fluctuations δρ φ = ρ 2 φ and δρ r = ρ 2 r at each time-slice. Finally, the inflaton and radiation perturbations are given by the definitions δ φ = δρ φ /ρ φ and δ r = δρ r /ρ r . Having calculated the energy density perturbation functions for the inflaton and radiation, we insert them into Eq. (12) to obtain the gravitational potential perturbation Φ. The results of this method are shown in Sec. IV.

F. Calculation of temperature anisotropies
In order to constrain the allowed parameter space for reheating, taking the effects of Higgs modulation/blocking into account, we must connect our results to the temperature inhomogeneities observed in the CMB. We expect those to directly depend on the gaugeinvariant comoving curvature perturbation, R. On superhorizon scales, R is equivalent to the Bardeen parameter ζ, defined as 6 [17,[65][66][67][68][69] Here Φ is the gravitational potential of Eq. (10), whilē ρ φ (ρ r ) and δ φ (δ r ) are the energy density background and perturbation on the background for the inflaton (radiation) respectively. In most models of inflation, the curvature perturbations are adiabatic and, thus, constant on superhorizon scales. However, in the case of modulated reheating, the situation is different. The curvature perturbations grow with time on superhorizon scales due to the spatial dependence of the inflaton decay. In order to gain intuition for this superhorizon growth, one can construct an unusual type of temporary isocurvature perturbation for superhorizon scales during reheating-the relative isocurvature perturbations between the inflaton and the radiation bath-as different amounts of energy are transferred from the inflaton to radiation in different Hubble patches, The associated time evolution of the Bardeen parameter is given by [70]ζ Note that these are not the usual isocurvature perturbations that lead to observables in the CMB (for example, the relative perturbations between radiation and matter that survive after the end of inflation). Instead, these are to be thought of as short-lived isocurvature perturbations that exist only during the reheating period. Furthermore, these isocurvature perturbations become rather ill-defined as the energy density of the inflaton vanishes in each Hubble patch. However, as mentioned above, it is informative to study their evolution during reheating to provide a better understanding for the way the Bardeen parameter grows on superhorizon scales.
In Fig. 3 we plot the time evolution of S (solid lines) andζ/H (dashed lines) calculated from respective Eqs. (29) and (30) for the case of perturbative inflaton decay. Both S andζ/H are shown as functions of the number of e-folds after the end of inflation (the beginning of the inflaton oscillations). One can see the rapid initial increase of S (near the characteristic time of perturbative inflaton decay described in the following section) followed by decay once energy is transferred to radiation. The time derivative of the Bardeen parameter re-scaled by the Hubble rate,ζ/H, rapidly decreases, approaching zero shortly after the end of inflation. As discussed further in Sec. III A, we have checked that taking the derivative of the Bardeen parameter in Eq. (28) gives the same result as the value ofζ shown in Fig. 3. From either perspective, we can see that we are left with only adiabatic perturbations once every Hubble patch has completely reheated after inflation. Thus, the isocurvature fluctuations play no role in the calculation of CMB temperature anisotropies. While for simplicity we have discussed the temporary isocurvature perturbations associated with Higgs-modulated reheating in the case of per-turbative inflaton decay, similar qualitative conclusions can be drawn in the case of reheating through resonant particle production.
The exact relation between the temperature fluctuations and the gravitational potential or Bardeen parameter must also account for dynamics taking place at later times and, in particular, during the decoupling of CMB photons from the primordial plasma. On scales that remain outside of the horizon at the time of last scattering, the geodesics of CMB photons are altered by the distortions of spacetime due to matter perturbations in what is know as the Sachs-Wolfe effect. Apart from the gravitational potential contributions, the full Sachs-Wolfe effect is calculated by taking into account perturbations intrinsic to the radiation plasma at the moment of photon decoupling. At linear order in the perturbations, the final result is expressed as [71,72] ∆T where Φ f and ζ f refer to the final values of the gravitational potential and Bardeen parameter respectively at the time of CMB decoupling. We are only interested in the largest scales observable in the CMB, which re-enter the horizon well into the matter dominated epoch. Thus, we approximate that the amplitude of superhorizon perturbations calculated through the end of reheating are conserved through last scattering.

III. RESULTS FOR THE CASE OF PERTURBATIVE INFLATON DECAY
In this section we discuss the temperature fluctuations produced by Higgs-modulated reheating in the case of perturbative inflaton decay, derived using the method described in Secs. II D and II F. Again, we assume a decay rate given by Eq. (21), corresponding to a Yukawa-type coupling between the inflaton and the fermion to which it decays. In Sec. III A we start with the simplest scenario in which we make the following two approximations: gauge bosons produced by the resonant decay of the Higgs condensate do not back-react on the Higgs evolution, as well as the assumption that the frequency of oscillations by Higgs the condensate is much slower than the Hubble rate. The effects of backreaction on the Higgs evolution and assuming the opposite limit of rapid Higgs oscillations will be investigated in Secs. III B and III C respectively. Finally, in Sec. III D we present constraints on a generalized parameter space for perturbative reheating by requiring that the amplitude of temperature fluctuations not exceed that observed in the CMB and then describe how our results could be extrapolated to models with a lower inflation scale in Sec III E.
Before proceeding, we will briefly comment on the value of Higgs self-coupling at the inflation scale λ I . Assuming no new physics couples to SM Higgs, RG evolution of λ between the electroweak and the inflation scale will yield λ I 10 −2 at the end of high scale inflation [73]. In the following sections we will primarily use λ I = 10 −3 since a smaller λ I causes the production of slightly larger adiabatic density perturbations and intensifies effects such as the backreaction of gauge bosons on the Higgs dynamics. After assuming λ I = 10 −3 in order to more clearly demonstrate various aspects of density perturbations produced by Higgs modulation/blocking, we will use λ I = 10 −2 as a benchmark when constraining the parameter space of reheating.
A. Simplest perturbative reheating case As mentioned above, we start with the simplest scenario in which we ignore the effect of gauge boson backreaction on the evolution of the Higgs condensate and assume the frequency of the Higgs oscillations is slower than the Hubble rate. In Fig. 4 we present the Bardeen parameters ζ as a function of the number of e-folds after the end of inflation. We have equated the Hubble scale at the end of inflation to the inflaton mass H I = m φ and used the value λ I = 10 −3 for the Higgs self-coupling.
As noted above, we see the growth of the Bardeen parameters coincide with that of the corresponding isocurvature perturbations in Fig. 3. More specifically, we can see the Bardeen parameter ζ of Fig. 4 stabilizes at the same time that its derivativeζ/H in Fig. 3 approaches zero. We have also checked that numerically differentiating the Bardeen parameters of Fig. 4 with respect to time yields the same values ofζ as those shown in Fig. 3. Thus, we conclude that the growth of adiabatic density perturbations on superhorizon scales can indeed be described as the temporary generation of isocurvature modes during the inhomogeneous reheating process.
We find two main results, which are clearly demonstrated in the simplest case of perturbative reheating but also generally hold under more complicated assumptions discussed in subsequent sections. First, we find that the perturbations caused by Higgs modulation/blocking, ∆T /T | H ∼ ζ/5, are larger for larger values of the Yukawa coupling y and the decay rate Γ 0 . In fact, they can exceed the temperature fluctuations observed in the CMB ∆T /T | CMB ∼ 10 −5 by several orders of magnitude for certain parameter combinations. For example, we find ). An extensive examination of the full parameter space (Γ 0 , y, λ I ) and bounds from the CMB will be presented in Sec. III D.
Second, we see that Higgs blocking of the inflaton decay into fermions takes place only for y 1. In Fig. 4, we observe that the Bardeen parameter in the case of y = 1 grows sharply only after N 3 e-folds, whereas for smaller values of y the growth happens much sooner and more gradually. This trend is consistent with the results of FSSV, where we showed that large Yukawa couplings are needed to cause a significant delay of the reheating process. In that paper we showed that the reheat temperature could be suppressed by up to an order of magnitude compared to the unblocked case.
In Fig. 4 one can see that larger values of the inflaton decay rate lead to more rapid increase of the Bardeen parameter. This trend can be understood by comparing the timescale relevant for Higgs blocking to that of reheating, which can only be completed when Γ φ H. For example, at large values of y = 1 and Γ 0 /m φ = 0.1, the decay rate of the inflaton is already larger than the Hubble rate by the time Higgs blocking has been lifted. Hence, after blocking is no longer an obstacle, reheating takes place instantaneously and ζ increases sharply. On the other hand, for y = 1 and smaller decay rates (right panel, with Γ 0 /m φ = 10 −3 , 10 −4 ), the expansion rate is still larger than the decay rate when Higgs blocking is lifted. As a result, reheating requires more time and the increase of ζ happens more gradually.
For cases where y 0.1, there is never any Higgs blocking; yet substantial density perturbations may still result simply due to Higgs modulation. In other words, the Yukawa couplings are not sufficiently high to make the argument of the Heaviside function in Eq. (21) negative, and thus do not block reheating (or block reheating only at an exponentially suppressed number of rare Hubble patches). However, whether or not Higgs blocking is ever significant, there is always a residual dependence of the phase space factor on h, which creates differences between decay rates at different Hubble patches. This Higgs modulation can thus still lead to the production of large density perturbations.
It might be possible that, due to the deviation from a pure de-Sitter spacetime during inflation, the equilibrium distribution of Higgs VEVs given by Eq. (17) is not a valid approximation and larger field values could be expected from the quantum fluctuations of the spectator Higgs field [39]. For larger Higgs VEVs across a significant number of Hubble patches, smaller Yukawa couplings would be sufficient to produce the same results calculated under the assumption of a de-Sitter equilibrium distribution. More specifically, the curves shown in Fig. 4 for a given combination of (Γ 0 , y, λ I ) can be approximately re-interpreted as corresponding to the evolution of the Bardeen parameter with (Γ 0 , y , λ I ), where the re-scaled Yukawa coupling is given by y = (h/h )y for the modified characteristic Higgs VEVh . While a detailed analysis of the changes to our results for deviations of the Higgs PDF from Eq. (17) are beyond the scope of this work, the possibility of larger Higgs VEVs across a sufficient number Hubble patches would imply that the constraints on the parameter space of reheating in this paper should be considered conservative.

B. Including the effects of backreaction on Perturbative Reheating
The results discussed in the previous section were obtained by neglecting the effects of the gauge bosons which are produced resonantly from the decay of the Higgs condensate [29]. These gauge bosons can potentially back-react on the dynamics of the Higgs boson. In this section, we summarize the treatment of backreaction discussed in FSSV and show that the associated effects on the density fluctuations produced during Higgs-modulated reheating may be neglected.
Gauge Boson Production: The induced mass of the SM W-bosons m W = g|h|/2 vanishes when the oscillating Higgs field crosses zero, leading to a substantial production of gauge boson particles [38,73]. During the oscillations of the Higgs field, the transverse components of the SM gauge fields oscillate with a time-dependent frequency ω k depending on the mode with wavenumber k, and with a corresponding occupation number n k [29,47]. The backreaction associated with the resonant production of the gauge bosons on the Higgs dynamics manifests as an effective mass squared term in the Higgs equation of motion, Eq. (18), given by where g is the coupling constant appearing in the covariant derivative of the Higgs to the gauge field. The resonant production of the gauge bosons is governed by the quantity q W ≡ g 2 /(4λ), see Appendix A for further discussion. We ignore the non-Abelian self-interactions of the gauge fields [74], which may change the Higgs condensate decay time but would not drastically affect our overall results. Backreaction: The backreaction from gauge bosons takes effect when the effective mass squared of the gauge bosons m 2 h(W ) is of the same order as the effective mass squared m 2 h(λ) of the Higgs field given by its selfcoupling [29], When the condition in Eq. (33) is met, we assume the Higgs field instantaneously decays away and the inflaton decay rate-no longer blocked-becomes equal to Γ 0 . The full decay rate of the inflaton, accounting for the effects of backreaction (BR), is therefore [33] The expression for Γ BR above imposes that Higgs blocking only affects the dynamics during the time period when backreaction can be neglected, m 2 h(W ) ≤ 3λ I h 2 . Following the analysis in FSSV, the time from the end of inflation when the Higgs condensate decays away due to backreaction effects reads where the factor µ 0 = 0.185 is obtained from a numerical fit to the exact solution for the time dependence of n k [33]  Relative difference (%) y = 0.1 y = 0.5 y = 1 Figure 5. For the case of perturbative inflaton decay, the fractional difference between perturbations, (ζ − ζBR)/ζ, with (ζBR) and without (ζ) including the effects of backreaction, as a function of the unblocked decay rate Γ0, for three different Yukawa couplings y = 0.1 (blue dotted), y = 0.5 (red dashed), and y = 1 (green solid). We assume HI = m φ and λI = 10 −3 . and the occupation number at wave number k = 0 is In the following we point out the conditions for which backreaction can impact the calculation of the density perturbations produced by Higgs-modulated reheating. First, the unblocked decay rate of the inflaton field Γ 0 must be low enough to ensure that the Higgs condensate will decay before reheating has been completed. In other words, since the "seed" of the perturbations from Higgsmodulated reheating is the difference in decay rates between different Hubble patches, only further modifications to the decay rates will cause the calculation of perturbations to change. The effects of backreaction can only modify the decay rate in a given Hubble patch if the Higgs field is given sufficient time to decay before reheating has been completed and the decay rate has become the same (equal to Γ 0 ) in every patch. More specifically, the relation between the unblocked inflaton decay rate Γ 0 and the Hubble parameter at the time of the Higgs condensate decay H dec , should be Γ 0 H dec .
The second condition for the effects of backreaction to be important requires that the Yukawa coupling y is high enough, such that the decay rate of the inflaton does not become immediately equal to Γ 0 after inflation has ended. As mentioned earlier, if Higgs blocking never occurs the only source of the perturbations is the phase space factor causing slight differences in the decay rates of different patches. The smaller the Yukawa coupling is, the smaller the deviations from Γ 0 caused by the phase space factors are. As a result, for Yukawa couplings small enough such that the transition fromΓ φ to Γ 0 is immediate once Higgs blocking is lifted, the effects of backreaction do not cause considerable modification of the decay rates. Therefore, the density perturbations calculated for small y when accounting for the effects of backreaction, are not considerably different compared to those calculated in Sec. III A ignoring the effects of backreaction.
We demonstrate the conditions under which the effects of backreaction are relevant in Fig. 5, where we plot the fractional difference between perturbations with and without backreaction (ζ − ζ BR )/ζ as a function of the unblocked decay rate Γ 0 , for the values of the Yukawa coupling y = 0.1 (blue dotted), y = 0.5 (red dashed), and y = 1 (green solid), while setting λ I = 10 −3 . It is evident that changes to the size of density perturbations due to the effects of backreaction increase for larger Yukawa couplings and decrease for larger decay rates Γ 0 , thus verifying our intuitive understanding described above. More specifically, deviations between the two methods approach zero for Γ 0 /m φ 10 −2 . Furthermore, the largest differences between the perturbations with and without backreaction, in the case of Yukawa couplings y 1, are well below an order of magnitude and require extremely small decay rates to become significant. Fig. 6 compares the evolution of density perturbations with and without including the effects of backreaction. In order for the effects to be visible by eye, we only show examples with relatively low decay rates Γ 0 = 10 −3 m φ and Γ 0 = 10 −4 m φ . We plot the ratio of the Bardeen parameters calculated when including the effects of backreaction to those in the right panel of Fig. 4 where backreaction was ignored. As expected from Fig. 5, the largest modification of the perturbations occurs for the combination of the largest Yukawa coupling y = 1 with the smallest decay rate Γ 0 = 10 −4 m φ . Since Yukawa couplings as large as y = 1 are only relevant for the case of the inflaton decaying into the SM through a coupling to the top quark and the associated perturbations are much larger than observed in that part of our reheating parameter space, we choose to ignore effects of backreaction in the following sections. 7 The conclusions about the effects of backreaction drawn so far have been based on calculations assuming a rather small value of the Higgs quartic coupling λ I = 10 −3 . However, the approximation we make by ignoring backreaction holds for larger self-couplings because the second term of Eq. (35) contains a contribution proportional to the logarithm ln(λ I ) [33]. Thus, larger couplings λ I > 10 −3 result in longer decay times for the Higgs condensate, t dec . The effects of backreaction for larger quartic couplings are, according to the first condition described above, only relevant for even smaller values of Γ 0 in order for the Higgs field to decay before reheating is completed. Choices of λ I > 10 −3 further limit the portion of parameter space where the effects of backreaction are relevant and we are therefore justified in ignoring the effects of backreaction for cases with larger quartic couplings subsequently considered in this paper.

C. The case of rapid Higgs oscillations in Perturbative Reheating
After inflation has ended, the Higgs experiences damped oscillations. For simplicity, we calculate the Higgs evolution in each Hubble patch by considering two limiting cases where the oscillation period of the Higgs field τ Higgs is either much longer or much shorter than the Hubble time, H −1 . In the work presented so far, we have assumed the limit τ Higgs H −1 , which allows for the resolution of the zero-crossings in the oscillations of the Higgs that gradually cause Higgs modulation/blocking to be lifted [33]. Under such an assumption, it is sufficient to directly sample values of the Higgs field from the distribution of Higgs VEVs computed at every time-slice, as described in Sec. II D.
If, however, many Higgs oscillations occur during one Hubble timescale τ Higgs H −1 , the oscillations cannot be accurately resolved. We thus define an effective value of the Higgs field as 7 As we will see below there is a small portion of allowed parameter space with Yukawa couplings equal to 1 for extremely small decay rates. However, considering the small size of this region in parameter space and the insignificance of the backreaction effect in any case, ignoring it remains a very good approximation.  37) and is given by δh rapid = h 2 eff . In the case of the latter, the width is calculated using the PDF of Higgs VEVs h presented in Fig. 1, as δh slow = h 2 . We use HI = m φ , y = 1, Γ0/m φ = 0.1 and λI = 10 −3 . in a patch denoted by j, where ρ h j is the energy density in the Higgs condensate at the particular patch. A comparison between the time-evolution of the Higgs PDF width δh under the assumption of rapid oscillations (τ Higgs H −1 ) and our standard scenario (where τ Higgs H −1 ) is shown in Fig. 7. For rapid oscillations, the width corresponds to the standard deviation of a PDF of effective Higgs values h eff , as defined in Eq. (37), and is given by In our standard scenario, the width is calculated using the PDF of Higgs VEVs h presented in Fig. 1, as δh slow = h 2 . Since h = h eff = 0 for a symmetrical Higgs distribution, this figure also describes the evolution of the characteristic Higgs values h = h +δh slow andh eff = h eff +δh rapid , which governs bothΓ φ and δ Γ . We can see that in the rapid oscillation case δh rapid and henceh eff decrease more slowly. As a result, Higgs blocking is lifted somewhat later under the assumption of rapid Higgs oscillations when compared to the scenario studied in Sec. II D with τ Higgs H −1 . More typically, the Higgs oscillations would take place in the intermediate regime τ Higgs H −1 , so that the actual perturbation values lie between the results derived in the two limiting scenarios mentioned above. In Fig. 8 we present a comparison between results derived in the two regimes for Yukawa couplings y = 0.01, y = 0.1 and y = 1. Similarly to the backreaction comparison, we plot the ratio of Bardeen parameters calculated in the rapid Higgs oscillation scenario (τ Higgs H −1 ) to those shown in Fig. 4 for the slowly-oscillating Higgs. For the case of y = 1 (green lines), the initial value of the ratio begins at 0 instead of 1 due to the effect of Higgs blocking 8 which is effective for a longer time when assuming rapid Higgs oscillations. Also when assuming τ Higgs H −1 , sampling h eff instead of h causes Higgs blocking to be lifted slightly later, as explained above. As a result, ζ becomes non-zero, while ζ rapid is still blocked, thus leading to the ratio ζ rapid /ζ smoothly increasing from zero.
We also observe in Fig. 8 that the importance of how we treat the Higgs oscillations is largely dependent on the decay rate of the inflaton. More specifically, smaller decay rates amplify the differences between the perturbations calculated in standard and rapid oscillation scenarios because reheating takes place later. As a result, for lower Γ 0 the Bardeen parameters increase the most between N = 2 and N = 5 e-folds, during which the widths δh differ the most between the two cases (cf. Fig 7). If, on the other hand, the decay rate is larger then the Bardeen parameters increase most during the first 2 efolds, when δh is the same in both cases. Therefore, the differences between the associated perturbations are smaller for larger Γ 0 .
More generally, we find that the Bardeen parameters calculated assuming rapid Higgs oscillations (τ Higgs H −1 ) do not differ from the standard case (τ Higgs H −1 ) by more than O(1) factors; i.e. ζ rapid /ζ ∼ 1 in Fig. 8. Since, as mentioned above, the Bardeen parameters calculated under the more realistic assumption of τ Higgs H −1 should lie in between those calculated in the two limiting cases, reasonably good agreement between both calculations allows us to confidently ascertain the corresponding constraints from temperature anisotropies in the CMB (see Sec. III D).

D. Comparison of Temperature Fluctuations from space-dependent Perturbative Reheating with CMB data
In this section we compare results of our calculations of the Bardeen parameter at the end of reheating with the value of the temperature anisotropy measured in the CMB, ∆T /T CMB ≈ 10 −5 . Since the predictions of the temperature fluctuations in our scenario depend on several parameters, we may use the CMB bounds to place constraints on the associated parameter space shown in Fig. 9. Specifically we examine the dependence of the Figure 9. For the case of perturbative (space-dependent) inflaton decay to SM particles, constraints on parameters obtained by requiring that temperature fluctuations do not exceed CMB observations. Here Γ0/m φ is the unblocked inflaton decay rate in units of inflaton mass; y is the Yukawa coupling of the SM particles to the Higgs; and λI is the Higgs self-coupling at the scale of inflation. In the top panels of the figure, we depict the three-dimensional parameter space (λI , y, Γ0), and in the bottom panels we depict two-dimensional slices at the value λI = 10 −2 , which is typical in the SM for inflation at high scales. The left (right) panels show results derived in the τHiggs H −1 (τHiggs H −1 ) regime of slow (rapid) Higgs oscillations. The red and green regions show Higgs-induced temperature inhomogeneities ∆T /T |H which are larger and smaller than those observed in the CMB, respectively. Hence the white and green regions are allowed as they satisfy ∆T /T |H ≤ 10 −5 . The hatched region at ∆T /T |H ≤ 10 −7 indicates the untested regime where the temperature fluctuations ∆T /T |H are smaller than Planck's O(1)% sensitivity (e.g. if we assume that the origin of the temperature anisotropies observed by Planck arise from the quantum fluctuations of the inflaton itself). For the intermediate regime of Higgs oscillations with τHiggs H −1 , the amplitude of temperature fluctuations would lie somewhere in between the values shown for the slow-(left panels) and rapid-oscillation (right panels) approximations. temperature fluctuations on the Higgs self-coupling λ I , the Yukawa coupling y of the fermions to the Higgs, and the unblocked decay rate of the inflaton Γ 0 . In the top panels of the figure, we depict the three-dimensional parameter space (λ I , y, Γ 0 ), and in the bottom panels we depict two-dimensional slices at the value λ I = 10 −2 which is typical for the SM in high scale inflation. We explore a wide range of different values for the unblocked decay rate 10 −5 ≤ Γ 0 /m φ ≤ 10 −1 and the Yukawa coupling 10 −3 ≤ y ≤ 1.
In all diagrams, red regions correspond to parameter combinations that lead to Higgs-modulated fluctuations in excess of what is observed in the CMB (∆T /T | H ∆T /T CMB ); and green regions correspond to Higgsinduced fluctuations that are below the CMB values (∆T /T | H ∆T /T CMB ). The white region indicates Higgs-induced fluctuations of the same size as those of CMB observations ∆T /T H ≈ 10 −5 . Hence both white and green regions are in principle allowed by the CMB, while the red region is certainly excluded.
We assume CMB observations are not sensitive to Higgs-induced perturbations with ∆T /T | H 0.01 ∆T /T CMB ∼ 10 −7 . The reason is that such small temperature anisotropies are below the O(1%) sensitivity of Planck. Then the observed anisotropies of the CMB ∆T /T CMB ≈ 10 −5 must be produced via some other mechanism, e.g. the standard density fluctuations arising from quantum fluctuations of the inflaton field. The regime for the Higgs-induced fluctuations ∆T /T | H 10 −7 is indicated by the hatched region of Fig. 9.
The left (right) panels of Fig. 9 are for the cases of slow (rapid) Higgs oscillations. Sampling the effective Higgs values of Eq. (37) (rapid oscillations, right panels of Fig. 9) leads to slightly larger overall perturbations than the case of slow oscillations. However, one can see by comparing the left and right panels that the observed differences are very small and much below an order of magnitude. As explained in Sec. III C, the true period of the Higgs oscillations could in fact be in the intermediate regime between the two limiting cases we investigate, with τ Higgs H −1 . In the intermediate case, the amplitude of temperature fluctuations would lie in between those calculated in the slow-and rapid-oscillation approximations, which are in reasonably close agreement.
Let us now summarize the basic parameter dependencies of Fig. 9. The Higgs quartic self-coupling λ I at the end of inflation only appears in the three-dimensional plots in the top panels of Fig. 9. It is evident that the size of temperature fluctuations decreases as λ I increases from 10 −3 to 10 −1 . This effect can be explained by the width of the Higgs PDF at the end of inflation, given by Eq. (17), decreasing for larger λ I . In other words, Higgs PDFs with a larger value of λ I are peaked around smaller values of the Higgs field h and give smaller probabilities for larger values of h to exist in a given Hubble patch. Since the Higgs modulation that causes the temperature fluctuations is more pronounced for larger values of the Higgs field h, a larger λ I leads to smaller overall ∆T /T | H .
The unblocked decay rate Γ 0 of the inflaton to fermions and of the Yukawa coupling y of the same fermions to the Higgs appear in both the three-and the two-dimensional plots of Fig. 9. We can see larger overall temperature fluctuations arise for larger y and for larger Γ 0 . Larger y leads to larger fermion masses and, as a result, to more significant Higgs modulation of the decay rate. Since Higgs modulation is the seed of the temperature fluctuations we examine, larger y will also lead to larger overall temperature fluctuations. The effects of increasing the decay rate are more subtle, but can be understood from Fig. 4. The inflaton decays at earlier times in each Hubble patch for larger Γ 0 . Values of the Higgs VEV in each patch, evolving according to Eq. 19, are larger at earlier times. Thus, larger Γ 0 leads to larger fermion masses when the inflaton decays, also causing larger overall temperature fluctuations.
Here we note several caveats in our calculations due to some of the other approximations we have made. First, when deriving the perturbations we only treat the largest observable scales in the CMB. This simplification is inevitable since our method uses the one-point PDF of the Higgs field, which lacks any information regarding scaledependence. Therefore, our results are consistent with the assumption of a scale-invariant power-spectrum of Higgs fluctuations and remain a very good approximation (up to factors of ∼ O(1)) for mild scale dependencies.
In a related approximation mentioned in previous sections, our calculations assume a purely de-Sitter spacetime when the largest observable scales exit the horizon during inflation and that the PDF of the superhorizon modes of the spectator Higgs field has reached its equilibrium. Without this assumption we would not be able to sample the equilibrium distribution function of Eq. (17) for the initial condition of the Higgs in each Hubble patch. We leave a more detailed treatment of the Higgs fluctuations for future work, noting in particular that calculations valid at smaller angular scales could allow for constraints in the (unhatched) green region of Fig. 9.
In order to more clearly demonstrate the parameter dependence of the temperature fluctuations ∆T /T | H produced by Higgs modulated reheating, we estimate a fitting function in terms of the Yukawa coupling y and the unblocked decay rate Γ 0 based on parameter space data shown in the bottom left panel of Fig. 9, Using this fitting function, we can constrain parame-ter values 9 for each SM particle separately based on its Yukawa coupling, provided that the inflaton decay channel we are examining is the dominant one during reheating. As a reminder, the inflaton can always decay into massless photons without the Higgs influencing the size of the produced temperature fluctuations. Furthermore, while there is a running of the fermionic Yukawa couplings between the Electroweak and the inflation scale, we do not expect them to be significantly modified and, thus, we choose to use their standard EW values. We compare the Higgs-induced temperature fluctuations (calculated from Eq. (38)) produced for different Yukawa couplings with the CMB value of ∆T /T | CMB ≈ 10 −5 and place upper bounds on the quantity Γ 0 /m φ where Γ 0 is the unblocked decay rate. For the cases with Yukawa coupling y < 1, we previously showed in FSSV that Higgs blocking is minimal and we can approximate the time of reheating by Γ 0 ∼ H. The reheat temperature T reh is then given directly in terms of Γ 0 as where g * ≈ 106.75 is the number of relativistic degrees of freedom of the radiation bath. The same relation applies even for cases with y = 1 where Higgs blocking is present, provided that the unblocked decay rate Γ 0 is sufficiently low. This condition ensures that blocking will be lifted before Γ 0 ∼ H. When this condition is not fulfilled, Eq. (39) gives a reheat temperature larger than its actual value since it does not take the delays due to Higgs blocking into consideration. In Sec. III A we explained how Fig. 4 shows that for y = 1 and Γ 0 /m φ 10 −2 Higgs blocking is lifted before Γ 0 ∼ H. In the following we will see that our upper bound for Γ 0 in the case of the topquark (y t = 1) is Γ 0 /m φ | max top 10 −2 and, thus, Eq. (39) still applies. 10 Hence, by setting upper bounds on the decay rate Γ 0 we can immediately constrain the reheat temperature for inflaton decays to all SM fermions.
We assume the maximum allowed value of the unblocked decay rate to be Γ 0 10 −1 × m φ in order to ensure that the decay of the inflaton does not arise from a strongly coupled theory. A corresponding lower bound does not exist, provided that reheating is complete by the time of Electroweak Symmetry breaking (EWSB) or, at the latest, Big Bang Nucleosynthesis (BBN), whose energy scale is much lower than the inflationary one [75]. According to Eq. (39), the temperature of EWSB T EWSB ≈ 160 GeV corresponds to 9 The values mentioned are subjected to the error of our fit. For more accuracy see Figs. 4 and 9. 10 Sec. III A uses λ I = 10 −3 , whereas here we use λ I = 10 −2 . However, for larger values of the self-coupling λ I > 10 −3 blocking is lifted even sooner and, thus, our statement remains accurate. Additionally we have made sure that blocking is lifted before Γ 0 | EWSB ≈ 10 −24 × m φ , which is much below any of the limits we are setting. Stronger lower bounds on T reh would apply if we wanted (vanilla) thermal leptogenesis to provide the origin of today's matter/antimatter asymmetry; then the bounds would be T reh > 10 9−11 GeV depending on the model (for example, see Ref. [76]). Let us consider the potentially most constrained case, i.e., where the inflaton primarily decays to the top quark, the most massive of the fermions with the largest SM Yukawa coupling y t ≈ 1. We find that the upper bound on the decay rate is Γ 0 /m φ | top < 3 × 10 −6 corresponding to a bound on the reheat temperature This bound may be in conflict with models of thermal leptogenesis, depending on the mass of the inflaton and the lightest right-handed neutrino, but will never be in conflict with EWSB or BBN.
For the electron and the muon, with y e ≈ 3 × 10 −6 and y µ ≈ 6 × 10 −4 respectively, any value of the unblocked decay rate yields temperature fluctuations which are unconstrained by CMB observations and lie within the hatched region of Fig. 9 (bottom panels). For inflaton decays to the the tau lepton, with y τ ≈ 10 −2 , we find the constraint Γ 0 /m φ |τ < 0.04 corresponding to

E. Lowering the scale of inflation
In the previous sections we have explored the temperature fluctuations arising from Higgs modulation/blocking of reheating after high-scale inflation, by setting the Hubble parameter at the end of the inflationary epoch equal to the mass of the inflaton H I = m φ . However, it is also useful to look into cases where the Hubble parameter H I takes much smaller values and, thus, investigate what our perturbations would look like if inflation were to occur at a lower scale.
Upon examination of the equations which govern the production of the gravitational and Bardeen parameters, as well as numerically solving the same system of equations for various different combinations of the input parameters, we see that inflation at a lower scale can produce the same perturbations as those of high-scale inflation. Re-scaling of our results for the amplitude of temperature fluctuations is possible provided that the free parameters are chosen such that where the subscripts "hs" and "ls" refer to high and lowerscale inflation respectively. Another necessary condition for the applicability of our results to inflation at lower  scales is that the unblocked decay rate Γ 0 should be smaller than the value of the Hubble parameter H I by the same amount as for inflation at high scales, Inserting Eq. (43) into Eq. (42) gives and, therefore, Eqs. (43) and (44) are the two independent conditions for high-scale and lower-scale inflation to result in the same perturbations.
The re-scaling of the results of our paper as a function of inflation scale can be understood by examining the role of these two conditions. Going back to Eqs. (2)-(3) for the background evolution, we see that the factor setting different Hubble patches apart is We notice that the two important quantities Γ 0 /H I and y h/m φ ∝ y H I /m φ are the same regardless of the inflationary scale if the parameters are chosen such that they match the conditions in Eqs. (43)- (44). In Eqs. (12)- (13), which govern the growth of perturbations in our model, the term of importance which essentially fuels the potential production, is We, therefore, conclude that a lower inflation scale H I will lead to a re-scaling of the parameter space introduced in Fig. 9 and allow larger values of the Yukawa couplings y. In other words, the green region of the same Figure will move further up to larger values of log(y) and our constraints will be relaxed. Table I presents the largest Yukawa couplings allowed by observations for certain choices of the unblocked decay rate Γ 0 and of the inflation scale H I . We can see that as the scale of inflation goes down the allowed values of y increase, provided that Γ 0 is chosen such that Eq. (43) is satisfied. For example, for inflation scale H I = 10 −4 m φ (four orders of magnitude lower than our canonical case), and taking decay rate Γ 0 /m φ = 10 −5 to match the condition in Eq. (43), the Yukawa coupling y can be as large as 100 without violating CMB constraints. 11 In order to derive these results we have taken m φ to be the same both in high-scale and in the case of a lower inflationary scale, even though it could in principle vary between the two.
A commonly used example in the literature for constructing models with a low inflationary scale consistent with the CMB are α-attractor models, like the Tmodel [77], whose single-field potential is The mass scale µ = O 10 −6 M Pl is chosen to set the amplitude of curvature perturbations ζ ∼ 10 −5 , in line with the results of the Planck collaboration [10]. The parameter α controls the tensor-to-scalar ratio and the scale of inflation, which scales as H ∝ √ α. This potential shows a simple realization of a model where lowering the inflationary scale does not affect the inflaton mass close to the minimum of the potential, m 2 φ = µ 2 /3. It has nevertheless been shown that the T-model at low α preheats efficiently both through self-resonance and through parametric resonance of a companion spectator field [48,50,78].
In general, keeping a large inflaton mass while lowering the energy scale of inflation will likely introduce significant non-linear terms to the potential that will lead to significant self-resonance. This is required, since low-scale single-field slow-roll inflation requires a flat "plateau" in the potential and the transition from a large inflaton mass at low field values to a flat plateau at large field values requires non-linear terms in the potential. Keeping this in mind, a more thorough investigation of Higgs blocking effects in low-scale inflation- 11 While values of the Yukawa coupling y 4π are larger than what is typically permitted by perturbative unitarity, as discussed in FSSV the relevant parameter is actually the mass of the fermion given by Eq. (22) and, thus, such large Yukawa couplings can be considered representative of cases with y 4π and proportionally larger Higgs field values (cf. discussion at the end of Sec. III A).
ary models must take self-resonance into account. Assuming that self-resonance occurs efficiently, the Universe at the end of low-scale inflation will be populated by scalar (inflaton) particles whose mass will be much larger than the Hubble scale, at least for the α-attractor potential of Eq. (47). The momentum of these particles will be of the order of their mass, thus they will be either non-relativistic or slightly relativistic. We leave a detailed (and somewhat model-dependent) analysis of Higgs-modulated reheating effects in low-scale inflation for future work.

IV. RESULTS FOR THE CASE OF RESONANT INFLATON DECAY
We now present the effects of Higgs blocking and modulation on the generation of perturbations during gauge preheating, as a representative case of resonant decay of the inflaton. In FSSV we explored the effects of Higgs blocking on gauge preheating and found that complete preheating is only possible when the mass of the gauge bosons is M H I , with the exact threshold depending on the Chern-Simons coupling strength. The upper panel of Fig. 10 shows the probability distribution of gauge field masses M = g|h|/2 for two values of the gauge coupling g = 0.8, g = 0.1 and two values of the Higgs self-coupling λ = 10 −2 , 10 −3 , based on the probability distribution of Higgs values, given in Eq. (17).
For the case of λ = 10 −3 , there is a wider range of gauge field masses in different patches, reaching up to M 2H I for g = 0.8. Those Hubble patches in which the gauge field mass satisfies M H I will successfully preheat resonantly; on the other hand those Hubble patches in which the gauge field masses are larger (M > H I ) will not have successful resonant preheating and must reheat later through e.g. perturbative decay of the inflaton through the same coupling to gauge bosons. In the case of efficient preheating, the reheat temperature will be T reh ∼ m Pl H I , while for perturbative decays the reheat temperature can be lowered by orders of magnitude [33]. This range of behaviors from one Hubble patch to another will lead to large temperature fluctuations ∆T /T ∼ 1 for the case of λ = 10 −3 .
We thus restrict our analysis to λ = 10 −2 , where the gauge field masses predominantly satisfy M H I (as seen in the upper panel of Fig. 10), and assume the Chern-Simons coupling is large enough such that the entirety of the observable Universe can completely preheat through parametric resonance of gauge bosons. Although these assumptions prevent the generation of manifestly large temperature fluctuations associated with incomplete preheating in a significant number of Hubble patches, each patch will nevertheless preheat at a slightly different time. We thus explore how the distribution of gauge boson masses leads to inhomogeneous preheating and determine if the amplitudes of the associated perturbations are in agreement with CMB observations. Following the discussion of Sec. II E we solve Eq. (27), the linearized evolution equation of the gauge field modes, for a grid of comoving wave numbers k, starting enough e-folds before the end of inflation, so that we can reliably initialize each mode in its Bunch-Davies vacuum state. The middle and lower panels of Fig. 10 show the results of using the linearized analysis of fluctuations until the point when the radiation energy density equals the inflaton background energy density in each patch. In order to extend the energy density evolution beyond the point of complete preheating, we assume that the entire energy density of each patch subsequently red-shifts like radiation, ρ ∝ a −4 . This introduces a sharp "knee" around N 1.2 e-folds, which is evident in the evolution of the Bardeen parameter ζ.
In order to validate the approximation which assumes an instantaneous transfer of energy into radiation, we repeat the calculation of the Bardeen parameter but numerically smooth the transition from the inflaton dominated epoch to radiation domination. The evolution of the Bardeen parameters are shown using both approximations in the middle and lower panels of Fig. 10. It is clear that the final amplitudes of the perturbations are largely independent of how we treat the transition between matter and radiation domination. The bulk of the evolution of the perturbations occurs within the first e-fold after inflation, when the gauge fields are amplified enough to account for a non-negligible part of the energy density of the Universe, at the percent level or above. Smoothing the transition between epochs does indeed make the evolution of ζ smoother, but altering the precise details of the energy densities will not significantly change the amplitudes of the perturbations. Furthermore, the middle and lower panels of Fig. 10 show that, in the case of parametric resonance, the effects of Higgs modulation lead to perturbations that are significantly larger than those observed in the CMB for both g = 0.8 and g = 0.1. Despite the uncertainty in the exact running of the SM couplings to the inflation scale, gauge couplings within the range 0.1 ≤ g ≤ 0.8 can represent a wide variety of different models. As a result, preheating into Higgsed gauge bosons cannot be the main source of reheating the Universe, at least given the fairly generic assumptions described in Sec. II. If the inflaton couples to both a massless gauge boson (photon) and to massive ones (W ± and Z), the problem becomes highly parameter dependent, as it depends crucially on the strength of the Chern-Simons coupling to each gauge boson. We expect the effect to still be present in the full electroweak sector, but suppressed compared to our current analysis. A key factor in this process will be the ratio of the energy density of the inflaton that ends up in massive and massless gauge bosons. We leave a detailed computation of such a scenario for future work.

V. DISCUSSION AND CONCLUSIONS
If the Higgs field effective VEV has large non-zero fluctuations during inflation, it could imprint considerable effects on the subsequent stages of reheating, namely resonant particle production (preheating) and perturbative decays from coherent oscillations of the inflaton field. The quantum oscillations in the Higgs field give it a location-dependent effective VEV, imparting mass to any Standard Model (SM) particles to which it couples. If the particle mass exceeds the inflaton mass in some Hubble patch, then reheating there may be delayed, a phenomenon known as Higgs Blocking [33].
Adiabatic fluctuations arise because the Universe exhibits a space-dependent reheat temperature, due to the correspondingly space-dependent Higgs-induced particle masses. Consequently, density perturbations are created that later source CMB fluctuations as well as potentially seed large scale structure. Our scenario differs from the standard paradigm for the generation of density fluctuations in the following way: Unlike the standard case where curvature perturbations for a given scale k are generated at one time and are constant on superhorizon scales, in our case inflaton decay continues to source the perturbations through the end of reheating, leading to growth of perturbations even on superhorizon scales. Here, we have considered two cases: (i) the case of a single-field inflation model in which the inflaton (not a SM field) decays into SM particles coupled through a Yukawa interaction and (ii) the effects of a non-zero Higgs effective VEV on the non-perturbative inflaton decay. For the case of non-perturbative decay, we considered an abelian gauge field coupling to the inflaton through a Chern-Simons term, as found in models of natural inflation [15,79,80].
For perturbative inflaton decay to SM particles (with masses determined by the Higgs VEV), we find that, for high scale inflation with m φ ∼ H I , fermions with SM Yukawa couplings larger than y 0.1 would overproduce temperature fluctuations in the CMB (see Fig. 9), unless one considers low values of the inflaton decay rate Γ 0 10 −4 m φ and correspondingly lowered values of the reheat temperature. For reheating into the top quark, the reheat temperature must be lowered to a point where tension can arise with certain thermal leptogenesis models. For scenarios in which the inflaton decays primarily to one of the majority of SM fermions with smaller Yukawa couplings, a variety of upper bounds can be set on the inflaton decay width similar to those indicated by Eq. (38). However, such constraints can vanish when reheating into the lightest SM fermions and can be significantly relaxed when the scale of inflation is lowered, as shown in Table I. In the case of parametric resonance, the effects of Higgs modulation lead to density perturbations that are significantly larger than those observed in the CMB, even when assuming relatively small couplings, g = 0.1, for SM gauge bosons at the inflation scale. As a result, preheating into Higgsed gauge bosons cannot be the main source of reheating the Universe. If the inflaton couples to both a massless gauge boson (photon) and to massive ones (W ± and Z), the problem becomes highly parameter dependent and is left for future work.
In summary, Higgs-modulated reheating can significantly constrain the parameter space for inflationary models where reheating occurs by inflaton decay to SM particles. Even though quantum fluctuations of the inflaton may produce the observed spectrum of temperature anisotropies in the CMB, any realistic model must also provide for a mechanism to reheat the Universe. Typically considered as independent challenges, we have demonstrated that Higgs modulated reheating could potentially ruin the spectrum of density perturbations produced by quantum fluctuations of the inflaton. We note that specific models of inflation which reheat preferentially into either photons or the lightest SM fermion species are unaffected. Without a concrete inflationary model which is demonstrably able to avoid the series of constraints we have calculated under a set of relatively generic assumptions, the most straightforward way to avoid the effects of Higgs modulated reheating is to introduce additional dynamics into the Higgs sector which stabilize its quantum fluctuations during inflation. Thus, Higgs modulated reheating can also be used as a window into the dynamics of the Higgs during inflation and, potentially, as a probe into physics beyond the SM. the term ω k in Eq. (A3). This is similar to the production of the massive modes from the resonant oscillations of the inflaton field discussed in Sec. II E.
We solve Eqs. (A2) and (A3) numerically during reheating. We find that various resonant bands exist where the W-bosons are produced resonantly, for different values of q W . The corresponding occupation number is then [47] from which we calculate the effective Higgs mass term induced by W-bosons using an approximate expression for the expectation value W 2 as in Eq. (32).