Higgs Vacuum Stability with Vector-like Fermions

We present the effects of vector-like fermions (VLF) on the stability of the Higgs electroweak vacuum, using the renormalization group improved Higgs effective potential. We review the calculation of the one-loop beta-functions of the standard model couplings paying particular attention to the fermion contributions. From this, we derive the VLF contributions to the beta-functions. Using these beta-functions we determine the scale at which the effective Higgs quartic-coupling becomes zero and goes negative, signaling vacuum instability. We find that for certain VLF masses and Yukawa couplings, the Higgs quartic stays positive for field values all the way up to the Planck scale, implying that the meta-stable vacuum of the standard model can be rendered absolutely stable if VLFs are present with certain parameters. For other values of VLF parameters, the Higgs vacuum is metastable as in the standard model. For cases where the vacuum is metastable, we compute the probability of quantum tunneling from the false electroweak vacuum into a deeper true vacuum in our Hubble volume by numerically solving for the bounce configuration in Euclidean space-time and computing the bounce action for it. We compare our numerical solution with the analytical approximation for the bounce action commonly used in the literature, and comment on when the latter may be used.


Introduction
The stability of the electroweak (EW) vacuum can be studied using the Higgs effective potential (for a review see Ref. [1]). Recent investigations (see for example Refs. [2,3,4]) at the NNLO level have revealed that the Higgs vacuum is metastable in the standard model (SM), with the life-time in the false (EW) vacuum being much much larger than the age of the universe. This situation arises because the Higgs quantum effective potential V eff (h) has a smaller value for h ∼ 10 10 GeV when compared to its value at the EW vacuum expectation value (VEV) v ≈ 246 GeV, i.e., V eff (h ∼ 10 10 GeV ) < V eff (v) for the SM.
There are many compelling reasons to expect physics beyond the standard model (BSM). These include theoretical reasons such as the gauge hierarchy problem, and observational reasons, such as neutrino mass generation, dark matter, and generation of the baryon asymmetry of the universe. A plethora of BSM extensions have been proposed to address these shortcomings of the SM. These inevitably add new particles to the SM particle content. In particular, resolution of the gauge hierarchy problem necessarily have new states coupled to the Higgs. In such cases, the above conclusions on Higgs EW vacuum stability must be revisited by including the effects of such new particles. Many such BSM extensions include vector-like fermions (VLF) that couple to the Higgs, and are often the lightest BSM states. They therefore have a significant effect on EW vacuum stability. Some examples of such models that include vector-like fermions are in the following contexts: the gauge hierarchy problem such as AdS-space/composite-Higgs models in Refs. [5,6,7,8], Higgs-portal dark matter models in Refs. [9,10,11,12,13], gauge-coupling unification in Refs. [14,15,16,17,18], neutrino mass generation and vacuum stability in Refs. [19,20], a combination of these in Ref. [21,22,23], and effective models in Refs. [24,25]. Motivated by these considerations, we study the effect of VLFs that are coupled to the Higgs on EW vacuum stability. 1 Many of these models may also contain new bosonic states apart from VLFs. In such a case, a full conclusion about the stability of the EW vacuum in that model can be reached only after including the contributions of these bosonic states also. However, fermions usually have the biggest role in destabilizing the EW vacuum, and so our analysis here addresses the most crucial ingredient in this problem. Hence, our goal here is to analyze model-independently the generic effects of VLFs on EW vacuum stability.
We set the stage for our analysis by writing the classical Higgs potential as Including quantum effects, we can write the quantum Higgs effective potential as where λ eff has a dependence on h of the form ln (h/M ) with M being a subtraction scale. For h m h (the physical Higgs mass m h ≈ 125 GeV), the mass term has a negligible effect, and thus, to an excellent approximation, we can write 2 Renormalization group improved Higgs effective potential We have in the SM the Lagrangian density, showing only the terms relevant to our analysis here, where the · represents the antisymmetric combination in SU (2) space, q L = (t L b L ) T is the SU (2) doublet, and, t = (t L t R ) T and b = (b L b R ) T are the top-quark and bottom-quark Dirac fermions. It is sufficient for our purposes to keep only the top Yukawa coupling y t in the SM as the others are much suppressed. Next, we present the SM β-functions and extend them to include the VLF contributions.

SM RGE
We first discuss the SM RGE β-functions at the 1-loop level and include some significant 2-loop effects. We use the SM RGE to find the Higgs field value at which λ(µ) becomes zero and compare our results with those in the literature. Denoting the relevant SM couplings generically as κ i = {λ, y t , g 3 , g 2 , g 1 }, the RGE are of the form We derive the fermion contributions to β κ in Appendix A since our goal in this work is to extend them to include VLF contributions. We take the other terms from the literature (see for example Ref. [3]). Putting these together, the 1-loop β-functions, β with b a = (−7, −19/6, 41/10) for g a = (g 3 , g 2 , g 1 ) respectively, and N c = 3 for a fermion in the fundamental representation of SU (3). For g 1 , we use the SU (5) normalization, i.e. the SM hypercharge gauge-coupling g is related to g 1 by g 1 = 5/3g . The precision of the full 2-loop (or higher order) calculations that are available in the literature are not required for our purposes since our goal is to analyze BSM physics contributions that involve as yet experimentally undetermined parameters. However, to help compare our numerical results to what has been obtained in the literature for the SM, we will include 2-loop SM contributions to the β-functions that depend on y t and g 3 as they are numerically the most significant. They are (see for example Ref. [3]) We use these RGE to determine the Higgs field value µ at which λ(µ) becomes zero, signalling vacuum instability. This will be discussed in Sec. 2.3. We discuss next the VLF contributions to the RGE.

VLF contributions to the RGE
We add an SU(2) doublet VLF χ = (χ 1 χ 2 ) T and an SU(2) singlet VLF ξ and couple it to the Higgs as follows where the · represents the antisymmetric combination in SU(2) space. 2 Extracting the Higgs interactions from this yields  [27]). For simplicity, in this work, we do not turn-on such mixed Yukawa couplings; an analysis including such mixed couplings will be the subject of future work.
We derive the 1-loop RGE contributions due to the VLF (see Appendix B for the derivation) and add them to the SM contributions given above. The VLF contributions to the RGE in Eqs. (7)-(9), and the SM and VLF contributions to the RGE for the new couplingỹ, are where n 3 is the number of colored VLF SU(3) triplets, i.e. VLQs, n 2 is the number of SU(2) doublets, n 1 is the number of SU(2) singlets, n F is the number of complete VLF families coupled 2 If another SU(2) singlet VLF ζ is added, we can add the terms After adding the ζ, the one doublet and two singlet VLF structure then mimics the SM quark or lepton structure in a generation. For keeping the field content minimal, we will omit the ζ in our work here, and therefore will not include theỹ2 term.
to the Higgs (a family is a doublet and a singlet both present), andn V LQ F = 1 if the VLF is a VLQ family or zero otherwise. For example, n F = 0 for either VLF singlets or doublets added (but not both), and n F = 1 for one SU(2) doublet and one singlet VLF added together such that a Yukawa couplingỹ can be written down with the Higgs. Only VLQs contribute to β g 3 and VLLs do not. For instance, for one VLQ family of χ and ξ, we have N c = 3, n 3 = 3, n 2 = 1, n 1 = 1, n F = 1 and n V LQ F = 1. Our goal in this work is to analyze the stability of the electroweak (EW) vacuum for which the behavior of V eff at large field values is most important. We have therefore not kept the finite mass effects in the RGE as they are small, being of the form (m/µ) for µ m where m collectively denotes the particle masses. We include the VLF contributions only for µ ≥ M V L , where M V L is the vector-like fermion mass.

RGE Numerical integration results
We take the input parameters as follows and as compiled in Ref. [3], with the renormalization point taken as the top mass scalem t : The EW VEV: v = 246. The U (1) coupling constant:g 1 = 5/3 g = 5/3 × 0.35937 (NLO).
In terms of these inputs, we set the top mass to bem t =ỹ t v/ √ 2 and the Higgs massm h = 2λ v. In this work, since our interest is in analyzing a new physics (VLF) model with unknown parameters, the full precision to which these are defined are not so important, and the above specification is more than adequate for our purposes.
The RGE are a coupled set of first order differential equations for the couplings λ(µ), y t (µ), g 3 (µ), g 2 (µ), g 1 (µ). We take the inputs given above atm t and integrate the RGE numerically, including both the SM contributions in Section 2.1 and the VLF contributions in Section 2.2. As already mentioned, we include the VLF contribution only for µ ≥ M V L .
In Figs. 1-4 we show the evolution of the different couplings for the SM and also for some representative VLF cases. In these figures, the dashed lines are for the SM with only the SM particle content with no VLFs, while the solid lines are for various VLF cases. As is evident, in the SM all the couplings decrease with field value h ≡ µ. The Higgs quartic coupling λ becomes zero and goes negative at about µ ∼ 10 10.5 GeV, while all the other couplings stay positive all the way up to M Pl . It is interesting that β λ approaches zero for large field values (cf Fig. 11). In the following, we discuss the evolution of the couplings in the presence of VLFs.
In Fig. 1 we show the evolution of the couplings with field value h ≡ µ for (n) degenerate SU (2) singlet VLQs for various M V L , where the M V L values are shown in the notation (nEm) ≡ n × 10 m GeV. When 3 or more degenerate singlet VLQs of mass 3 TeV are added, interestingly, λ never goes negative, unlike in the SM. When we add only singlet VLQs, SU(2) invariance forbids a coupling of the Higgs to such VLQs, (we do not turn-on mixed Yukawa couplings between SM fermions and VLFs as we noted earlier). However, these VLQs contribute to β g 3 and also possibly     Figure 1: The evolution of λ, y t , g 3 with Higgs field value µ for (n) number of degenerate SU (2) singlet VLQs of mass 3 TeV (first three plots), and λ with 3 degenerate singlet VLQs of mass 3 × 10 3 GeV, 10 5 GeV and 10 7 GeV shown respectively as 3E3, 1E5 and 1E7 (last plot).    Figure 2: The evolution of λ, y t , and g 3 with Higgs field value µ for (n) number of degenerate SU (2) doublet VLQs of mass 3 TeV (first three plots), and λ with a doublet VLQ of mass 3 × 10 3 GeV, 10 5 GeV and 10 7 GeV shown respectively as 3E3, 1E5 and 1E7 (last plot).
β g 1 if the VLQ has hypercharge, and because of the coupled nature of the RGEs, λ(h) does see the effect of the VLQ. y t and all g a decrease as µ increases. We discuss the implications of this to vacuum stability in Sec. 4. In Fig. 2 we show the evolution of the couplings with field value h ≡ µ for (n) degenerate SU (2) doublet VLQs for various M V L . We observe that when we add one or more doublet VLQs of mass 3 TeV, λ never goes negative for the same reasons as above. We also see that adding a doublet VLQ with mass up to about 10 5 GeV will have this feature. If we add five or more doublets with 3 TeV mass, we find that y t becomes negative at around 10 16 GeV invalidating the theory beyond that scale. y t and all g a decrease as µ increases. We discuss in Sec. 4 the implications to vacuum stability of λ remaining positive.
In Fig. 3 we show the evolution of the couplings with field value h ≡ µ for a degenerate family of one SU (2) doublet VLL and one singlet VLL for various M V L andỹ. The VLL Yukawa coupling values are shown as (ỹ) and the M V L values are shown in the notation (nEm) ≡ n × 10 m GeV. For y(M V L ) = 1, we see thatỹ increases as µ increases. y t eventually starts increasing at large µ, which is a behavior unlike in the SM. We notice that the scale at which λ becomes negative decreases as y increases, or as M V L decreases.
In Fig. 4 we show the evolution of the couplings with field value h ≡ µ for a degenerate family of one SU (2) doublet VLQ and one singlet VLQ, for various M V L andỹ. We see that forỹ(M V L ) = 0.5, y decreases as µ increases. For M V L = 3 TeV, ifỹ > 0.35, the scale at which λ becomes negative decreases asỹ increases or as M V L decreases and is lesser than in the SM, while ifỹ < 0.35, the scale at which λ becomes negative is larger than in the SM. In fact, forỹ < 0.3, λ stays positive all the way up to M P l . We see that whenỹ = 0.1 for example, λ stays positive all the way up to M P l for M V L up to about 10 5 GeV.      These examples illustrate a range of effects on the evolution of the couplings due to VLFs. For the cases when λ does go negative, the EW vacuum is not the absolute minimum but is meta-stable. There is then a non-zero probability that the EW vacuum will tunnel quantum mechanically away to those (large) field values where V eff < 0. We turn next to an analysis of this possibility and a computation of the tunneling probability.

Tunneling away from EW vacuum
In this section, we compute the tunneling probability from the metastable electroweak Higgs vacuum into a deeper true vacuum via quantum mechanical barrier penetration. We do this by computing the Euclidean action for the Bounce configuration of the Higgs field (for a review see for example Ref. [28]). We compute the bounce configuration using the running couplings that we presented in Sec. 2. From the Bounce action S B , we compute the tunneling probability P tunl .

Method of computing the tunneling probability
We briefly review here how to compute the bounce configuration and the tunneling probability (for details see for example Refs. [1,28] and references therein). In Sec. 5 we discuss in detail why we do not use the approximation commonly used in the literature, but resort to actually solving the bounce EOM numerically as described in this section.
Let us recall that the Lagrangian density L and the action S for the Higgs field in Minkowski coordinates are with the effective potential defined in Eq. (2). We define the Euclidean time τ = it, the Euclidean and write the Euclidean action as where ∂ i is with respect to the Euclidean coordinates ρ i . As mentioned earlier, taking V eff (v) = 0 at the EW minimum, the central question of interest in our work here is whether the EW vacuum is absolutely stable or if there is a possible transition to other field values, which would be possible only if V eff (h) < 0 for some h typically much larger than v. If there exists field values for which V eff (h) < 0, the EW vacuum at h = v is typically separated from this by a (large) barrier and a vacuum transition can only occur via quantum tunneling. In such a situation we would like to know the time-scale of the tunneling in comparison to the age of the universe and see if we can gain an understanding of why the universe has not tunneled away to the true vacuum with h v, but is in the EW vacuum today. If for some large field value, h ≡ σ say, suppose V eff (σ) = 0, and suppose V eff is negative for h σ. (V eff may turn around and have a second minimum (or not) for h > σ depending on other BSM contributions in the RGE.) Equivalently, from our definition of the field dependent coupling in Eq. (3), suppose λ(σ) = 0, and that λ(h) is negative for h σ. The vacuum configuration defined to have total energy E = 0 at h = v can quantum mechanically tunnel to h ≥ σ with V eff (h) < 0. If the vacuum were to tunnel so, the field then runs down the potential classically toward large field values h > σ. The tunneling probability is given in terms of the bounce configuration [28] which satisfies δS = 0, This configuration is a solution of the equation of motion (EOM). In Euclidean coordinates, the EOM reads We look for an O(4) symmetric solution [29], which implies that it depends on ρ, i.e. h(ρ i ) = h(ρ).
The EOM then reads with the boundary conditions (BC) (dh/dρ)(ρ = 0) = 0 and h(ρ → ∞) = v. We must also have h(ρ = 0) = h 0 ≥ σ for this to represent tunneling. This EOM is identical to that of a classical particle moving in a potential −V eff with a "friction" term present (second term in Eq. (25)) that dies-off as 1/ρ as ρ increases.
In Euclidean space-time, the bounce configuration h B (ρ), will have the feature of a fairly sharp transition in ρ from h 0 to v. In Minkowski space-time, this configuration looks like an expanding bubble with the bubble-wall separating a region of true-vacuum inside and the false EW vacuum outside. The bubble nucleation probability per unit 4-volume is given by [ where we have included a prefactor of M 4 on dimensional grounds with M an appropriate mass scale, ∆V 4 is a unit 4-space-time volume, and S B is the Euclidean action for the bounce configuration h B (ρ) given by We make the choice M 4 = V eff (h 0 ) since h 0 is typically the largest scale in the problem and gives the largest tunneling rate, and hence the most conservative bound on the allowed VLF parameter-space from vacuum tunneling. If a bubble bigger than the critical size is nucleated anywhere in our past light-cone it would have engulfed us by now and we would not find ourselves in the EW vacuum now. The (dimensionless) volume of our past light cone is about V 4 ∼ (1/m 4 t ) exp (404), which is nothing but our Hubble 4-volume in 1/m 4 t units, and we choose this unit since our starting point for the running is at the m t scale. Thus, the total probability that we would have nucleated a bubble in our Hubble volume and tunneled into the true vacuum by now is P tunl ∼ (dP tunl /dV 4 )V 4 , which gives [1] If P tunl 1 for the given V eff , we deem this as an acceptable situation. In other words, if P tunl 1, we take this to mean that the probability that we would have tunneled into the true vacuum due to a bubble nucleating in our past light-cone is essentially unity, and therefore the model that generated that V eff we consider is disfavored. Evidently, the larger S B is, the smaller P tunl is, and the latter is exponentially suppressed by S B . In the following section, we solve numerically the bounce EOM to get the bounce configuration h B (ρ), compute S B for this h B , and compute P tunl . We do this for the SM and some VLF extensions.

Tunneling probability numerical results
Here we describe the method we use to solve the bounce EOM numerically, obtain the bounce configuration h B (ρ), the Bounce action S B and the tunneling probability P tunl for the SM and various VLF extensions.
The value of S B largely depends on the behavior of V eff at large field values h where the h 4 term dominates, which is why we included only that term in V eff . Nevertheless, for completeness, we insert an EW minimum at v by including it in λ(h) for h ∼ v as follows. Keeping the h 2 term, we have for . This definition implies that we have the effective quartic coupling given by Slightly above the scale v, we match this to the λ(h) obtained by solving the RGE. In this way we effectively obtain a minimum at h = v while the larger field value evolution is governed by the RGE. We add a constant term to V eff and make V eff (v) = 0.
With theĥ( ) and (dĥ/dρ)| got above at the initial BC, we then numerically integrate the EOM in Eq. (25) forρ ∈ ( ,ρ end ) and obtainĥ B (ρ) over this domain. The large values of the fields and the presence of the friction term complicates the numerical implementation. A further challenge is that satisfying the required end-condition requires an extremely sensitive tuning of the starting valueĥ 0 . By an iterative search algorithm we are able to obtain the bounce configuration h B (ρ) using Mathematica.
Piecing together the analytical solution above and the numerical solution, we obtain the bounce configuration over the complete domainρ ∈ (0,ρ end ). Following this procedure, we present below the bounce configuration, the bounce action evaluated for this bounce, and the tunneling probability for the SM and various VLF extensions.
For the SM, the λ(h ≡ µ), V eff (h ≡ µ) and the bounce configuration h B (ρ) obtained numerically are shown in Fig. 5. The V eff (µ) is positive for smaller µ, crosses zero at about µ ≈ 10 10.75 GeV, and is negative for larger µ. The blue dot shows the starting field value (h 0 ) of the bounce and the red dot the ending field value (v). For this bounce, we find by numerical integration of Eq. (26) that the value of the Euclidean bounce action is S B = 2866 (in = 1 units). From this, we compute the tunneling probability into the true vacuum in our Hubble volume from Eq. (27) to be P tunl ∼ 10 −1013 , which is an incredibly small probability. This, and many other comparisons we have done for the SM, are in excellent agreement with the results obtained in Ref. [3].
Next, we solve the bounce EOM, and compute the S B and P tunl for various VLF representations. We start with a (color singlet) VLL family with SM-like hypercharge assignment present, i.e. an SU(2) singlet with hypercharge −1 and an SU(2) doublet VLL with hypercharge −1/2 both present,   Fig. 6. The V eff (µ) is positive for smaller µ, crosses zero at about µ ≈ 10 6.5 GeV, and is negative for larger µ. The blue dot shows the starting field value (h 0 ) of the bounce and the red dot the ending field value (v). For this bounce configuration, we find S B = 467 and P tunl ∼ 10 −3 . This parameter-space point is thus acceptable as the tunneling probability into the true vacuum is sufficiently small for us to understand why the electroweak vacuum has still not tunneled away into the true vacuum within the age of the universe. That is, for this model with VLL present, the probability of a true vacuum bubble having nucleated in our past light-cone is sufficiently small, although this probability is much much larger than in the SM. We see that the presence of a VLL increases the tunneling probability dramatically compared with the SM. As another example, we consider a VLL family with M V L = 10 3 GeV andỹ = 0.61, for which we find S B = 429 and P tunl ∼ 10 14 . This large value implies that the EW vacuum could have tunneled into the true vacuum in our Hubble volume essentially with unit probability. Therefore, this parameter-space point can be considered severely disfavored. These two examples also show that P tunl is extremely sensitive toỹ, with a small change of 0.01 inỹ between the two cases resulting in a change of S B of 38, which in turn results in a P tunl 17 orders of magnitude different because of its exponential dependence on S B . As another example, for a VLL family with M V L = 10 5 GeV andỹ = 0.59, the bounce configuration is shown in Fig. 7. For this bounce configuration, we find S B = 482 and P tunl ∼ 10 −1 which is acceptable. For a VLL family with M V L = 10 7 GeV and y = 0.6, we find S B = 560 and P tunl ∼ 10 −32 . For a VLL family with M V L = 10 7 GeV andỹ = 0.62, the bounce configuration is shown in Fig. 8. For this bounce configuration, we find S B = 491 and     Figure 9: For a VLQ family with M V L = 3 × 10 3 GeV andỹ = 0.57, the running coupling and the effective potential as a function of the field value h = µ, and the bounce configuration h B (ρ). The blue (red) dot shows the starting (ending) value of the bounce. P tunl ∼ 10 −1 which is acceptable.
The hidden sector Higgs-portal dark matter model of Refs. [10] essentially behaves like a VLL family considered above, for the following reason. Although in this model the VLF dark matter is a singlet and does not couple directly to the Higgs, due to the Higgs mixing with a hiddensector scalar, a coupling with the Higgs is induced with sizeỹ ≡ κs h , where the right-hand-side is in the notation of that paper and involve the parameters of that model. As can be inferred from the analysis in Ref. [10], we requireỹ 1 to keep the direct-detection rate small in order to honor experimental constraints. Thus, from the results above, we infer that EW vacuum stability constraints are not too severe in such models.
Next, we compute S B and P tunl with a color triplet VLQ family with SM-like hypercharge assignment present, consisting of an SU(2) singlet VLQ with hypercharge 2/3 and an SU(2) doublet VLQ with hypercharge 1/6 both present, for various common mass M V L and variousỹ. For a VLQ family with M V L = 3 × 10 3 GeV andỹ = 0.5, we find S B = 806 and P tunl ∼ 10 −147 , which is comfortably acceptable. For a VLQ family with M V L = 3 × 10 3 GeV andỹ = 0.57, the bounce configuration is shown in Fig. 9. For this bounce configuration, we find S B = 469 and P tunl ∼ 10 −4 , which is allowed. For a VLQ family with M V L = 10 5 GeV andỹ = 0.54, we find S B = 513 and P tunl ∼ 10 −19 , which is comfortably acceptable. For a VLQ family with M V L = 10 7 GeV and y = 0.6, we find S B = 474 and P tunl ∼ 10 −2 , which is acceptable.

Second minimum in V ef f
Thus far, we have investigated the situation when only the VLF is present and the effective potential has only a minimum at v and no second minimum at large field values but rather runs off in a bottomless manner. If the VLF is accompanied by other states, presumably in a UV completion that it is a part of, one can contemplate the possibility of the potential being turned around due to the contributions of the extra states and the appearance of a second minimum at large field values. We encode this possibility by adding a second minimum in the effective potential as shown in Fig. 10 for the case of a VLQ family with M V L = 3 × 10 3 GeV andỹ = 0.75. The λ(h ≡ µ), V eff (h ≡ µ), and the bounce configuration for this modified potential are shown in Fig. 10. The V eff (µ) is positive for smaller µ, crosses zero at about µ ≈ 10 4.6 GeV and becomes negative, obtains a minimum at about µ ≈ 10 5.1 GeV, crosses zero again and becomes positive for larger µ. The blue dot shows the starting field value (h 0 ) of the bounce and the red dot the ending field value (v). For this bounce configuration, we find S B = 3107 and P tunl ∼ 10 −1166 , which is an incredibly tiny tunneling probability and very comfortably acceptable. The reason why this parameter-space point which was excluded if no second minimum was present as found earlier, is now allowed if a second minimum present is as follows. In the bounce configuration for this situation, the field value starts close to the second minimum, and stays there for a substantial amount of time (i.e. ρ) near the minimum since dφ/dρ is small there, and as a result the friction term starts reducing in significance due to its 1/ρ behavior. Since the friction term becomes small, the field value can overcome the barrier and reach v. (It can even overshoot v leading to an imaginary solution if the initial value of the field is chosen too big.) Therefore, the starting field value is much lesser compared to the earlier case without the second minimum, and we find that the resulting S B is much larger and P tunl much smaller, now allowing a parameter-space point that was excluded earlier.

Absolute stability of the EW vacuum
We have seen that the EW vacuum is meta-stable in the SM as there is a deeper minimum below the EW vacuum albeit shielded by a potential barrier, and due to the tunneling probability being incredibly small, the life-time of the metastable vacuum is extremely large compared to the age of the universe. In Section 3 we added VLFs and analyzed regions ofỹ and M V L parameter-space for which there again is a deeper minimum making the EW vacuum metastable. We computed the tunneling probability and found that in some regions of parameter space, P tunl is acceptably small while in others it is unacceptably large. In this section, we highlight VLF cases where the addition of VLFs makes the EW vacuum the global minimum, rendering it absolutely stable.
Consider first adding some number of either SU(3) singlet VLQs or doublet VLQs, but not both. For instance, we showed in Section 2.2, Fig. 1, that when 3, 4 or 5 SU(2) singlet VLQs all with 3 TeV mass are added, λ(h) never goes negative, implying that the EW minimum is the global minimum and absolutely stable, unlike the SM situation. As we show in Fig. 2, the same conclusion holds also when we add one to four SU(2) doublet VLQs with a 3 TeV mass, or one doublet with mass up to about 10 5 GeV. When both singlet and doublet VLFs are present, i.e. when a VLF family is added, the situation changes since a Yukawa coupling (ỹ) with the Higgs can be written down. Nevertheless, whenỹ is small, the behavior is similar to the above two cases. For a VLQ family with one singlet and one doublet VLQ added, as can be seen in Fig. 4, for a smallỹ = 0.1 and for M V L 10 5 GeV the EW minimum becomes absolutely stable. Thus, as we see in these examples, the presence of either SU(2) singlet VLQs, or doublet VLQs, or a full family with a small enoughỹ, allows the intriguing possibility that the EW vacuum is rendered absolutely stable.
For example, the hidden-sector dark matter model in Ref. [30] contains a singlet VLQ mediating loop-level couplings between the hidden-sector dark matter and the SM. Such models can also be written down with a doublet VLQ. For proper choices of the number of VLQs and masses, it is interesting that the Higgs vacuum could be absolutely stable in such models, unlike in the SM in which it is meta-stable.

Comparison with the analytical approximation of S B
Here we compare our numerical results for S B obtained in Sec. 3.2 with an analytical approximation developed in Refs. [31,32], which is, where t is a typical scale at which the bounce makes the transition from large field values to v. This approximation can yield a reasonably good estimate of S B when the bounce transition happens at a fairly constant value of λ(t), i.e. when h 0 is close to where β λ (h 0 ) ≈ 0. Furthermore, when S B is so large that errors due to the transition not happening at a constant λ(t) are small compared to S B , this approximation yields a good enough estimate. When these conditions are not realized, one has to be cautious in using the expression in Eq. (28). We elaborate on this statement below with many examples. In Fig. 11 in the left-panel we show β λ (µ) vs. µ where µ ≡ h(ρ). In the right-panel we show the (absolute value of) integrand of Eq. (26) made dimensionless by multiplying the integrand by 1/m 4 t and denoted as |Î(ρ)|, vs. λ(h(ρ)), with ρ being the parameter (not shown). As shown in the topmost panel in Fig. 11, for the SM, it is evident that most of the contribution to the integral comes from when λ takes a specific value. For the SM, we can compare the S B computed numerically in Sec. 3.2, which is 2866, with the S B got from the approximation in Eq. (28) with the t taken to be at the scale at which β λ = 0 where λ = −0.009, which gives S approx B = 2848. This is in excellent agreement with our numerical computation of S B , and as discussed earlier, this is because β λ = 0 does get satisfied for the SM, presenting a natural choice for t. That this approximation works is also borne out by the plot showing |Î(ρ)| for the SM, where most of the contribution to the S B integral is indeed coming for λ = −0.009, where the bounce spends most of its time (ρ). Indeed, Eq. (28) was put-forth for the SM, where it can be safely applied.
As we see from the last three panels from top in Fig. 11, with VLF present, β λ is not close to zero anywhere, and thus there is no clear choice of t that is suggested. In such a situation, we cannot use Eq. (28) but have to compute S B numerically. Indeed as the |Î(ρ)| for these cases show, the integral gets its contributions for a range of λ. These show the inadequacy of the approximate formula, and that a numerical evaluation is necessary. We have therefore computed the bounce EOM numerically and the bounce action for it, from which we computed the tunneling probability.

Conclusions
We study the stability of the electroweak vacuum in the presence of new vector-like fermions. We work with the one-loop renormalization group improved Higgs effective potential, identifying the Higgs field value h ≡ µ. We first review the computation of the beta-functions in the SM, paying particular attention to the SM fermion contributions. We use dimensional regularization for our computation. We then derive the VLF contributions to the 1-loop beta-functions which can be applied to various SU (3) and SU (2) representations, namely VLQs and VLLs. We apply this to a few example cases with singlet VLFs, doublet VLFs, and a family consisting of one doublet VLF and one singlet VLF coupled to the Higgs via the Yukawa couplingỹ. We numerically run the RGE to determine the scale at which λ(µ) becomes zero and goes negative. The Higgs effective quartic-coupling λ(µ) becoming negative signals that the EW vacuum is a false vacuum and is unstable, and can tunnel away quantum mechanically via barrier penetration to (large) field values that have a lower effective potential. We compute the probability P tunl that the EW vacuum would have tunneled away by a true-vacuum bubble nucleating in our Hubble four-volume. Computing P tunl requires computing the bounce configuration in Euclidean spacetime, and the value of the Euclidean action S B for the bounce configuration. We solve the bounce configuration EOM numerically and compute S B for it.
We compare our numerical evaluation with the approximation commonly used for S B , which is written in terms of λ at a single scale where β λ (µ) is approximately zero. This is because the bounce transition is mostly completed when λ(µ) has this value. For the SM, there is such a scale which is about 10 16 GeV, and we verify by comparing with our numerical evaluation that the approximation is perfectly adequate. When VLFs are present, there is no scale at which β λ (µ) is close to zero, and so the approximation cannot be applied. A numerical evaluation is then required, which we resort to.
We take example cases where a single VLL or VLQ family is added and show the bounce transition, compute S B for it and obtain P tunl . We find that P tunl is extremely sensitive toỹ as it exponentially depends on S B . Interestingly, we find that for some VLF representations and parameters, adding only singlet VLFs, a doublet VLFs, or a full family with a small enoughỹ, λ stays positive to arbitrarily large scales, i.e. the EW vacuum is rendered absolutely stable, unlike in the SM in which it is meta-stable. For other parameters, adding VLFs still keeps the EW vacuum meta-stable, either with a larger P tunl than in the SM, or a smaller P tunl .
In summary, our work here helps us get an idea of what the impact of VLFs is on the stability of the Higgs electroweak vacuum.
Here, we review the calculation of the fermion contributions to the 1-loop β-functions in the SM so that we can extend this to include VLF contributions in the next section. Since we are interested in field values h much larger than the particle masses, we neglect the particle masses.
We expand the SM Lagrangian density shown in Eq. (5) by writing the Higgs doublet as H = GeV is the EW VEV, h is the physical Higgs boson and φ i are the Goldstone bosons. The Lagrangian density in terms of the bare fields {h, t 0 } and bare coupling y 0 is where for notational brevity we denote y t just as y, and the covariant derivatives are in the usual notation. In terms of the renormalized fields and counter-terms, we have for the {t, h} sector where the renormalized fields h, t are defined byh = √ Z h h, t 0L,R = Z t L,R t L,R , and the renormalized Yukawa coupling y is defined by y 0 = Z y y. Expanding as a perturbation series in y, we define a's to leading order as (Z h − 1) ≡ a h y 2 /2, (Z t L,R − 1) ≡ a t L,R y 2 /2, (Z y − 1) ≡ a y y 2 /2, and we also define (Ẑ y − 1) ≡ (Z y Z h Z t L Z t R − 1) = (a y + a h /2 + a t L /2 + a t R /2)y 2 /2 ≡â y y 2 /2. Similarly the renormalized Lagrangian density for the Goldstone fields can be written down.
We give next the one-loop corrections involving the t, which we compute using dimensional regularization. We present only the SMF t contributions, since our goal is to use these to derive the VLF contributions to the β-functions later.
We find the Goldstone boson contributions proportional to y to be as follows: the φ 3 contribution to Σ t L,R is equal to the h contribution, the φ 1,2 contribution to Σ t R is each equal to the h contribution, and to Σ t L is proportional to y b which we drop and take to be zero, the φ 3 contribution to the ht LtR and ht RtL vertices V is each negative of the h contribution to the same, and the φ 1,2 contribution to V are proportional to y b and hence we take them to be zero.
One way to extract the β-function is from the divergent part of the bare coupling. 3 From the contributions computed above, we find the contributions proportional to y to be y 0 ⊃ y + (y 3 /(16π 2 ))((3 + 2N c )/2)(1/ ), from which we obtain the fermionic contribution to β y as after including the t, h, φ 1,2,3 contributions. This is in agreement with the results in Refs. [33,3], for example. Interestingly, the φ 1,2,3 contribute zero after including all their contributions. The yg 2 3 and yg 2 2 contributions to β y can be written as [3] β y ⊃ y which are included in Eq. (8). To derive the yg 2 1 contribution, we start by extracting the relevant Feynman rules for the hypercharge gauge boson B µ interactions. With all momenta going into the vertex, with Y L,R being the hypercharges of ψ L,R and Y H = 1/2 being the hypercharge of the Higgs doublet, we have the Feynman rules: Computing the B µ contribution at 1-loop order in the 't Hooft-Feynman ξ = 1 gauge, we obtain the following divergent pieces: the ψ RψL h vertex correction due to B µ exchange gives iV (Bµ) ⊃ −i8Y L Y R yg 2 /( √ 2 16π 2 ); the Higgs 2-point function correction due to φ 3 − B µ exchange gives iΣ ; and the fermion 2-point function corrections due to B µ exchange gives iΣ L,R p /(16π 2 ). We include in the counter-terms a piece to cancel these divergences at the subtraction scale p 0 , i.e., i(Z h − 1) . Thus, since the bare coupling is y 0 = Z y y, we get the contribution where we make use of Y H = Y R − Y L required for U (1) Y invariance of the Yukawa term in the Lagrangian, and g = 3/5g 1 . For the top, using Y L = 1/6, Y R = 2/3, we obtain the contribution shown in the last term of Eq. (8).
We compute the fermion loop contribution to the Higgs 4-point vertex that is proportional to y, from which we can write the bare coupling as λ 0 ⊃ λ − N c y 4 /(8π 2 ) (1/ + finite). From this, we infer that this contribution leads to Writing the Higgs four-point vertex as iλ eff , the contribution to its evolution, i.e. β λ , due to the fermion loop in the h-leg is just four times √ Σ h . Thus, for this contribution we have iλ eff ⊃ (−iλ)(i/p 2 )(4×iΣ φ /2), and from Eq. (31), we have λ(µ) = λ(M )+N c λy 2 /(4π 2 ) ln(µ/M )+..., where µ is the renormalization scale, and M is a subtraction scale. From this, and since β λ = dλ(µ)/d ln µ, we have We turn next to the β-functions of the gauge couplings g a = {g 3 , g 2 , g 1 }, focusing on the SM fermion contribution. We recall the definition β a = g 3 a b a /(16π 2 ). For β g 3 we have the well-known result (see for example Ref. [26]) where the second term is due to fermions, with n 3 as the number of colored fermions in the fundamental representation of SU (3). Note that the top-quark is vector-like with respect to the SU(3). In the SM, at large µ, we have n 3 = 6 for three generations of quarks, which implies b 3 = −7.
Similarly, for β g 2 we have where we have taken N = 2 for SU(2) in the first term, the second term is the fermion SU (2) doublet contribution with n 2 being the number of doublet fermions, and the last term is the Higgs doublet contribution. Since the SM fermions are chiral under SU (2) with only the L chirality contributing, we include an extra factor of 1/2 in the second term (since we are neglecting the effects of masses). Thus in the SM, at large µ, n 2 = (3N c + 3) for the three generations of quark and lepton doublets, which yields b 2 = −19/6. Lastly, for β g 1 we have where the sum in the first term is over all fermions f with hypercharge Y f , and the sum in the second term is over all complex scalars φ with hypercharge Y φ . We recall that we use SU (5) normalization for g 1 , i.e. the SM hypercharge gauge coupling g is related to g 1 by g 1 = 5/3g . Thus in the SM, at large µ, for three generations, f Y 2 f = 3 × (10/3) = 10, and for one Higgs doublet containing two complex fields with Y H = 1/2, φ Y 2 φ = 2(1/2) 2 = 1/2, we get b 1 = 41/10. This agrees with, for example, Ref. [34].
We complete our derivation of the SM fermion contributions to the 1-loop β-functions. After adding the other contributions, the complete 1-loop β-functions are as given in Eqs. (7)-(9).

B VLF contributions to the RGE
Here we extend the β-functions derived in Sec. 2.1 and App. A to include VLF contributions, which we denote as β V LF κ . We first derive the β V LF ga . The β V LF g 3 is got easily from Eq. (37). Since the SM quark is vectorlike with respect to SU (3), we have an identical contribution for a VLQ, and we obtain the result shown in Eq. (16). For obtaining β V LF g 2 , we note that this is similar to the SU (3) contribution owing to the fact that for a VLF SU (2) is also vector-like just as the SMF was for SU (3). Thus taking twice the second term in Eq. (38) will give us β V LF g 2 as given in Eq. (17). Since for a VLF both L and R chiralities contribute, we take twice the first term in Eq. (39) to obtain β V LF g 1 as given in Eq. (18); the 2n 2 is just the number of fermions in n 2 doublets having hypercharge Y χ .
Next, we derive the contributions present only when a full family is added, i.e. when theỹ operator of Eq. (13) can be written down.
Let us recall that the SM top (and bottom) sector have the following Feynman rules for the couplings with the Higgs-doublet fields {h, φ 1,2,3 } (all vertices have an overall (−iy t / √ 2)) written in terms of the Dirac spinors t = (t L t R ) T and b = (b L b R ) T : Now, when a full VLF family is present, we have the SU (2) doublet VLF χ = (χ 1 χ 2 ) T and a singlet ξ. We can assemble the following Dirac spinors: ψ 1 = (χ 1L ξ R ) T , ψ 2 = (ξ L χ 1R ) T , ξ = (ξ L ξ R ) T and χ 2 = (ξ 2L ξ 2R ) T . Using these, we can write the Feynman rules with the Higgs-doublet fields {h, φ 1,2,3 } as follows (all vertices have an overall (−iỹ/ √ 2)): We write it this way to bring forth the analogy between the SMF and VLF, with the realization that for the VLF we have two Dirac sets that each have a similarity with the SM couplings. The first Dirac fermion ψ 1 has identical couplings, while the second Dirac fermion ψ 2 has couplings that is similar but not identical, with a change i → −i in the φ 3 couplings, and P L → P R in the φ 1,2 couplings. We observe that all the diagrams that contribute to the β-functions are immune to both of these changes, and therefore each of them give the same SM contribution as for the t quark. Furthermore, the Goldstone bosons with each Dirac fermion contributes zero to the β-function as in the SM. Thus, to obtain the VLF contributions to β V LF λ , we just multiply the SM contribution after combining Eqs. (35) and (36) by a factor of 2 and obtain Eq. (19). Next, the VLF contribution to β yt is due to only the wavefunction renormalization contribution to h, i.e. Σ h , and this contribution can be got from the second term in Eq. (32), but multiplied by 2 since there are two VLF Dirac sets as argued above, and changing the coupling to y tỹ 2 instead of y 3 t , which then gives us the β V LF yt in Eq. (20). Next, consider the evolution of either the hψ 1ψ1 coupling, or the hψ 2ψ2 coupling, either of which isỹ. The VLF contribution to βỹ is due to these three contributions: (i) the vertex contribution proportional to 3ỹ 3 as in the first term in Eq. (32), (ii) the VLF contributions in Σ h which yields twice the second term in Eq. (32) proportional to 2 × 2N cỹ 3 since each of the ψ 1 and ψ 2 contribute as in the SM, and (iii) the top-quark contribution in Σ h which yields 2N c y 2 tỹ as in the second term in Eq. (32). Adding these three contributions then gives the first part of βỹ in Eq. (21). We write theỹg 2 a contributions to βỹ following Eqs. (33) and (34), which gives the last part in Eq. (21). We thus complete the derivation of the VLF contributions to the β-functions given in Eqs. (16)- (21).