Higgs Vacuum Stability with Vector-like Fermions

TL;DR

This paper investigates how vector-like fermions influence the stability of the Higgs vacuum, showing that certain parameters can make the vacuum absolutely stable up to the Planck scale, extending the standard model's predictions.

Contribution

It provides a detailed calculation of the VLF effects on Higgs vacuum stability, including one- and two-loop beta-functions, and analyzes the conditions for absolute stability or metastability.

Findings

VLFs can stabilize the Higgs vacuum up to the Planck scale.

The paper computes tunneling probabilities for metastable vacua.

It compares numerical and analytical bounce action calculations.

Abstract

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. We also include the significant two-loop 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.

Equations

$$\mathcal{V}=\frac{m_h^2}{2}{h}^2+\frac{\lambda }{4}{h}^4.$$

$$V_{\mathrm{eff}}(h)=\frac{m_{h\mathrm{eff}}^2}{2}{h}^2+\frac{{\lambda }_{\mathrm{eff}}(h)}{4}{h}^4,$$

$$V_{\mathrm{eff}}(h)=\frac{{\lambda }_{\mathrm{eff}}(h)}{4}{h}^4.$$

$$\frac{d\lambda (\mu )}{d\ln \mu }={\beta }_{\lambda }\left(\lambda (\mu ),y_t(\mu ),g_3(\mu ),g_2(\mu ),g_1(\mu ),...\right).$$

$$\mathcal{L}\supset \bar{t}ito0.0ptto5.30904pt\text{hfil}/\text{hfil}\text{hss}\partial t-\lambda (H^{†}H{)}^2-(y_t{\bar{q}}_L\cdot H^{\ast }{t}_R+h.c.),$$

$$\frac{d{\kappa }_i(\mu )}{d\ln \mu }={\beta }_{{\kappa }_i}({\kappa }_j(\mu )).$$

$${\beta }_{\lambda }^{(1)}$$

$${\beta }_{y_t}^{(1)}$$

$${\beta }_{g_a}^{(1)}$$

$${\beta }_{\lambda }^{(2)}$$

$${\beta }_{y_t}^{(2)}$$

$${\beta }_{g_3}^{(2)}$$

$$\mathcal{L}\supset -M_{\chi }\bar{\chi }\chi -M_{\xi }\bar{\xi }\xi -(\tilde{y}\bar{\chi }\cdot H^{\ast }\xi +h.c.),$$

$$\mathcal{L}\supset -\frac{\tilde{y}}{\sqrt{2}}h({\bar{\chi }}_1\xi +\bar{\xi }{\chi }_1).$$

$${\beta }_{g_3}^{(1)\mathrm{VLF}}$$

$${\beta }_{g_2}^{(1)\mathrm{VLF}}$$

$${\beta }_{g_1}^{(1)\mathrm{VLF}}$$

$${\beta }_{\lambda }^{(1)\mathrm{VLF}}$$

$${\beta }_{y_t}^{(1)\mathrm{VLF}}$$

$${\beta }_{\tilde{y}}^{(1)}$$

$${\beta }_{g_3}^{(2)\mathrm{VLF}}$$

$${\beta }_{\lambda }^{(2)\mathrm{VLF}}$$

$${\beta }_{y_t}^{(2)\mathrm{VLF}}$$

$${\beta }_{\tilde{y}}^{(2)}$$

$$\mathcal{L}=\frac{1}{2}{\partial }^{\mu }h{\partial }_{\mu }h-V_{\mathrm{eff}}(h);S=\int d^{4}x\mathcal{L},$$

$$S_E=\int d^4\rho \left[\frac{1}{2}({\partial }_{i}h{)}^2+V_{\mathrm{eff}}(h)\right],$$

$${\partial }_i^{2}h=\frac{\partial V_{\mathrm{eff}}}{\partial h}.$$

$$\frac{d^{2}h}{d{\rho }^2}+\frac{3}{\rho }\frac{dh}{d\rho }=\frac{\partial V_{\mathrm{eff}}}{\partial h},$$

$$S_B=2{\pi }^2{\int }_0^{\infty }d\rho {\rho }^3\left[\frac{1}{2}{\left(\frac{d{h}_B}{d\rho }\right)}^2+V_{\mathrm{eff}}(h_B)\right].$$

$$P_{\mathrm{tunl}}=(h_0/m_t{)}^4{e}^{(404-S_B)}.$$

$$S_B^{\mathrm{approx}}=\frac{8{\pi }^2}{3(-\lambda (t))},$$

$$\mathcal{L}$$

$$i{\Sigma }_h(p_h)+i{a}_h\frac{y^2}{2}{p}_h^2=-i\frac{N_c{y}^2}{8{\pi }^2}{p}_h^2\ln \left(\frac{p_h}{p_0}\right).$$

$$iV(p_h)=i{y}^3/(16\sqrt{2}{\pi }^2)\left(1/ϵ-\ln p_h/2+\ln \sqrt{4\pi }-\gamma /2+1/2+\mathcal{O}(ϵ)\right),$$

$$a_y=(2N_c+3)/(16{\pi }^2)\left(1/ϵ-\ln p_0/2+\ln \sqrt{4\pi }-\gamma /2+1-\ln 2/(2N_c+3)\right).$$

Full text

1 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, 5]) 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_{\rm eff}(h)$ has a smaller value for $h\sim 10^{10}~{}$GeV when compared to its value at the EW vacuum expectation value (VEV) $v\approx 246$ GeV, i.e., $V_{\rm eff}(h\!\sim\!10^{10}~{}GeV)<V_{\rm 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. [6, 7, 8, 9], Higgs-portal dark matter models in Refs. [10, 11, 12, 13, 14], gauge-coupling unification in Refs. [15, 16, 17, 18, 19], neutrino mass generation and vacuum stability in Refs. [20, 21, 22, 23, 24, 25], universal extra dimension model in Ref. [26], SM extensions with an additional $U(1)$ gauge symmetry in Refs. [27, 28, 29, 30], models with extended scalar sector in Refs. [31, 32], a combination of these in Ref. [33, 34, 35], models of inflation in Refs. [36, 37], and effective models in Refs. [38, 39, 40]. Motivated by these considerations, we study the effect of VLFs that are coupled to the Higgs on EW vacuum stability.111 New chiral (4th generation) fermions that get their mass from the Higgs are severely disfavored with a single Higgs doublet from the recent LHC Higgs cross-section and couplings measurements. In contrast, VLFs tend to have milder constraints on them owing to their nice decoupling property.

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

[TABLE]

Including quantum effects, we can write the quantum Higgs effective potential as

[TABLE]

where $\lambda_{\rm eff}$ has a dependence on $h$ of the form $\ln{(h/M)}$ with $M$ being a subtraction scale. For $h\gg m_{h}$ (the physical Higgs mass $m_{h}\approx 125$ GeV), the mass term has a negligible effect, and thus, to an excellent approximation, we can write

[TABLE]

Denoting the field value as $h\equiv\mu$ and denoting $\lambda_{\rm eff}(\mu)$ as just $\lambda(\mu)$, it can be shown (see for example Ref. [41]) that $\lambda(\mu)$ obeys a renormalization group equation (RGE) of the form

[TABLE]

The RGE is interpreted now as an evolution with field value $h=\mu$, and the $\beta$-function $\beta_{\lambda}$ is the usual $\beta$-function for the coupling $\lambda$, governed by the RGE. The $\lambda(\mu)$ obtained by integrating the RGE has the leading logs of the form $\log^{n}\!{(\mu/M)}$ resummed. $\beta_{\lambda}$ is shown as a function of $\lambda$ itself, and also of the other couplings that contribute significantly, which, in the SM, are the top Yukawa coupling $y_{t}$ and the $SU(3),SU(2),U(1)$ gauge couplings $g_{a}=\{g_{3},g_{2},g_{1}\}$. All these couplings also evolve with $\mu$ via analogous RGE equations with their corresponding $\beta$-functions $\beta_{y_{t}},\beta_{g_{a}}$. We neglect the contributions of the other SM couplings to the $\beta$-functions as they contribute insignificantly. From Eq. (3) we see that for $h\gg m_{h}$, the instability is signalled by the Higgs quartic effective coupling $\lambda(\mu)$ becoming negative.

As we show explicitly later, $\beta_{\lambda}$ obtains a negative contribution from $y_{t}$, while it obtains a positive contribution from $\lambda(\mu)$ itself and from gauge couplings. Thus, the top quark has the important effect of decreasing $\lambda(\mu)$, and for $y_{t}$ as large as in the SM, for the observed $m_{h}$, it drives $\lambda(\mu)$ negative at higher energies, signaling vacuum instability. The effect of fermions coupled to the Higgs is generally to destabilize the electroweak vacuum, although in this work we show that this statement is not so definite. Many extensions beyond the SM (BSM) include new fermions, and the question we address in this work is what the effects of new fermions might be on Higgs vacuum stability in light of the observation made above.

The subject of this work is to include VLF contributions to the $\beta$-functions and find the consequences for EW vacuum stability and how it is changed from the SM. We ask if $\lambda(\mu)$ still becomes zero with VLF present, and if so at what $\mu$, and compare it with the SM case. If the vacuum is unstable, we compute the tunneling probability to ascertain if it decays within the age of the Universe, in which case it is unacceptable. On the contrary, if the lifetime in the EW vacuum is comparatively much larger than the age of the Universe, it is metastable and phenomenologically acceptable. To this end, we study some simple VLF extensions of the SM, where the VLFs are either in the trivial or fundamental representations of $SU(3),SU(2)$ and $U(1)$, and demonstrate their effects on Higgs vacuum stability. In particular, the VLFs we add are of two kinds, namely, $SU(3)$ triplet vector-like quarks (VLQ) and $SU(3)$ singlet vector-like leptons (VLL).

The paper is organized as follows: In Sec. 2 we list the 1-loop RGE in the SM and include some significant 2-loop corrections from the literature. We present a derivation of the fermionic contributions to the RGE in Appendix A. We then derive the 1-loop VLF contributions to the RGE and add these to the SM RGE. We present the calculational details of the VLF contributions in Appendix B. We add the significant two-loop contributions to the beta-functions. We integrate the RGE numerically and show the evolution of the couplings as a function of the field value $h\equiv\mu$. In Sec. 3 we compute the probability that our electroweak vacuum would have tunneled into a deeper true vacuum in our Hubble volume in the case where the EW vacuum is metastable. We do so by solving for the bounce configuration numerically and computing the Euclidean action for this. In Sec. 4 we make some remarks for the case when the VLFs render the EW vacuum absolutely stable. In Sec. 5 we compare our numerical evaluation of the bounce action to an approximation commonly used in the literature, and provide a cautionary note on when the approximation can be applied. We offer our conclusions in Sec. 6.

2 Renormalization group improved Higgs effective potential

We have in the SM the Lagrangian density, showing only the terms relevant to our analysis here,

[TABLE]

where the $\cdot$ 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 $\beta$-functions and extend them to include the VLF contributions.

2.1 SM RGE

We first discuss the SM RGE $\beta$-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 $\lambda(\mu)$ becomes zero and compare our results with those in the literature. Denoting the relevant SM couplings generically as $\kappa_{i}=\{\lambda,y_{t},g_{3},g_{2},g_{1}\}$, the RGE are of the form

[TABLE]

We derive the fermion contributions to $\beta_{\kappa}$ 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. [4]). Putting these together, the 1-loop $\beta$-functions, $\beta^{(1)}_{\kappa}$, are

[TABLE]

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^{\prime}$ is related to $g_{1}$ by $g_{1}=\sqrt{5/3}g^{\prime}$.

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 $\beta$-functions that depend on $y_{t}$ and $g_{3}$ as they are numerically the most significant. They are (see for example Ref. [4])

[TABLE]

where $n_{3}^{\rm(SM)}=6$ is the number of $SU(3)$-triplets (i.e. quarks) in the SM.

We use these RGE to determine the Higgs field value $\mu$ at which $\lambda(\mu)$ becomes zero, signalling vacuum instability. This will be discussed in Sec. 2.3. We discuss next the VLF contributions to the RGE.

2.2 VLF contributions to the RGE

We add an $SU(2)$ doublet VLF $\chi=(\chi_{1}\ \chi_{2})^{T}$ and an $SU(2)$ singlet VLF $\xi$ and couple it to the Higgs as follows

[TABLE]

where the $\cdot$ represents the antisymmetric combination in $SU(2)$ space.222If another $SU(2)$ singlet VLF $\zeta$ is added, we can add the terms ${\cal L}\supset-M_{\zeta}\bar{\zeta}\zeta-(\tilde{y}_{2}\bar{\chi}H\zeta+h.c.)$. After adding the $\zeta$, 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 $\zeta$ in our work here, and therefore will not include the $\tilde{y}_{2}$ term.

Extracting the Higgs interactions from this yields

[TABLE]

If $\chi$ and $\xi$ have color $N_{c}^{\prime}=3$ we call them vector-like quarks (VLQ) and if they are trivial under $SU(3)$, i.e. $N_{c}^{\prime}=1$, we call them vector-like leptons (VLL). The $SU(3),SU(2)$ and $U(1)$ gauge interactions are standard and we do not show them explicitly. We denote the hypercharge of $\chi$ as $Y_{\chi}$, that of $\xi$ as $Y_{\xi}$, and that of the Higgs doublet is $Y_{H}=1/2$ as in the SM. If $N_{c}^{\prime}=3$ the VLQ have gluon interactions, while if $N_{c}^{\prime}=1$ the VLL do not have gluon interactions.

For SM-like choices of $Y_{\chi}$ and $Y_{\xi}$, mixed Yukawa couplings between the VLF and the standard model fermions (SMF) can be written down. However, collider, flavor changing neutral current and other precision constraints restrict how large such couplings can be (for details, see for example Ref. [42]). 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 $\tilde{y}$, are

[TABLE]

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 to the Higgs (a family is a doublet and a singlet both present), and $\hat{n}_{F}^{VLQ}=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 $\tilde{y}$ can be written down with the Higgs. Only VLQs contribute to $\beta_{g_{3}}$ and VLLs do not. For instance, for one VLQ family of $\chi$ and $\xi$, we have $N_{c}^{\prime}=3$, $n_{3}=3$, $n_{2}=1$, $n_{1}=1$, $n_{F}=1$ and $\hat{n}_{F}^{VLQ}=1$.

To improve precision, we include the dominant 2-loop VLF contributions to the $\beta$-functions obtained from the package ‘SARAH’ [43, 44], which are,

[TABLE]

where, $n^{VLQ}_{F}$ is the number of colored families, and as noted earlier, $n_{3}^{\rm(SM)}=6$. We have explicitly checked that the above dominant contributions closely reproduce numerically the full 2-loop running from SARAH.

Our goal in this work is to analyze the stability of the electroweak (EW) vacuum for which the behavior of $V_{\rm 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/\mu)$ for $\mu\gg m$ where $m$ collectively denotes the particle masses. We include the VLF contributions only for $\mu\geq M_{VL}$, where $M_{VL}$ is the vector-like fermion mass.

2.3 RGE Numerical integration results

We take the input parameters as follows and as compiled in Ref. [4], with the renormalization point taken as the top mass scale $\tilde{m}_{t}$:

• The EW VEV: $v=246.2~{}$GeV,
• The Higgs quartic: $\tilde{\lambda}=0.12710$ (NNLO),
• The top Yukawa coupling: $\tilde{y}_{t}=0.93558$ (partial 3-loop),
• The $SU(3)_{c}$ coupling constant: $\tilde{g}_{3}=1.1666$ (partial 4-loop),
• The $SU(2)_{L}$ coupling constant: $\tilde{g}_{2}=0.64755$ (NLO),
• The $U(1)$ coupling constant: $\tilde{g}_{1}=\sqrt{5/3}\,g^{\prime}=\sqrt{5/3}\times 0.35937$ (NLO).

In terms of these inputs, we set the top mass to be $\tilde{m}_{t}=\tilde{y}_{t}v/\sqrt{2}$ and the Higgs mass $\tilde{m}_{h}=\sqrt{2\tilde{\lambda}}\,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 $\lambda(\mu)$, $y_{t}(\mu)$, $g_{3}(\mu)$, $g_{2}(\mu)$, $g_{1}(\mu)$. We take the inputs given above at $\tilde{m}_{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 $\mu\geq M_{VL}$.

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\equiv\mu$. The Higgs quartic coupling $\lambda$ becomes zero and goes negative at about $\mu\sim 10^{10.5}~{}$GeV, while all the other couplings stay positive all the way up to $M_{\rm Pl}$. It is interesting that $\beta_{\lambda}$ approaches zero for large field values (cf Fig. 9). 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\equiv\mu$ for $(n)$ degenerate $SU(2)$ singlet VLQs for various $M_{VL}$, where the $M_{VL}$ values are shown in the notation $(rEm)\equiv r\times 10^{m}~{}$GeV.

When 3 or more degenerate singlet VLQs of mass 3 TeV are added, interestingly, $\lambda$ 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 $\beta_{g_{3}}$, and also to $\beta_{g_{1}}$ if the VLQ has hypercharge, and because of the coupled nature of the RGEs, $\lambda(h)$ does see the effect of the VLQ. In particular, even if $\tilde{y}$ is very small, the VLQ contribution to $\beta_{g_{3}}$ given in Eq. (15) still remains, which, being positive, results in $g_{3}(\mu)$ being larger for larger $\mu$ as compared to the SM. A larger $g_{3}$ means that the second term in Eq. (8) is more negative, causing the $y_{t}(\mu)$ to be smaller in comparison to the SM case. A smaller $y_{t}$ implies a less negative contribution to $\beta_{\lambda}$ from the third term of Eq. (7), which means that the $\lambda(\mu)$ is larger with VLQ present. For large enough $n_{3}$ this results even in a turn-around to a positive $\beta_{\lambda}$ allowing for the possibility of $\lambda$ never going negative. 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\equiv\mu$ for $(n)$ degenerate $SU(2)$ doublet VLQs for various $M_{VL}$.

We observe that when we add one or more doublet VLQs of mass 3 TeV, $\lambda$ 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. We discuss in Sec. 4 the implications to vacuum stability of $\lambda$ remaining positive. If we add five or more doublets with 3 TeV mass, we find that due to the large positive VLF contribution to $\beta_{g_{2}}$ given in Eq. (16), $g_{2}$ grows and becomes non-perturbative at $\mu\approx 10^{16}~{}$GeV, invalidating this perturbative analysis at around that scale. In Fig. 2, we have restricted the five doublet curves to the region $g_{2}<10$ so that our perturbative analysis is reliable. The negative contribution proportional to $g_{2}^{2}$ in $\beta_{y_{t}}$ given in Eq. (8) becomes significant as $g_{2}$ becomes large and leads to a smaller $y_{t}$. Also, $\beta_{g_{3}}$ gets a large positive VLF contribution from Eq. (15) causing $g_{3}$ to increase with $\mu$.

In Fig. 3 we show the evolution of the couplings with field value $h\equiv\mu$ for a degenerate family of one $SU(2)$ doublet VLL and one singlet VLL for various $M_{VL}$ and $\tilde{y}$. The VLL Yukawa coupling values are shown as $(\tilde{y})$ and the $M_{VL}$ values are shown in the notation $(rEm)\equiv r\times 10^{m}~{}$GeV.

For $\tilde{y}(M_{VL})\!=\!1$, we see that $\tilde{y}$ increases as $\mu$ increases. $y_{t}$ eventually starts increasing at large $\mu$, which is a behavior unlike in the SM. We notice that the scale at which $\lambda$ becomes negative decreases as $\tilde{y}$ increases, or as $M_{VL}$ decreases.

In Fig. 4 we show the evolution of the couplings with field value $h\equiv\mu$ for a degenerate family of one $SU(2)$ doublet VLQ and one singlet VLQ, for various $M_{VL}$ and $\tilde{y}$.

We see that for $\tilde{y}(M_{VL})=0.5$, $\tilde{y}$ decreases as $\mu$ increases. For $M_{VL}=3~{}$TeV, if $\tilde{y}>0.35$, the scale at which $\lambda$ becomes negative decreases as $\tilde{y}$ increases or as $M_{VL}$ decreases and is lesser than in the SM, while if $\tilde{y}<0.35$, the scale at which $\lambda$ becomes negative is larger than in the SM. In fact, for $\tilde{y}<0.3$, $\lambda$ stays positive all the way up to $M_{Pl}$. We see that when $\tilde{y}=0.1$ for example, $\lambda$ stays positive all the way up to $M_{Pl}$ for $M_{VL}$ 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 $\lambda$ 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_{\rm eff}<0$. We turn next to an analysis of this possibility and a computation of the tunneling probability.

3 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. [45]). 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}$.

3.1 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, 45] 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 ${\cal L}$ and the action $S$ for the Higgs field in Minkowski coordinates are

[TABLE]

with the effective potential defined in Eq. (2). We define the Euclidean time $\tau=it$, the Euclidean coordinates $\rho^{i}=\{\tau,x^{1},x^{2},x^{3}\}$ with $i=\{1,2,3,4\}$, and the invariant $\rho^{2}=\tau^{2}+(x^{1})^{2}+(x^{2})^{2}+(x^{3})^{2}$, and write the Euclidean action as

[TABLE]

where $\partial_{i}$ is with respect to the Euclidean coordinates $\rho^{i}$.

As mentioned earlier, taking $V_{\rm 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_{\rm eff}(h)<0$ for some $h$ typically much larger than $v$. If there exists field values for which $V_{\rm 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\gg v$, but is in the EW vacuum today. If for some large field value, $h\equiv\sigma$ say, suppose $V_{\rm eff}(\sigma)=0$, and suppose $V_{\rm eff}$ is negative for $h\gtrsim\sigma$. ($V_{\rm eff}$ may turn around and have a second minimum (or not) for $h>\sigma$ depending on other BSM contributions in the RGE.) Equivalently, from our definition of the field dependent coupling in Eq. (3), suppose $\lambda(\sigma)=0$, and that $\lambda(h)$ is negative for $h\gtrsim\sigma$. The vacuum configuration defined to have total energy $E=0$ at $h=v$ can quantum mechanically tunnel to $h\geq\sigma$ with $V_{\rm eff}(h)<0$. If the vacuum were to tunnel so, the field then runs down the potential classically toward large field values $h>\sigma$. The tunneling probability is given in terms of the bounce configuration [45] which satisfies $\delta S=0$, starting with $h(t\!=\!-\infty)=v$, attaining a value $h(t\!=\!0)\equiv h_{0}\geq\sigma$ and returning to $h(t\!=\!\infty)=v$. This configuration is a solution of the equation of motion (EOM). In Euclidean coordinates, the EOM reads

[TABLE]

We look for an $O(4)$ symmetric solution [46], which implies that it depends on $\rho$, i.e. $h(\rho^{i})=h(\rho)$. The EOM then reads

[TABLE]

with the boundary conditions (BC) $(dh/d\rho)(\rho\!=\!0)=0$ and $h(\rho\!\rightarrow\!\infty)=v$. We must also have $h(\rho\!=\!0)=h_{0}\geq\sigma$ for this to represent tunneling. This EOM is identical to that of a classical particle moving in a potential $-V_{\rm eff}$ with a “friction” term present (second term in Eq. (28)) that dies-off as $1/\rho$ as $\rho$ increases.

In Euclidean space-time, the bounce configuration $h_{B}(\rho)$, will have the feature of a fairly sharp transition in $\rho$ 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 [1] $\Delta P_{\rm tunl}/\Delta V_{4}=M^{4}\exp{(-S_{B})}$, where we have included a prefactor of $M^{4}$ on dimensional grounds with $M$ an appropriate mass scale, $\Delta V_{4}$ is a unit 4-space-time volume, and $S_{B}$ is the Euclidean action for the bounce configuration $h_{B}(\rho)$ given by

[TABLE]

We make the choice $M^{4}=V_{\rm 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}\sim(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_{\rm tunl}\sim(dP_{\rm tunl}/dV_{4})V_{4}$, which gives [1]

[TABLE]

If $P_{\rm tunl}\ll 1$ for the given $V_{\rm eff}$, we deem this as an acceptable situation. In other words, if $P_{\rm tunl}\gtrsim 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_{\rm eff}$ we consider is disfavored. Evidently, the larger $S_{B}$ is, the smaller $P_{\rm 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}(\rho)$, compute $S_{B}$ for this $h_{B}$, and compute $P_{\rm tunl}$. We do this for the SM and some VLF extensions.

3.2 Tunneling probability numerical results

Here we describe the method we use to solve the bounce EOM numerically, obtain the bounce configuration $h_{B}(\rho)$, the bounce action $S_{B}$ and the tunneling probability $P_{\rm tunl}$ for the SM and various VLF extensions.

The value of $S_{B}$ largely depends on the behavior of $V_{\rm eff}$ at large field values $h$ where the $h^{4}$ term dominates, which is why we included only that term in $V_{\rm eff}$. Nevertheless, for completeness, we insert an EW minimum at $v$ by including it in $\lambda(h)$ for $h\sim v$ as follows. Keeping the $h^{2}$ term, we have for $h\sim v$ the potential $V=-(\mu_{h}^{2}/2)h^{2}+(\lambda/4)h^{4}$, which we define to be $V_{\rm eff}=(\lambda(h)/4)h^{4}$ as in Eq. (3). This definition implies that we have the effective quartic coupling given by $\lambda(h)=\lambda(v)(1-2v^{2}/h^{2})$ for $h\sim v$. Slightly above the scale $v$, we match this to the $\lambda(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_{\rm eff}$ and make $V_{\rm eff}(v)=0$.

We obtain a solution of the EOM in Eq. (28) numerically, subject to the BC $h(\rho\!\!=\!\!0)=h_{0}$, $(dh/d\rho)(\rho\!\!=\!\!0)=0$. $h_{0}$ is unknown and so we iteratively search for that $h_{0}$ that will lead to $h(\rho_{\rm end})=v$ and $(dh/d\rho)(\rho_{\rm end})=0$. Although in theory $\rho_{\rm end}\rightarrow\infty$, in practice it can be picked finite but large enough that the bounce has completed the transition from $h_{0}$ to $v$. In our numerical implementation, we work with the dimensionless quantities $\hat{\rho}\equiv\tilde{m}_{t}\rho$, $\hat{h}\equiv h/\tilde{m}_{t}$ and $\hat{V}_{\rm eff}(h)=V_{\rm eff}(h)/{\tilde{m}_{t}}^{4}$.

The “friction” term that goes like $1/\rho$ in Eq. (28) will be problematic numerically near $\rho\rightarrow 0$, and we therefore obtain an analytical solution in this region, valid for $\hat{\rho}\in(0,\epsilon)$ for $\epsilon\ll 1$, and match this onto a numerical solution of the EOM for $\hat{\rho}\geq\epsilon$. We now give the solution valid in $\hat{\rho}\in(0,\epsilon)$. For small $\hat{\rho}$, we expand as $(3/\hat{\rho})d\hat{h}/d\hat{\rho}\equiv s(\hat{\rho})=s_{0}+s_{1}\hat{\rho}+(s_{2}/2)\hat{\rho}^{2}+{\cal O}(\hat{\rho}^{3})$, and require all the $s_{i}$ to be finite so that the friction term is finite as $\hat{\rho}\rightarrow 0$. Integrating this, we find $\hat{h}(\hat{\rho})=\hat{h}_{0}+(s_{0}/6)\hat{\rho}^{2}+{\cal O}(\hat{\rho}^{3})$. Differentiating the earlier equation, we have $d^{2}\hat{h}/d\hat{\rho}^{2}=s_{0}+(4s_{1}/3)\hat{\rho}+(5s_{2}/6)\hat{\rho}^{2}+{\cal O}(\hat{\rho}^{3})$. We find $\lambda(\hat{h}(\hat{\rho}))=\lambda_{0}+({\beta_{\lambda}}_{0}/\hat{h}_{0})(d^{2}\hat{h}/d\hat{\rho}^{2})_{0}\hat{\rho}^{2}/2+{\cal O}(\hat{\rho}^{3})$, and $\beta_{\lambda}(\hat{h}(\hat{\rho}))={\beta_{\lambda}}_{0}+({\partial\beta_{\lambda}}/\partial\hat{h})_{0}(d^{2}\hat{h}/d\hat{\rho}^{2})_{0}\hat{\rho}^{2}/2+{\cal O}(\hat{\rho}^{3})$, where $\lambda_{0}\equiv\lambda(h_{0})$, ${\beta_{\lambda}}_{0}\equiv\beta_{\lambda}(h_{0})$, $(d^{2}\hat{h}/d\hat{\rho}^{2})_{0}\equiv(d^{2}\hat{h}/d\hat{\rho}^{2})(\hat{\rho}=0)$ and $({\partial\beta_{\lambda}}/\partial\hat{h})_{0}\equiv({\partial\beta_{\lambda}}/\partial\hat{h})({\hat{h}=\hat{h}_{0}})$. Also, $\partial\hat{V}_{\rm eff}/\partial\hat{h}=[\lambda(\hat{h})+\beta_{\lambda}(\hat{h})/4]\hat{h}^{3}$. Substituting these into the EOM in Eq. (28), we get by matching powers of $\hat{\rho}$, $s_{0}=(\lambda_{0}+{\beta_{\lambda}}_{0}/4)\hat{h}_{0}^{3}/2$, $s_{1}=0$ and $s_{2}=(3/8)(d^{2}\hat{h}/d\hat{\rho}^{2})_{0}[3(\lambda_{0}+{\beta_{\lambda}}_{0}/4)\hat{h}_{0}+({\beta_{\lambda}}_{0}/\hat{h}_{0}+({\partial\beta_{\lambda}}/\partial\hat{h})_{0}/4)\hat{h}_{0}^{3}]$. This is the solution valid for $\hat{\rho}\in(0,\epsilon)$, and the solution at $\hat{\rho}=\epsilon$ is got by substituting $\hat{\rho}=\epsilon$ in this.

Taking the $\hat{h}(\epsilon)$ and $(d\hat{h}/d\hat{\rho})|_{\epsilon}$ obtained as above at the point $\epsilon$ as a BC, we numerically integrate the EOM in Eq. (28) for $\hat{\rho}\in(\epsilon,\hat{\rho}_{\rm end})$ and obtain $\hat{h}_{B}(\hat{\rho})$ 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 $\hat{h}_{0}$. By an iterative search algorithm we are able to obtain the bounce configuration $h_{B}(\hat{\rho})$ using Mathematica.

Piecing together the analytical solution above and the numerical solution, we obtain the bounce configuration over the complete domain $\hat{\rho}\in(0,\hat{\rho}_{\rm 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 $V_{\rm eff}(h\!\equiv\!\mu)$ and the bounce configuration $h_{B}(\rho)$ obtained numerically are shown in Fig. 5.

The $V_{\rm eff}(\mu)$ is positive for smaller $\mu$, crosses zero at about $\mu\approx 10^{10.75}~{}$GeV, and is negative for larger $\mu$. 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. (29) that the value of the Euclidean bounce action is $S_{B}=2866$ (in $\hbar=1$ units). From this, we compute the tunneling probability into the true vacuum in our Hubble volume from Eq. (30) to be $P_{\rm tunl}\sim 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. [4].

Next, we solve the bounce EOM, and compute the $S_{B}$ and $P_{\rm 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, for various common mass $M_{VL}$ and various $\tilde{y}$.

For a VLL family with $M_{VL}=10^{3}~{}$GeV and $\tilde{y}=0.6$, the $V_{\rm eff}(h\!\equiv\!\mu)$, and the bounce configuration are shown in Fig. 6 (top row).

The $V_{\rm eff}(\mu)$ is positive for smaller $\mu$, crosses zero at about $\mu\approx 10^{6.5}~{}$GeV, and is negative for larger $\mu$. 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}=472$ and $P_{\rm tunl}\sim 10^{-6}$. 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_{VL}=10^{3}~{}$GeV and $\tilde{y}=0.61$, for which we find $S_{B}=422$ and $P_{\rm tunl}\sim 10^{17}$. 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_{\rm tunl}$ is extremely sensitive to $\tilde{y}$, with a small change of $0.01$ in $\tilde{y}$ between the two cases resulting in a change of $S_{B}$ of 50, which in turn results in a $P_{\rm tunl}$ 23 orders of magnitude different because of its exponential dependence on $S_{B}$. As another example, for a VLL family with $M_{VL}=10^{5}~{}$GeV and $\tilde{y}=0.57$, the bounce configuration is shown in Fig. 6 (middle row). For this bounce configuration, we find $S_{B}=498$ and $P_{\rm tunl}\sim 10^{-7}$ which is acceptable. For a VLL family with $M_{VL}=10^{7}~{}$GeV and $\tilde{y}=0.6$, the bounce configuration is shown in Fig. 6 (bottom row). For this bounce configuration, we find $S_{B}=500$ and $P_{\rm tunl}\sim 10^{-4}$ which is acceptable.

The hidden sector Higgs-portal dark matter model of Ref. [11] essentially behaves like a VLL family considered above, for the following reason. Although in the model of Ref. [11], the VLF dark matter is a singlet and does not couple directly to the Higgs, due to the Higgs mixing with a hidden-sector scalar, a coupling with the Higgs is induced with size $\tilde{y}\equiv\kappa 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. [11], we require $\tilde{y}\ll 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_{\rm 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_{VL}$ and various $\tilde{y}$. With the addition of a VLQ family, in Fig. 7 we show the regions of stability, meta-stability and instability as a function of $M_{VL}$ (in GeV) and $\tilde{y}$. In the region marked “stable” the Higgs electroweak minimum is the absolute minimum and is discussed further in Section 4; in the region marked “metastable” there is a lower minimum at large field values with $P_{\rm tunl}\lesssim{\cal O}(1)$, and in the region marked “unstable” $P_{\rm tunl}\gtrsim{\cal O}(1)$.

We find that for $\tilde{y}\gtrsim 0.5$ the $P_{\rm tunl}\gtrsim{\cal O}(1)$ quite independently of $M_{VL}$. This parameter-space leads to an unstable vacuum and we consider this region disfavored from the vacuum stability point of view.

3.2.1 Second minimum in $V_{eff}$

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. 8 for the case of a VLQ family with $M_{VL}=3\times 10^{3}~{}$GeV and $\tilde{y}=0.75$. The $V_{\rm eff}(h\!\equiv\!\mu)$, and the bounce configuration for this modified potential are shown in Fig. 8.

The $V_{\rm eff}(\mu)$ is positive for smaller $\mu$, crosses zero at about $\mu\approx 10^{4.65}~{}$GeV and becomes negative, obtains a minimum at about $\mu\approx 10^{5}~{}$GeV, crosses zero again and becomes positive for larger $\mu$. 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}=3071$ and $P_{\rm tunl}\sim 10^{-1150}$, 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. $\rho$) near the minimum since $d\phi/d\rho$ is small there, and as a result the friction term starts reducing in significance due to its $1/\rho$ 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_{\rm tunl}$ much smaller, now allowing a parameter-space point that was excluded earlier.

4 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 $\tilde{y}$ and $M_{VL}$ 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_{\rm 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.3, Fig. 1, that when 3, 4 or 5 $SU(2)$ singlet VLQs all with 3 TeV mass are added, $\lambda(h)$ never goes negative, implying that the EW minimum is the global minimum and absolutely stable, unlike the SM situation. The reason for this behavior is explained in detail in Section 2.3. 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 less than 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 ($\tilde{y}$) with the Higgs can be written down. Nevertheless, when $\tilde{y}$ 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 $\tilde{y}=0.1$ and for $M_{VL}\lesssim 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 $\tilde{y}$, allows the intriguing possibility that the EW vacuum is rendered absolutely stable.

For example, the hidden-sector dark matter model in Ref. [47] 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.

5 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. [48, 49], which is,

[TABLE]

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 $\lambda(t)$, i.e. when $h_{0}$ is close to where $\beta_{\lambda}(h_{0})\approx 0$. Furthermore, when $S_{B}$ is so large that errors due to the transition not happening at a constant $\lambda(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. (31). We elaborate on this statement below with many examples.

In Fig. 9 in the left column we show $\beta_{\lambda}(\mu)$ vs. $\mu$ where $\mu\equiv h(\rho)$. In the right column we show the (absolute value of) integrand of Eq. (29) made dimensionless by multiplying the integrand by $1/m_{t}^{4}$ and denoted as $|\hat{I}(\hat{\rho})|$, vs. $\lambda(h(\rho))$, with $\rho$ being the parameter (not shown).

As shown in the topmost row in Fig. 9, for the SM, it is evident that most of the contribution to the integral comes from when $\lambda$ 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. (31) with the $t$ taken to be at the scale at which $\beta_{\lambda}=0$ where $\lambda=-0.009$, which gives $S_{B}^{\rm approx}=2848$. This is in excellent agreement with our numerical computation of $S_{B}$, and as discussed earlier, this is because $\beta_{\lambda}=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 $|\hat{I}(\hat{\rho})|$ for the SM, where most of the contribution to the $S_{B}$ integral is indeed coming for $\lambda=-0.009$, where the bounce spends most of its time ($\hat{\rho}$). Indeed, Eq. (31) was put-forth for the SM, where it can be safely applied.

As we see from the last three rows in Fig. 9, with VLF present, $\beta_{\lambda}$ 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. (31) but have to compute $S_{B}$ numerically. Indeed as the $|\hat{I}(\hat{\rho})|$ for these cases show, the integral gets its contributions for a range of $\lambda$. 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.

6 Conclusions

We study the stability of the electroweak vacuum in the presence of new vector-like fermions. We work with the renormalization group improved Higgs effective potential, identifying the Higgs field value $h\equiv\mu$. 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 $\tilde{y}$. We include the significant 2-loop contributions in the beta-functions. We numerically integrate the RGE to determine the scale at which $\lambda(\mu)$ becomes zero and goes negative.

The Higgs effective quartic-coupling $\lambda(\mu)$ 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_{\rm tunl}$ that the EW vacuum would have tunneled away by a true-vacuum bubble nucleating in our Hubble four-volume. Computing $P_{\rm tunl}$ requires computing the bounce configuration in Euclidean space-time, 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 $\lambda$ at a single scale where $\beta_{\lambda}(\mu)$ is approximately zero. This is because the bounce transition is mostly completed when $\lambda(\mu)$ 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 $\beta_{\lambda}(\mu)$ 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_{\rm tunl}$. We find that $P_{\rm tunl}$ is extremely sensitive to $\tilde{y}$ 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 $\tilde{y}$, $\lambda$ 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_{\rm tunl}$ than in the SM, or a smaller $P_{\rm 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.

Acknowledgements: We thank Romesh Kaul for valuable discussions.

Appendix A SM $\beta$-functions (fermion contributions)

Here, we review the calculation of the fermion contributions to the 1-loop $\beta$-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=(1/\sqrt{2})((\phi^{1}+i\phi^{2})\ \ (v+h+i\phi^{3}))^{T}$, where $v\approx 246$ GeV is the EW VEV, $h$ is the physical Higgs boson and $\phi^{i}$ are the Goldstone bosons. The Lagrangian density in terms of the bare fields $\{\tilde{h},t_{0}\}$ and bare coupling $y_{0}$ is

[TABLE]

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

[TABLE]

where the renormalized fields $h,t$ are defined by $\tilde{h}=\sqrt{Z_{h}}h$, ${t_{0}}_{L,R}=\sqrt{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)\equiv a_{h}y^{2}/2$, $(Z_{t_{L,R}}-1)\equiv a_{t_{L,R}}y^{2}/2$, $(Z_{y}-1)\equiv a_{y}y^{2}/2$, and we also define $(\hat{Z}_{y}-1)\equiv(Z_{y}\sqrt{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\equiv\hat{a}_{y}y^{2}/2$. Similarly the renormalized Lagrangian density for the Goldstone fields can be written down.

The Feynman vertices in momentum-space are as follows:

• propagator $(t_{L,R}t_{L,R})(p)$ is $i\hbox to0.0pt{\hbox to5.03125pt{\hfil/\hfil}\hss}p/(p^{2}+i\epsilon)$ with the counter-term $ia_{t_{L,R}}\hbox to0.0pt{\hbox to5.03125pt{\hfil/\hfil}\hss}py^{2}/2$,
• propagator $(hh)(p)$ is $i/(p^{2}+i\epsilon)$ with the counter-term $ia_{h}p^{2}y^{2}/2$,
• Yukawa coupling $ht\bar{t}$ is $-iy/\sqrt{2}$ with the counter-term $-i\hat{a}_{y}y^{3}/(2\sqrt{2})$,
• for the Goldstone bosons, the Yukawa coupling $\phi^{3}t\bar{t}$ is $-y\gamma^{5}/\sqrt{2}$ with the counter-term $\hat{a}_{y}y^{3}\gamma^{5}/(2\sqrt{2})$; and $\phi^{1,2}b\bar{t}$ is $iyP_{L}(1,i)/\sqrt{2}$ with the counter-term $-i\hat{a}_{y}y^{3}P_{L}(1,i)/(2\sqrt{2})$.

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 $\beta$-functions later.

The 1-loop correction to the Higgs 2-point function (wave-function renormalization of the $h$) due to the fermion is $i\Sigma_{h}(p_{h})=i(N_{c}y^{2}/(8\pi^{2}))p_{h}^{2}\left(1/\epsilon-\ln{p_{h}}-i\pi/2+\ln\sqrt{4\pi}-\gamma/2+{\cal O}(\epsilon)\right)$, after including a factor of $(-1)$ for the fermion loop, where $\gamma\approx 0.5772$ is the Euler-Mascheroni constant. To cancel this divergence, we fix the counter-term from the condition $i\Sigma_{h}(p_{0})+ia_{h}p_{0}^{2}y^{2}/2=0$ at the subtraction scale $p_{0}$. As mentioned above, since our goal is to derive the $t$ contributions to the $\beta$-functions, we show only the $y$ dependent terms in the counter-terms, and we omit other terms. This yields $a_{h}=-N_{c}/(4\pi^{2})\left(1/\epsilon-\ln{p_{0}}-i\pi/2+...\right)$, and we have

[TABLE]

The 1-loop correction to the fermion 2-point function (wave-function renormalization of the $t_{L,R}$) proportional to $y$ is $i\Sigma_{t_{L,R}}(p_{t})=iy^{2}/(32\pi^{2})\hbox to0.0pt{\hbox to5.03125pt{\hfil/\hfil}\hss}p_{t}\left(1/\epsilon-\ln{p_{t}/2}+\ln\sqrt{4\pi}-\gamma/2+{\cal O}(\epsilon)\right)$, and to cancel the divergence we fix the counter-term from the condition $i\Sigma_{t}(p_{0})+ia_{t}\hbox to0.0pt{\hbox to5.03125pt{\hfil/\hfil}\hss}p_{0}y^{2}/2=0$, which yields $a_{t_{L,R}}=-1/(16\pi^{2})\left(1/\epsilon-\ln{p_{0}/2}+\ln\sqrt{4\pi}-\gamma/2\right)$, and we have $i\Sigma_{t_{L,R}}(p_{t})+ia_{t_{L,R}}\hbox to0.0pt{\hbox to5.03125pt{\hfil/\hfil}\hss}p_{t}y^{2}/2=-iy^{2}/(32\pi^{2})\hbox to0.0pt{\hbox to5.03125pt{\hfil/\hfil}\hss}p_{t}\ln{p_{t}/p_{0}}$. The vertex 1-loop correction proportional to $y$ is

[TABLE]

where we take the Higgs momentum as $p_{h}$, and the fermion momenta as $-p_{h}/2$ and $p_{h}/2$. To cancel the divergence, we fix the vertex counter-term from the condition $iV(p_{0})-i\hat{a}_{y}=0$, which yields $iV(p_{h})-i\hat{a}_{y}y^{3}/(2\sqrt{2})=-iy^{3}/(16\sqrt{2}\pi^{2})\ln{p_{h}/p_{0}}$, and we have

[TABLE]

We discuss next the Goldstone boson contributions proportional to $y$. Starting with the self-energy corrections, we have the $\phi^{3}$ contribution to $\Sigma_{t_{L,R}}$ is equal to the $h$ contribution, the $\phi^{1}$ and $\phi^{2}$ contributions to $\Sigma_{t_{R}}$ are equal to the $h$ contribution, and, the $\phi^{1}$ and $\phi^{2}$ contributions to $\Sigma_{t_{L}}$ are proportional to $y_{b}$ which we drop and take to be zero. Turning next to the vertex corrections, we have the $\phi^{3}$ contribution to the $ht_{L}\bar{t}_{R}$ vertex ($V_{LR}$) is negative of the $h$ contribution to this vertex, the $\phi^{3}$ contribution to the $ht_{R}\bar{t}_{L}$ vertex ($V_{RL}$) is again negative of the $h$ contribution to this vertex, and, the $\phi^{1,2}$ contribution to $V_{LR,RL}$ is proportional to $y_{b}$ and hence we take it to be zero.

One way to extract the $\beta$-function is from the divergent part of the bare coupling. 333We briefly summarize here the method to obtain the $\beta$-function from the bare coupling, following ’t​ Hooft’s method as described in Ref. [41]. With $\kappa_{0}$ the bare-coupling, we write in $d=4-\epsilon$ dimensions, $\kappa_{0}\mu^{-\Delta(d)}=\kappa(\mu,d)-b(\kappa(\mu,d))/\epsilon$, $\Delta(d)\equiv\Delta-\rho\epsilon$, $\kappa_{0}$ being $\mu$ independent, and $\kappa(\mu,d)$ is the renormalized coupling. Then, we write $d\kappa(\mu,d)/d\ln{\mu}\equiv\beta_{\kappa}-\alpha_{\kappa}\epsilon$ with $\beta_{\kappa}$ being the $\beta$-function, and by matching powers of $\epsilon$, we obtain $\beta_{\kappa}=-\Delta\kappa-\rho b+\rho\kappa\partial b/\partial\kappa$. We generalize this to a system of many couplings $\kappa_{i}$ by writing ${\kappa_{0}}_{i}\mu^{-\Delta_{i}(d)}=\kappa_{i}(\mu,d)-\sum_{j}b_{ij}(\kappa(\mu,d))/\epsilon$ with $\Delta_{i}(d)\equiv\Delta_{i}-\rho_{i}\epsilon$. Then, the $\beta$-functions are $\beta_{\kappa_{i}}=-\Delta_{i}\kappa_{i}-\sum_{j}[\rho_{i}b_{ij}-(\partial b_{ij}/\partial\kappa_{j})\rho_{j}\kappa_{j}]$. For the couplings encountered here, we have: $\Delta_{y}=0$, $\rho_{y}=-1/2$; $\Delta_{\lambda}=0$, $\rho_{\lambda}=-1$; and $\Delta_{g_{a}}=0$, $\rho_{g_{a}}=-1/2$ (for $a=\{1,2,3\}$).

From the contributions computed above, we find the contributions proportional to $y$ to be $y_{0}\supset y+(y^{3}/(16\pi^{2}))((3+2N_{c})/2)(1/\epsilon)$, from which we obtain the fermionic contribution to $\beta_{y}$ as

[TABLE]

after including the $t,h,\phi^{1,2,3}$ contributions. This is in agreement with the results in Refs. [50, 4], for example. Interestingly, the $\phi^{1,2,3}$ contribute zero after including all their contributions.

The $yg_{3}^{2}$ and $yg_{2}^{2}$ contributions to $\beta_{y}$ can be written as [4]

[TABLE]

which are included in Eq. (8). To derive the $yg_{1}^{2}$ contribution, we start by extracting the relevant Feynman rules for the hypercharge gauge boson $B_{\mu}$ interactions. With all momenta going into the vertex, with $Y_{L,R}$ being the hypercharges of $\psi_{L,R}$ and $Y_{H}=1/2$ being the hypercharge of the Higgs doublet, we have the Feynman rules:

• $\phi^{3}(p_{3})h(p_{h})B_{\mu}$: $-g^{\prime}Y_{H}(p_{3}^{\mu}-p_{h}^{\mu})$;     $hhB_{\mu}B_{\nu}$: $2ig^{\prime 2}Y_{H}^{2}g_{\mu\nu}$;     $hB_{\mu}B_{\nu}$ : $2ig^{\prime 2}Y_{H}^{2}vg_{\mu\nu}$;
• $\psi_{L,R}\bar{\psi}_{L,R}B_{\mu}$: $ig^{\prime}Y_{L,R}\gamma^{\mu}$.

Computing the $B_{\mu}$ contribution at 1-loop order in the ’t​​ Hooft-Feynman $\xi\!\!=\!\!1$ gauge, we obtain the following divergent pieces: the $\psi_{R}\bar{\psi}_{L}h$ vertex correction due to $B_{\mu}$ exchange gives $iV^{(B_{\mu})}\supset-i8Y_{L}Y_{R}yg^{\prime 2}/(\sqrt{2}\,16\pi^{2}\epsilon)$; the Higgs 2-point function correction due to $\phi^{3}\--B_{\mu}$ exchange gives $i\Sigma_{h}^{(B_{\mu})}\supset-i4g^{\prime 2}Y_{H}^{2}p_{h}^{2}/(16\pi^{2}\epsilon)$; and the fermion 2-point function corrections due to $B_{\mu}$ exchange gives i\Sigma_{\psi_{L},\psi_{R}}^{(B_{\mu})}\supset i2g^{\prime 2}Y_{L,R}^{2}\mbox{\not p }/(16\pi^{2}\epsilon). We include in the counter-terms a piece to cancel these divergences at the subtraction scale $p_{0}$, i.e., $i(Z_{h}-1)p_{0}^{2}/2\supset-i\Sigma_{h}^{(B_{\mu})}(p_{0})$; i(Z_{\psi_{L},\psi_{R}}-1)\mbox{\not p }_{0}\supset-i\Sigma_{\psi_{L},\psi_{R}}^{(B_{\mu})}(p_{0}); and $i(y/\sqrt{2})(\hat{Z}_{y}-1)\supset iV^{(B_{\mu})}$. From these we determine $(Z_{y}-1)=-g^{\prime 2}(8Y_{L}Y_{R}+4Y_{H}^{2}-Y_{L}^{2}-Y_{R}^{2})/(16\pi^{2}\epsilon)$. Thus, since the bare coupling is $y_{0}=Z_{y}y$, we get the contribution

[TABLE]

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^{\prime}=\sqrt{3/5}g_{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 $\lambda_{0}\supset\lambda-N_{c}y^{4}/(8\pi^{2})\left(1/\epsilon+{\rm finite}\right)$. From this, we infer that this contribution leads to

[TABLE]

Writing the Higgs four-point vertex as $i\lambda_{\rm eff}$, the contribution to its evolution, i.e. $\beta_{\lambda}$, due to the fermion loop in the $h$-leg is just four times $\sqrt{\Sigma_{h}}$. Thus, for this contribution we have $i\lambda_{\rm eff}\supset(-i\lambda)(i/p^{2})(4\times i\Sigma_{\phi}/2)$, and from Eq. (34), we have $\lambda(\mu)=\lambda(M)+N_{c}\lambda y^{2}/(4\pi^{2})\ln(\mu/M)+...$, where $\mu$ is the renormalization scale, and $M$ is a subtraction scale. From this, and since $\beta_{\lambda}=d\lambda(\mu)/d\ln\mu$, we have

[TABLE]

We turn next to the $\beta$-functions of the gauge couplings $g_{a}=\{g_{3},g_{2},g_{1}\}$, focusing on the SM fermion contribution. We recall the definition $\beta_{a}=g_{a}^{3}b_{a}/(16\pi^{2})$. For $\beta_{g_{3}}$ we have the well-known result (see for example Ref. [41])

[TABLE]

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 $\mu$, we have $n_{3}=6$ for three generations of quarks, which implies $b_{3}=-7$. Similarly, for $\beta_{g_{2}}$ we have

[TABLE]

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 $\mu$, $n_{2}=(3N_{c}+3)$ for the three generations of quark and lepton doublets, which yields $b_{2}=-19/6$. Lastly, for $\beta_{g_{1}}$ we have

[TABLE]

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 $\phi$ with hypercharge $Y_{\phi}$. We recall that we use $SU(5)$ normalization for $g_{1}$, i.e. the SM hypercharge gauge coupling $g^{\prime}$ is related to $g_{1}$ by $g_{1}=\sqrt{5/3}g^{\prime}$. Thus in the SM, at large $\mu$, for three generations, $\sum_{f}Y_{f}^{2}=3\times(10/3)=10$, and for one Higgs doublet containing two complex fields with $Y_{H}=1/2$, $\sum_{\phi}Y_{\phi}^{2}=2(1/2)^{2}=1/2$, we get $b_{1}=41/10$. This agrees with, for example, Ref. [51].

We complete our derivation of the SM fermion contributions to the 1-loop $\beta$-functions. After adding the other contributions, the complete 1-loop $\beta$-functions are as given in Eqs. (7)-(9).

Appendix B VLF contributions to the RGE

Here we extend the $\beta$-functions derived in Sec. 2.1 and App. A to include VLF contributions, which we denote as $\beta^{VLF}_{\kappa}$.

We first derive the $\beta_{g_{a}}^{VLF}$. The $\beta_{g_{3}}^{VLF}$ is got easily from Eq. (40). Since the SM quark is vector-like with respect to $SU(3)$, we have an identical contribution for a VLQ, and we obtain the result shown in Eq. (15). For obtaining $\beta_{g_{2}}^{VLF}$, 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. (41) will give us $\beta_{g_{2}}^{VLF}$ as given in Eq. (16). Since for a VLF both $L$ and $R$ chiralities contribute, we take twice the first term in Eq. (42) to obtain $\beta_{g_{1}}^{VLF}$ as given in Eq. (17); the $2n_{2}$ is just the number of fermions in $n_{2}$ doublets having hypercharge $Y_{\chi}$.

Next, we derive the contributions present only when a full family is added, i.e. when the $\tilde{y}$ 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,\phi^{1,2,3}\}$ (all vertices have an overall $(-iy_{t}/\sqrt{2})$) written in terms of the Dirac spinors $t=(t_{L}\ t_{R})^{T}$ and $b=(b_{L}\ b_{R})^{T}$:

• $\{h,\phi^{3}\}t\bar{t}$ : $\{1,-i\gamma^{5}\}$ and $\{\phi^{1},\phi^{2}\}b\bar{t}$ : $\{-P_{L},-iP_{L}\}$.

Now, when a full VLF family is present, we have the $SU(2)$ doublet VLF $\chi=(\chi_{1}\ \chi_{2})^{T}$ and a singlet $\xi$. We can assemble the following Dirac spinors: $\psi_{1}=({\chi_{1}}_{L}\ \xi_{R})^{T}$, $\psi_{2}=(\xi_{L}\ {\chi_{1}}_{R})^{T}$, $\xi=(\xi_{L}\ \xi_{R})^{T}$ and $\chi_{2}=({\xi_{2}}_{L}\ {\xi_{2}}_{R})^{T}$. Using these, we can write the Feynman rules with the Higgs-doublet fields $\{h,\phi^{1,2,3}\}$ as follows (all vertices have an overall $(-i\tilde{y}/\sqrt{2})$):

• $\{h,\phi^{3}\}\psi_{1}\bar{\psi}_{1}$ : $\{1,-i\gamma^{5}\}$ and $\{\phi^{1},\phi^{2}\}\chi_{2}\bar{\psi}_{1}$ : $\{-P_{L},-iP_{L}\}$,
• $\{h,\phi^{3}\}\psi_{2}\bar{\psi}_{2}$ : $\{1,i\gamma^{5}\}$ and $\{\phi^{1},\phi^{2}\}\chi_{2}\bar{\psi}_{2}$ : $\{-P_{R},-iP_{R}\}$.

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 $\psi_{1}$ has identical couplings, while the second Dirac fermion $\psi_{2}$ has couplings that is similar but not identical, with a change $i\rightarrow-i$ in the $\phi^{3}$ couplings, and $P_{L}\rightarrow P_{R}$ in the $\phi^{1,2}$ couplings. We observe that all the diagrams that contribute to the $\beta$-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 $\beta$-function as in the SM. Thus, to obtain the VLF contributions to $\beta_{\lambda}^{VLF}$, we just multiply the SM contribution after combining Eqs. (38) and (39) by a factor of 2 and obtain Eq. (18). Next, the VLF contribution to $\beta_{y_{t}}$ is due to only the wavefunction renormalization contribution to $h$, i.e. $\Sigma_{h}$, and this contribution can be got from the second term in Eq. (35), but multiplied by 2 since there are two VLF Dirac sets as argued above, and changing the coupling to $y_{t}\tilde{y}^{2}$ instead of $y_{t}^{3}$, which then gives us the $\beta_{y_{t}}^{VLF}$ in Eq. (19). Next, consider the evolution of either the $h\psi_{1}\bar{\psi}_{1}$ coupling, or the $h\psi_{2}\bar{\psi}_{2}$ coupling, either of which is $\tilde{y}$. The VLF contribution to $\beta_{\tilde{y}}$ is due to these three contributions: (i) the vertex contribution proportional to $3\tilde{y}^{3}$ as in the first term in Eq. (35), (ii) the VLF contributions in $\Sigma_{h}$ which yields twice the second term in Eq. (35) proportional to $2\times 2N^{\prime}_{c}\tilde{y}^{3}$ since each of the $\psi_{1}$ and $\psi_{2}$ contribute as in the SM, and (iii) the top-quark contribution in $\Sigma_{h}$ which yields $2N_{c}y_{t}^{2}\tilde{y}$ as in the second term in Eq. (35). Adding these three contributions then gives the first part of $\beta_{\tilde{y}}$ in Eq. (20). We write the $\tilde{y}g_{a}^{2}$ contributions to $\beta_{\tilde{y}}$ following Eqs. (36) and (37), which gives the last part in Eq. (20). We thus complete the derivation of the VLF contributions to the $\beta$-functions given in Eqs. (15)-(20).