Eternal Inflation and Reheating in the Presence of the Standard Model Higgs Field

We study the details of eternal inflation in the presence of a spectator Higgs field within the framework of the minimal Standard Model. We have recently shown that in the presence of scalar field(s) which allow inflation only within a finite domain of field values, the universe reaches a steady state where the normalized distribution for the field(s) converges to a steady state distribution [1]. In this paper, we analyze this eternal inflation scenario with the renormalized Standard Model Higgs potential, since it also allows inflation in a finite domain, but turns over at high scales due to the running of the self-coupling, marking an exit from inflation. We compute the full steady state distribution for the Higgs using an integral evolution technique that we formulated in [1] and the fractal dimension of the universe. We then obtain a bound on the inflationary Hubble scale in order to have a large observable universe contained within the instability scale of $H\lesssim\mathcal{O}(10^9)$GeV depending on the top mass. Upon reheating of the universe, thermal fluctuations in the Higgs field could potentially pose another problem; however, we compute the rate of thermal bubble production and find that the probability of tunneling in the post-inflationary era is negligibly small even for very high reheat temperatures.


I. INTRODUCTION
The discovery of the Higgs boson [2], which is required to maintain unitarity of W, Z boson scattering, completes the Standard Model (SM) of elementary particles. For the first time we now have a theory that does not violate tree-level unitarity. When combined with gravity, we know that the theory does violate unitarity at the Planck scale, but one wonders at what scale, if any, does new physics enter before the Planck scale. There are reasons to suspect all sorts of possible new physics associated with dark matter, grand unification, strong CP problem, hierarchy problem, baryogenesis, inflation, etc. However, with the lack of discovery of new physics at the LHC, we have to keep an open mind to the possibility that such new physics may only enter at extremely high energies.
We can therefore enquire: what is the absolute highest scale at which new physics can enter? The Higgs sector itself may suggest a potential problem when studied beyond tree level: It is well known that when loop effects are taken into account, the Higgs self-coupling λ becomes negative at energies around 10 10 − 10 11 GeV for the latest top mass bounds [3] (mostly due to its interactions with top quarks and W, Z bosons). This makes the renormalized Higgs potential turn over at some field value h ht and then become negative, meaning there is a potential instability at high field values (see ahead to Fig. 4 (lower panel)). The lifetime of the electroweak vacuum is estimated to be much higher than the current lifetime of the universe [4][5][6][7], and so it may not appear as an immediate problem. However, this instability may be very important in the very early universe when energy densities were extremely high and the Higgs field could explore large field values.
In particular, this instability has potentially huge ramifications for evolution during early universe inflation [8][9][10][11]. This is because during inflation the Higgs field on super-horizon scales roughly undergoes quantum diffusion with step size ∼ H/(2π) per Hubble time, and drifts due to the renormalized potential V . For sufficiently large Hubble scale during inflation H, the quantum diffusion can lead to the Higgs field h becoming larger than the instability scale h ht . In this case the renormalized potential then causes the field to drift down the potential to even large field values, leading to a potential cosmic catastrophe and terminating inflation. In any case, it leads to field values on the "wrong" side of the potential, which is evidently not our observed electroweak vacuum.
This leads to several important questions to be addressed: (i) How frequently does this allow for a large observable universe, that is at-least ∼ 60 e-foldings across in each spatial direction? (ii) What bounds does this imply for the scale of inflation? (iii) What is the global structure of the universe in the context of eternal inflation? (iv) What are the effects of reheating? In this paper we will make progress in addressing these important questions.
Some of these questions have been the subject of considerable investigation in the literature, see Refs. [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. However, the existing works seem to be partially incomplete. For instance, usually a Gaussian ansatz is used in the estimations of the probability distribution of Higgs during inflation which is used to obtain bounds on inflationary Hubble by requiring a large observable universe (∼ e 180 H −3 ) to be born at the end. More importantly, it is often thought that the number of e-foldings provided by the inflaton (typically taken to be N ≈ 60) is essential to estimate the fate of the observable universe, since otherwise for large enough e-foldings and high enough Hubble value the distribution grows indefinitely. However, as we will show in this paper, neither the Gaussian ansatz, nor the requirement of finite number of e-foldings (irrespective of inflationary scales), is needed to describe the physics.
Indeed we will study the system much more systematically: we will consider the structure of the universe in the context of eternal inflation. This leads to results that are independent of the arbitrary choice of initial conditions for the Higgs. Other works in the literature often assume the Higgs begins near the origin and spreads out, then is truncated at N ≈ 60. But as we will see, this is not necessary, as a steady state is established leading us to consider the possibility of much longer phases of inflation, including eternal inflation. The assumption of N ≈ 60 is a conservative one and plausible. However, it is of significant interest to study the possibility that inflation lasted a lot longer and may very well have in fact been eternal. To our knowledge this is the first time such an investigation has taken place. We will show that the universe can obtain a steady state distribution and from this distribution we can determine the probability of obtaining a large observable universe. Along the way we will compute global properties, such as the fractal dimension of the universe. Furthermore, we will move beyond the Gaussian ansatz and carefully track the distribution of different Hubble patches more accurately than in the standard treatments.
At the time of reheating, temperature corrections to the Higgs potential can be important and potentially alter the constraints. For instance, in Ref. [17] it is claimed that high scale reheating can save the universe from the instability; that high scale inflation is allowed, because even though the Higgs can end up on the wrong side of the potential during inflation, it is rescued during reheating due to thermal corrections to the potential. These claims are potentially problematic as it requires extraordinarily fast reheating to prevent the Higgs field from rolling down to a catastrophe, and such problems were not properly taken into account. At the other end of the spectrum, Refs. [15,16] claim that an upper bound on the reheat temperature can be obtained, even if the scale of inflation is relatively low, because they claim the large thermal fluctuations in the Higgs field can cause it to be thrown over to the wrong side, unless the temperature is sufficiently low. However, the authors only considered the Higgs distribution at a single point in space. Instead a correct analysis would be to compute the probability of bubble formation of the Higgs at finite temperature. We will do this and show that even for extremely high reheat temperatures, the tunneling probability is small enough for a true vacuum bubble to form in the observable universe. In this sense, finite temperature SM poses no problems for the observable universe.
In typical models of inflation, it is well known that the universe on large scales undergoes eternal inflation (see [9,11]) and generates infinite numbers of patches in which inflation has ended. In this work we assume that inflation is driven by some other field in some potential, characterized by some scale H. In addition the Higgs field can act as a spectator field (to be clear, we are not assuming the "Higgs-inflation" scenario often studied in the literature, where the Higgs drives inflation, especially because the latest top mass bounds suggests this is unlikely to be viable). Nevertheless the Higgs can act as a field to cause inflation to terminate. In the case of the SM Higgs, inflation is allowed to continue in regions that contain Higgs within its slow-roll end point h end (on the unstable side of its potential), while other regions exit inflation 1 . Then it is important to note that the distribution of Higgs in each Hubble patch ultimately converges to a constant steady state distribution [1] and any initial transient behavior washes out. Therefore, there is no dependence on the amount of e-foldings that the inflaton may provide. All interesting statistical quantities, like the average sizes of regions containing Higgs within h ht , the bound on H to ensure the average size to be at least that of the observable universe, the fractal dimensions of the eternally inflating universe, etc, must be computed at steady state. Finally, in order for our observable universe to be achieved, inflation must end (atleast) locally and there must exist a large (∼ e 180 H −3 ) volume which thermalizes and ultimately settles to the electroweak vacuum 2 . To describe this post-inflationary era, a finite temperature field theory analysis is required and one must compute tunneling probability of the Higgs towards it's true vacuum by forming a thermal bubble.
In this paper, apart from pointing out subtleties and assumptions made in this subject in the literature, our primary objective is two-fold: First, to compute the above mentioned statistical quantities for the eternally inflating universe in the presence of the SM Higgs. For this we make use of an integral (kernel) evolution technique [1] for the distribution of the Higgs (radial degree of freedom in a 4-D Euclidean field space). Second, to describe finite temperature effects after inflation by carefully tracking the rate of bubble formation. Throughout this work, we assume no new physics beyond the minimal SM and work with the (2-loop) renormalized effective Higgs potential with the latest values of relevant parameters [3]; we take central values of all parameters, and scan the top mass in the allowed range: m top = 172.9 ± 0.4 GeV. Also for the sake of simplicity, we assume H = const. throughout this work; this will be justified because the primary bounds will come at rather low scale inflation values, in which H must be very slowly varying (it requires a very tiny ε sr = −Ḣ/H 2 to achieve the correct amplitude of density fluctuations).
The organization of the paper is as follows: In section II, we first derive the kernel that propagates the probability distribution for the magnitude (or "radius") of a multi-component scalar field (relevant to the Higgs doublet) in any Hubble patch. Then to illustrate deviations from Gaussian approximation, we analyze a toy model of 4 fields with a rotationally symmetric "M" shaped potential and we also simulate a 1+1-dimensional universe. We show comparisons of various statistical quantities obtained from actual simulation, with that from Gaussian approximated and actual steady state distribution obtained using the radial field's kernel. In Section III, in order to have a concrete understanding of the situation with SM Higgs in 3+1 dimensional space-time, we simulate 1+1-dimensional universes with the SM Higgs for different choices of H and central top mass. We provide plots of various statistical quantities, and compare them with ones obtained through steady state distribution using radial field's kernel, and also Gaussian approximated ones. Having obtained confidence in the kernel propagation method, we then use this method to obtain the steady state distribution for our primary subject of interest: the SM Higgs in 3+1 dimensions for various H and different top masses, and use them to get the average size of regions with Higgs within h ht , bounds on inflationary Hubble, and fractal dimension of the eternally inflating universe. In section IV, we assume the scale of inflation is low enough that the a large universe can be achieved, and then study the post-inflationary era in a simplified way: we account for temperature corrections to the Higgs potential and calculate thermal tunneling probabilities of the Higgs at different temperatures, and show that the SM provides no instabilities even for very high reheating temperatures. In section V we summarize our results. Finally, in the Appendix we include the relevant SM beta functions.

II. EVOLUTION OF RADIAL FIELD AND STEADY STATE
We are interested in the evolution of a spectator Higgs field during inflation. On super-horizon scales it acquires de Sitter fluctuations, with standard deviation (or "kick") κ that acts as a diffusion per Hubble time and is determined by the Hubble scale (we have κ = H/(2π) in the important case of 3 spatial dimensions, but has different units in other dimensions). We also assume that the field is influenced by its potential V .
To model the evolution of the Higgs field on superhorizon scales, we make use of the usual Fokker-Planck equation. We consider the slow-roll approximation and for now we study general spatial dimensions D. We allow for a multi-component field of length n (the Higgs is a complex doublet with n = 4) ϕ, and call the corresponding probability distribution for the field p( ϕ, N ), where N = Ht are the number of e-foldings. The corresponding Fokker-Planck equation for the evolution of the distribution in any Hubble patch is: This holds only within a field domain that allows infla-tionary energy to dominate and hence inflation to continue. For our purposes we have V = V (| ϕ|) ≡ V (ρ), and the inflationary domain is 0 ≤ ρ ≤ ρ end where ρ end is the slow-roll end point. Also, there exists some instability scale ρ ht in the potential function such that once ρ gets past it, it's classical evolution is to roll down to higher and higher field values which quickly acquire negative energy densities and ultimately around ρ ≃ ρ end inflation ends. As we shall discuss, a useful region to focus on is in fact 0 ≤ ρ ≤ ρ ht , since when we are near the critical Hubble values of interest, the field value ρ ht is actually quite close to ρ end (see ahead to Fig. 9). Solving this PDE requires precise boundary conditions p(0, N ) and p(ρ end , N ), which can be highly nontrivial [30,31], unless the form of potential is highly simplified. Furthermore, it seems impossible to obtain exact solutions for any initial condition p(ρ, 0). This can be dealt with by means of numerical integral evolution with support only within the specified domain [1].

A. Kernel for Radial Field
To make progress, it is useful to discretize the time variable and step the evolution through time by use of a kernel. Upon time-discretization, the kernel K that evolves the distribution through one time step ǫ can be determined by the following chain of reasoning. Firstly, the kernel is implicitly defined by the condition The field itself evolves in a stochastic fashion, governed by the Langevin equation with η i being a set of n Normal random variables i.e. having the following probability density Then one can show the appropriate kernel is with support only in 0 ≤ | ϕ| i ≤ | ϕ| end , and is trivially obtained by replacing η a i for ϕ a i . Also, the derivative of potential (or any function of only the radial field ρ for that matter) is Since it is overly complicated to track all individual components of the field, we can focus on its radial (or magnitude) ρ = | ϕ|. This is all we really need as the Higgs potential is symmetric, and it is only the radial value that determines its central dynamics (i.e., if ρ is small then the field is safe, or if ρ is too large then it can be unstable; see ahead to Fig. 4 We emphasize that in the case of the SM, even though different components of the Higgs doublet are related by gauge redundancies, one must still use a volume measure on the 4 component field space, as we explained in Ref. [32]. In order to obtain the kernel for ρ then, we express the field variables in angular coordinates. With θ and α µ as the polar angle and the other n − 2 angles respectively, the integration measure is n a=1 dϕ a = dρ ρ n−1 dΩ(θ, α).
where dΩ is the differential solid angle. Throughout this paper, we will absorb the radial measure ρ n−1 and total solid angle Ω into the probability distribution p, so that it is normalized as ∞ 0 dρ p(ρ) = 1. Now since the whole dynamics must only depend on ρ at any time step, we can always align ϕ i and ϕ i−1 such that ϕ i−1 points along the polar axis for ϕ i . With this, we have the following simplification: and we can integrate over the solid angle for any n in general. For the relevant case of 4 fields (Higgs doublet), we have with 0 ≤ θ ≤ π, and therefore the radial field's kernel for a 4-field system is with support only in 0 ≤ ρ i ≤ ρ end , and Here I 1 is the modified Bessel function of the 1 st kind. Therefore we have the following evolution equation for the radial field's distribution This is an iterative matrix multiplication technique (upon discretization of field space) and will converge to the dominant eigenstate of the kernel, which after normalization will give the steady state distributionp(ρ)

B. Gaussian Approximation and its Limitations
A concrete way to analyze the Fokker-Planck equation is to study its moments. For example, the equation for the variance σ 2 can be obtained by multiplying the Fokker-Planck eq. (1) with φ 2 and integrating over the fields φ: Similarly one can obtain equations of motion for higher moments.
To simplify things however, very often in the literature a Gaussian approximation is made: . (15) where, as mentioned earlier, we have absorbed the radial integration measure ρ n−1 in the distribution. The only quantity then, is the variance σ 2 which obeys the above simple ODE. An important point to note however, is that for potentials like that of the SM Higgs and especially under a Gaussian ansatz, there exists a critical Hubble H  cr , p (0) grows indefinitely. In any case, people often simplify the analysis by evolving it for some finite e-foldings (typically ∼ 60 in the case of SM Higgs for our universe) and use this to calculate various quantities, e.g. probability of obtaining large volumes with field within the domain 0 ≤ ρ ≤ ρ ht etc. [12][13][14][15][16]. However in either case, the simplification made is significant. First of all it is obvious that due to a finite ρ end , for any H there exists a steady state distributionp(ρ) which cannot have support outside of ρ > ρ end . Then, to leading order the Gaussian approximation can be valid in principle, only as long as σ(N ) ρ end . For the latter case H > H (0) cr , this eventually breaks completely and therefore the approximation is not meaningful for all N . For the former case H < H (0) cr , even though it could be that σ(∞) < ρ end , the non-Gaussianity in the actual steady state distribution could be significant enough for there to be corrections in the estimated quantities. To demonstrate this, we consider the following "M" shaped potential: (ironically, this is upside down from the classical Higgs potential, but it qualitatively captures the RG corrected Higgs potential since it in fact rises, then goes negative at large values, just like this). This leads to the following equation for the variance in the Gaussian approximation which can be re-written as For the sake of comparison with actual 1-D simulations, we work in 1 spatial dimension (although the argument holds for any spatial dimension in general) for which we have κ = 1/ √ 2π [1]. This gives with H (0) cr the above mentioned critical Hubble. The solution is trivially obtained (with σ(0) ≈ 0): where Now for H ≤ H cr , we have an attractive fixed point σ − and thus σ(∞) = σ − dictating a constant distribution 3 . In this case, the Gaussian approximation is valid in principle since we have σ − < ρ end where ρ end marks the slow-roll end point 3 This is true for any σ(0) < σ + and neglecting the effect of the finite-ness of ρ end is warranted, at least in principle, to first order in approximations 4 . However on the other hand for H > H (0) cr , there is no real fixed point and the distribution grows indefinitely σ(∞) = ∞. The Gaussian approximation therefore completely breaks. In reality a finite ρ end would lead to a steady state distributionρ(ρ) which we can easily obtain numerically using the radial field's kernel derived earlier.

C. Comparisons
We now present some numerical results for various important quantities. For a 1+1-dimensional simulation, we have used our network-I as we outlined in detail in Ref. [1]. We begin with one cell with 4 field values all at zero and double the number of cells at each step. The field values in the daughter cells are dictated by Langevin equation (3) with ǫ = ln 2, and are produced only if the radial field ρ in the parent cell was less than the slow-roll end value ρ end . For the kernel evolution, we begin with a delta distribution p(ρ, 0) = δ(0) and evolve it until the normalized distributions at the final step and one previous step (with ǫ = ln 2) differed from each other by only 1 part in 10 3 . In order to converge to steady state in a few e-foldings in our simulations for this case of M potential, we carefully choose our parameters: For H ≪ H cr , we can use Gaussian approximation to estimate the scaling of the time of convergence. From eq. (20), the number of e-foldings to converge to σ − goes like On the other hand for H ≫ H cr , we can use classical equation of motion to estimate the time it takes for the field to roll to ρ end . A straightforward calculation gives the following We model different scenarios keeping λ fixed at 0.06 H 2 while varying ρ ht (in units of κ) such that we always have N * O(10) . For any given scenario, our results are very insensitive to increments in the slow-roll end values. c.f. eq. (22). Note that the Gaussian approximation is so simple that we carry it out in the continuum limit, while the full simulations and kernel evolution are more complicated, so they have a discrete time step, as mentioned above.
To begin, we first show plots of fied distributions for the two cases of H. Fig. 1 (upper panel) compares the In order to force a finite bounded distribution for the latter (otherwise spreading indefinitely), we show two illustrations of truncating the standard deviation of the radial field at the hilltop and the slow-roll end point, i.e. σ(∞) ≡ ρ ht /2 and σ(∞) ≡ ρ end /2 respectively. Fig. 1 (lower panel), on the other hand, compares simulation result against the steady state and Gaussian approximated (σ(∞) = σ − ) distributions for an H < H that we must compute, is the average size of regions that contain the radial field ρ within the hilltop value ρ ht = m/ √ λ. As can be seen from above distribution plots, the non-Gaussianity in the actual (steady state) distribution is important to track. Fig. 2 shows the average length of inflating chains in 1+1-dimensions, computed from steady state distribution obtained from kernel evolution (green curve), and Gaussian approximated distributions (cyan and blue dashed curves), as compared with the actual simulation result (red curves). The former two are obtained assuming completely uncorrelated patches, giving the following Here again for H > H (0) cr , we have the above two considerations of truncating the width of Gaussian approximated distribution at ρ ht /2 and ρ end /2. It is apparent from the figure that not only is the Gaussian approximated result further from the steady state result, the latter is much closer to the actual simulation result especially for smaller H, and gets even closer for larger lengths.
Another interesting quantity to consider is the fractal dimension of chains that are within the instability scale i.e. ρ < ρ ht , defined as where L <ρ ht (N ) is the length of inflating regions with field values within the instability point. With time discretization ǫ and total number of steps s such that  Fig. 3 compares the fractal dimension with ρ < ρ ht obtained from 1+1-dimensional simulations after 12, 15, and 18 2-foldings, respectively, with that from kernel propagation at steady state. The convergence of the simulation towards steady state result obtained from kernel evolution is apparent from this plot.
Having laid out a quantitative analysis for the M potential, all the while elucidating some subtleties and simplifying assumptions made in the literature, and also establishing validity of radial kernel (10) to compute various quantities of interest, we now analyze the eternal inflation scenario with the SM Higgs.

III. ANALYSIS WITH THE STANDARD MODEL HIGGS
The SM potential for the Higgs V = −µ 2 h 2 + λ h 4 /4 is known to pick up corrections from loops. These lead to logarithmically slow changes in λ as we go to high energies. For m top = 172.9 ± 0.4 GeV, the Higgs selfcoupling λ becomes negative at around ∼ 10 10 − 10 11 GeV within the minimal SM.  By replacing the energy scale by E → h, the corresponding renormalized effective Higgs potential can be written as (with G related to wave-function renormalization). The potential is plotted in Fig. 4 (lower panel) on a loglog scale. The point at which this reaches a maximum V ′ (h) = 0, we call the "hill top" value h ht . It is found to be given by There exist ideas in the literature as to how to "cure" the potential by adding new degrees of freedom that alter the Higgs beta functions (such as [33,34]), however we will not pursue that here. Our focus is to put the minimal SM on trial to ever increasing energies and see what constraints on inflation, etc, follow.

A. Gaussian Approximation and its Limitations
For Gaussian approximation, we can proceed similarly as before. Under Gaussian ansatz (15) and the renormalized Higgs potential, we have the following simple differential equation for the variance For all cases when the variance asymptotes to a constant eventually, i.e. H ≤ H (0)) cr , we may set dσ 2 /dN equal to zero giving Thus, we have and for all other Hubble values smaller than this critical value, we have two roots (fixed points) as before with the smaller being an attractor. Numerically, we obtain the following values of critical Hubble for the three values of top mass and D = 3: cr there is no fixed point. However, as we will again see, when studied more precisely beyond the Gaussian approximation, a steady state is realized for any H. However, the sizes of habitable patches will be strongly dependent on the value of H. On the other hand, for H < H (0) cr , the Gaussian approximation does give a steady state, though in the literature the distribution is often just cut off at N ≈ 60 instead.

B. 1+1 Dimensional Simulation with Higgs and Comparisons
Now let us move beyond the overly simple Gaussian approximation towards a full treatment of the problem. In order to validate our analysis concretely, we begin by simulating (using our network-I as mentioned earlier) 1+1dimensional inflating universes with the SM Higgs, for the central value of top mass and different H. With random walk step κ, we scale the field(s) and the potential as giving the following (discretized-) Langevin equation for the fields' evolution in any Hubble patch and the slow-roll end point is obtained as before To mimick some of the aspects of the desired 3+1dimensional scenario, we simulate universes with different numerical values of H ranging from 0.1 h ht to 10 h ht , and simultaneously κ ranging from 0.1 h ht /2π to 10 h ht /2π. Note that with these choices of parameters, it is practically impossible to achieve steady state in simulations since we would require a large number of 2-foldings. Therefore we obtain results from kernel propagation both for 18 2-foldings (which is the maximum amount we go in our simulations), and steady state. To obtain steady state, we propagate the distribution with the appropriate kernel until they start to deviate from each other only by 1 part in 10 3 , as mentioned before. Also for Gaussian approximation and D = 1, we have H (0) cr = 4.54 × 10 9 GeV for the central value of top mass. We again note that the simulation and kernel methods have a discrete time step of ǫ = ln 2, while the Gaussian approximation is in the continuum limit.
First, we present field distributions for both cases H > H The average length of regions containing Higgs within it's hilltop value, is shown in Fig. 6 (upper panel). As before, under the assumption of uncorrelated patches, average length obtained using steady state distribution is given as before (c.f. eq. (25)) with the appropriate kernel for the Higgs potential. Note that even though the actual average size (red curves) is bigger in reality (due to correlations between nearby Hubble patches), the assumption of uncorrelated patches (solid green curve) gets closer to the former especially for very large sizes. Finally, the fractal dimension of these regions containing Higgs within the hilltop given as before (c.f. eq. (40) and (26)), is shown in Fig. 6 (lower panel).

C. Full Analysis in 3+1 Dimensions
Having shown the validity of our kernel evolution technique concretely, we now move on the important case cr ) in 1+1dimensions. Solid red and green curves are from actual simulation and kernel propagation after 18 2-foldings respectively, dashed green curve is the steady state distribution obtained from kernel propagation, while dashed blue curve is the Gaussian approximation with σ(∞) ≈ 2.39 × 10 9 GeV. of 3 + 1-dimensions. Using the radial kernel (10) with D = 3 and κ = H/2π, and ǫ = ln 1.5, we numerically obtain steady state distributionsp(h) for different Hubble values and the three different top masses. Fig. 7 shows some of these distributions for the central value of top mass with different H, on a linear (upper panel) and loglog (lower panel) scale, respectively.
We can use the kernel method then to again compute average volumes, trusting that the assumption of uncorrelated Hubble patches is not a terrible approximation. In general as a function of Hubble H, Fig. 8  as before for the three different values of top mass. We would like the average volume to be at least as big as our observable universe. Now note that we need ≈ 60 e-foldings of inflation to produce our observable universe. This leads to a post-inflationary era made up of ∼ e 180 inflationary Hubble patches. If we then expand that co-moving volume to today it will indeed be on the size of our observable universe. Hence the patches created in this steady state distribution should typically be at least this large, and ideally the average < Volume > should be at least this large. We find that the bound on inflationary Hubble values, in order to achieve this is which is one of our primary findings.
We note that these values are less than the Gaussian approximated critical Hubble values H and N = ǫ s. As mentioned earlier, in any realistic scenario of our universe, after inflation it must reheat rendering thermal corrections extremely important in the analysis. Having dealt with the eternal inflation scenario at zero temperature, in the next section we will analyze these thermal corrections.

IV. THERMAL EFFECTS AFTER INFLATION
It is known at finite temperature, the Higgs potential is partially improved, in the sense that at high temperatures the new effective potential has an instability scale h ht that is pushed to even higher values. In the literature it has sometimes been assumed that this can relax the bounds on the inflationary Hubble scale completely, i.e., any H is now allowed, since even though the Higgs may be on the wrong side during inflation, it can be thrown back after inflation in the altered potential. Often authors assume "instant" reheating, and some have used related ideas to build interesting models to explain various physical phenomenon (e.g. see [17]).
However there are reasons to be skeptical that the Higgs field can be "saved" during the reheating era. In Fig. 9 we plot a measure of slow-roll, namely V ′′ /H 2 (which is on the order of the second slow-roll parameter η sr ). This is the effective mass-squared in units of Hubble. This is evaluated at H = H max that we found above, and a similar curve applies for nearby value of H too. The dashed lines are the turn-over scale h ht at zero temperature (during inflation). We see that for h > h ht then |V ′′ |/H 2 Max is appreciable, so the field is expected to start to roll-quickly. This gives very little time for reheating to occur and "save" the Higgs field. For H > H max , the field does not roll as fast, but this is an even more precarious situation since so much of the Higgs field is on the wrong side, and so it would be rather non-trivial for it all to be saved. In summary, it seems safe to assume that H < H max that we computed above. Although a full exploration of this issue is left for future work.
On the other hand, it was claimed in Ref. [15] that even though h ht is pushed to higher values, thermal fluctuations in the Higgs can be so large and can throw the Higgs to the wrong side, even if H < H max . This may be a particularly interesting approach when one is considering the formation of AdS regions where gravity will ultimately play a significant role. However, it is also of importance to note that no physical spatial scales were accounted for in this analysis or a direct inclusion of local gravitational effects. Thermal fluctuations were calculated at a point in space, which was then used to construct a Gaussian probability distribution in the Higgs. This was then required to have small enough support for h > h ht such that a large ∼ e 180 H −3 sized region can survive, and a bound on the reheat temperature was found. However, this does not ensure an instability. The correct analysis is to consider the formation of thermal bubbles; these have finite spatial extent (roughly given by the inverse effective Higgs mass), and one needs to take into account the spatial dependence of the two-point correlation function. The appropriate calculation is the corresponding finite temperature Euclidean bounce solutions (instantons). Taking this into account, we will be led to different conclusions from Ref. [15]. However, local gravitational effects could still alter the situation, and a full inclusion of such effects is an interesting topic.  We also mention there is important earlier work on this thermal tunneling (including Refs. [6,7,[35][36][37][38]), which is connected to our work here.
We ignore the possibility of a direct coupling between the inflaton and the Higgs, leaving this for future work, and assume a generic reheating scenario that takes at least a few e-foldings. One may imagine reheating takes place by coupling the inflaton to gluons or other SM particles by higher dimension operators. Once the universe gets in thermal equilibrium with all the SM degrees of freedom, we can then attach the thermal free energy to the effective Higgs potential: The free energies due to bosons and fermions in the thermal bath are respectively (we have subtracted out the radiation piece ∼ T 4 , since, although it is important to the total energy of the universe, it does not directly affect the Higgs dynamics as it is h-independent). The only significant contributions come from the W and Z bosons (6 and 3 total degrees of freedom respectively), and the top quark (12 total degrees of freedom), each of which have the fol-lowing masses and where g, g ′ , and y t are the SU(2), U(1) and top Yukawa couplings respectively. The effective temperature corrected Higgs potential is therefore With this new effective potential, we can calculate bounce actions for tunneling of the Higgs to true vacuum at finite temperature. In the next section, we numerically obtain bounce actions for different temperatures with the above temperature corrected effective Higgs potential, and show that there is no bound on reheat temperature.

A. Tunneling Probability at Finite Temperature
For a field theory at finite temperature in 3+1dimensions, the relevant bounce action associated with a thermal bubble is provided by the Euclidean O(3) solution [39] S B (T ) = 4π dr r 2 1 2 where h B (r, T ) is the bounce solution for the equation of motion with boundary conditions Then, the probability of tunneling via thermal bubble formation is roughly (we suppress the pre-factor here as it is of sub-leading importance) With V eff (h, T ), we numerically obtain bounce actions for different temperatures and hence tunneling probabilities.  This probability, when properly normalized, is per unit volume per unit time. If we take the inverse of this probability, it is seen to be far greater than the number of bubble sized patches in our comoving universe (that is very roughly ∼ e 180 ∼ 10 78 , or so). It is clear that not only is there no bound on reheat temperature, even extremely high temperatures approaching Planck scale (though they would not realistically be reached in any post-inflationary era) carry insignificant tunneling probabilities for there to be a thermal Higgs bubble formation in our observable universe.

V. SUMMARY
In this paper, we have for the first time analyzed the eternally inflating scenario in the context of the minimal SM. Since the SM Higgs field develops an instability at high field values, thereby allowing an inflationary phase only in a finite domain, the universe gets to an eternally inflating steady state in which any initial transient behavior washes out, and all interesting quantities like average volumes of regions, fractal dimensions, etc, are therefore independent of the number of e-foldings provided in the final stages of slow-roll inflation.
The corresponding steady state distribution in the field has significant deviations from Gaussianity depending on the value of inflationary Hubble H and is therefore important to keep track of. In order to obtain it for different choices of parameters, we derived a kernel that propagates the radial field's distribution in time, which upon normalization gives the desired steady state distribution. In order to compare it against Gaussian approximation, we began with a toy model of 4 scalar fields with a spherically symmetric M shaped potential in 1+1 dimensions and showed that the usual Gaussian approximated (radial) field distribution in any Hubble volume deviates significantly from the former. For a concrete analysis, we also compared them against actual 1+1-dimensional simulations, finding the kernel method to be much more accurate. We then used the radial field's kernel to analyze this eternal inflation scenario with the 2-loop renormalized SM Higgs potential. As before, we first compared kernel propagation analysis with actual 1+1-dimensional simulations, and showed that the former is accurate in calculating various statistical quantities. Then, we applied it to 3+1-dimensions and along with obtaining these statistical quantities, we obtained upper bounds on the inflationary Hubble scale H to obtain an observable universe (of size ∼ e 180 H −3 ) for upper, central, and lower value of top quark mass.
Finally, after inflation ends locally in various regions, they must reheat and thermal corrections (due to W, Z and top quark) become important, leading to the instability scale of Higgs getting pushed to higher values for larger temperatures. Assuming a generic reheating scenario, where the inflaton takes at least a few Hubble times to reheat the universe, and therefore thermal effects are unlikely to rescue the Higgs if H > H max , we then calculated tunneling probabilities from the electroweak vacuum at finite temperature to the true vacuum. We found that even for extremely high temperatures, the bounce actions are large enough that thermal bubble creation is negligibly rare in our comoving volume. Therefore, high reheating of the SM poses no threat to the electroweak vacuum metastability all the way to M pl , but a high scale of inflation does. Modified reheating scenarios will be the subject of future work.