Leading QCD-induced four-loop contributions to the $\beta$-function of the Higgs self-coupling in the SM and vacuum stability

We present analytical results for the leading top-Yukawa and QCD contribution to the $\beta$-function for the Higgs self-coupling $\lambda$ of the Standard Model at four-loop level, namely the part $\propto y_t^4 g_s^6$ independently confirming a result given in [1]. We also give the contribution $\propto y_t^2 g_s^6$ of the anomalous dimension of the Higgs field as well as the terms $\propto y_t g_s^8$ to the top-Yukawa $\beta$-function which can also be derived from the anomalous dimension of the top quark mass. We compare the results with the RG functions of the correlators of two and four scalar currents in pure QCD and find a new relation between the anomalous dimension $\gamma_0$ of the QCD vacuum energy and the anomalous dimension $\gamma_m^{SS}$ appearing in the RG equation of the correlator of two scalar currents. Together with the recently computed top-Yukawa and QCD contributions to $\beta_{g_s}$ [2,3] the $\beta$-functions presented here constitute the leading four-loop contributions to the evolution of the Higgs self-coupling. A numerical estimate of these terms at the scale of the top-quark mass is presented as well as an analysis of the impact on the evolution of $\lambda$ up to the Planck scale and the vacuum stability problem.


Introduction
The evolution of the Higgs self-coupling and of the Higgs field are important ingredients for the Renormalization Group (RG) improved Higgs potential and the study of vacuum stability. A precise determination of the Higgs self-coupling in the Standard Model extended up to the Planck scale is important because this parameter is close to zero at the Planck scale and the question whether the SM vacuum state is stable or not can only be answered definitively by reducing the uncertainties. 1 The largest source of uncertainty is the experimentally measured top mass M t . At a future linear e + e − collider it could, however, be measured with a precision which matches that of the theory input to the vacuum stability analysis (see figure 5 in [13]).
On the theory side there are three sources of uncertainty. The first is the difference between the effective Higgs potential V eff (Φ cl ) [22] and the approximation of the RG-improved potential where Φ cl = 0|Φ|0 is the classical field strength of the scalar SU(2) doublet Φ, γ Φ the anomalous dimension of Φ and µ 0 the scale where we start the evolution of fields and couplings, e. g. µ 0 = M t . This uncertainty is negligible at large values of Φ cl , e. g. close to

JHEP06(2016)175
the Planck scale [23][24][25][26][27]. In this approximation the SM vacuum is stable up to the scale Λ ∼ M Planck if λ(µ) > 0 for µ ≤ Λ. Another source of theoretical uncertainties is the matching of experimental parameters, e. g. M t , α s (M Z ), M H , to the parameters of the SM Lagrangian, y t (µ 0 ), g s (µ 0 ), λ(µ 0 ), . . . at some initial scale µ 0 renormalized in the MS-scheme. State of the art is the full numerical two-loop matching [16,28]. In order to improve precision here three-loop calculations might be attempted and different mass definitions than the pole mass could be used for the top mass.
In this paper we improve on the third source of uncertainty, namely we increase the precision in the β-functions, calculated in the MS-scheme. The β-function for a coupling X is defined as (1.2) and the anomalous dimension of a field f as where Z f is the field strength renormalization constant, where X, X 1 , X 2 , . . . are the couplings of the theory which we want to include in the analysis. The RG functions of the SM were computed at three-loop accuracy during the last years [10,[29][30][31][32][33][34][35]. The four-loop β-function for the strong coupling g s was first computed in pure QCD [36,37] and recently extended to the gaugeless limit of the SM, namely to include the dependence on the top-Yukawa coupling y t and the Higgs self-coupling λ [2,3]. The leading four-loop contribution to the Higgs self-coupling β-function was first presented in [1] and is independently confirmed in this paper.
The paper is structured as follows: in the next section we briefly describe the technical details of the calculation. Then the leading four-loop terms for β λ , β yt and γ Φ 2 are given and the relevance of the four-loop terms numerically determined at the scale of the top quark mass. Finally, we investigate the impact of the new contributions on the evolution of λ in order to estimate the uncertainty reduction due to four-loop β-functions.

The model: QCD plus minimal top-Yukawa contributions
For this calculation we start with the SM Lagrangian in the broken phase where . (2.1) The UV renormalization constants in the MS-scheme do not depend on masses and are the same as in the unbroken phase. Hence we can use all renormalization constants determined up to three-loop level in previous calculations [10] and set all masses to zero. There is no γ 5 in the ttH-vertex ∝ y t as opposed to the ttΦ 2 -vertex of the unbroken phase.  only appears in the Yukawa vertices with χ and Φ ± . As at low scales (where we start the evolution of couplings and fields) the strong coupling is the largest we take as the leading contribution to the vertex and self-energy corrections those where H only appears as an external field. This means that no Higgs or Goldstone propagators appear. The electroweak gauge-couplings as well as λ and all Yukawa couplings except y t are neglected. For the top-Yukawa vertex and the top self-energy these are the pure QCD corrections (see figure 1). For the Higgs self-energy and the quartic Higgs-vertex these are gluon insertions into the one-loop diagram ∝ y 2 t (see figure 3(a)) and ∝ y 4 t (see figure 2(a)) as well as diagrams with two fermion loops (see figure 3(b) and figure 2(b)). Thus we get the fourloop contributions which are numerically most significant to the evolution of λ avoiding γ 5 and its treatment in D = 4 dimensions completely.

JHEP06(2016)175
We compute the field strength renormalization constants Z (HH) 2 from the Higgs and Z (tt) 2 from the top self-energies. The quartic Higgs-vertex is renormalized with Z (4H) 1 and the top-Yukawa-vertex with Z (ttH) 1 . 2 From these we compute and All divergent integrals are regularized in D = 4 − 2ε space time dimensions and the renormalization constants are defined as Z = 1 + δZ in the MS-scheme.

Calculation with massive tadpole integrals
For the computation of the four-loop terms we use the setup described in detail in [3].
The generation of all necessary Feynman diagrams was done with QGRAF [38]. The C++ programs Q2E and EXP [39,40] are then used to identify the topology of the diagram. The Taylor expansion in external momenta, the fermion traces and the insertion of counterterms in lower loop diagrams was performed with FORM [41,42]. All colour factors were computed with the FORM package COLOR [43].
In the momentum space part of the diagrams we introduce the same auxiliary mass parameter M 2 in every propagator denominator. The self-energy diagrams are then expanded to second order in the external momentum q after applying a projector ∝ / q to the top self-energy diagrams and taking the trace over the external fermion line. Then we divide by q 2 before q is set to zero. In all vertex correction diagrams we can set q → 0 from the beginning. This is allowed as MS renormalization constants do not depend on external momenta. After this we are left with tadpole integrals. Subdivergences ∝ M 2 are canceled by counterterms computed from and inserted in lower loop diagrams. This is the same method for computing UV renormailzation constants as in our previous calculations [3,10,32]. It was first introduced in [44] and then further developed in [45]. A detailed explanation of the calculation of Z-factors with an auxiliary mass can be found in [14]. Up to three-loop order the tadpole integrals were computed with the FORM-based package MATAD [46]. The four-loop tadpoles are reduced to Master integrals using FIRE [47,48]. The needed four-loop Master integrals can be found in [37]. 2 In the notation of [10] these renormalization constants are Z (HH)

JHEP06(2016)175 3 Analytical results
In this section we give our results which can be found in machine readable format on http://www-ttp.particle.uni-karlsruhe.de/Progdata/ttp16/ttp16-008/ For a gerneric SU(N c ) gauge group the colour factors are expressed through the quadratic Casimir operators C F and C A of the fundamental and the adjoint representation of the corresponding Lie algebra. The dimension of the fundamental representation is called d R . The adjoint representation has dimension n g and the trace T F is defined by T F δ ab = Tr T a T b with the group generators T a of the fundamental representation. Higher order invariants are constructed from the symmetric tensors from the generators of the fundamental representation and analogously d abcd 3) The number of active fermion flavours is denoted by n f . The leading four-loop contributions to the β-functions for the Higgs self-coupling and the top-Yukawa coupling are found to be

Comparison with available QCD results
All diagrams discussed in the previous section are special in one aspect: they comprise just the minimal number of non-QCD vertexes, that is one for β yt , two for γ Φ 2 and four for β λ .

JHEP06(2016)175
Even more, these non-QCD vertexes are of one and the same type, namely the insertion of the scalar top-quark currentt t. This means that the corresponding anomalous dimensions should be related to some RG functions in pure QCD describing the QCD evolution of the scalar current(s). The corresponding "effective" Lagrangian implies, obviously, the following identities valid in all orders in α s (dots below stand for terms which have a different dependence on the SM coupling constants than y t α n s (first line) and y 2 t α n s (second line) correspondingly): Here γ m (α s ) is the quark mass anomalous dimension and the function γ SS q (α s ) appears in the evolution equation for the scalar correlator (B marking bare quantities) which is renormalized as The anomalous dimensions are found to be Needless to say that a comparison with available results for γ m (α s ) [49,50] and γ SS q (α s ) [51] confirms the relations (4.2) and (4.3) at four-loop order. 4 Finally, let us consider in some detail the last (and somewhat more complicated) case, viz. the β-function for the Higgs self-coupling β λ . The corresponding renormalization constant coincides up to a factor y 4 t 4·4! to that which renormalizes the (1PI) Green's function of the T-product of four scalar currents:

JHEP06(2016)175
With m t = 0 we could nullify all external momenta p i in eq. (4.9) and, thus, consider all j s -operators on the r.h.s. of (4.9) as insertions of the scalar quark current at zero momentum transfer. As is well-known such insertions can be generated by multiple differentiations of QCD Green functions w.r.t. a quark mass. 5 This means that the corresponding anomalous dimensions should be related to some pure QCD RG functions. The well-known way to construct the corresponding relations is to use the (renormalized) Quantum Action Principle [58,59]. Let us briefly outline the main points.
The Quantum Action Principle relates properties of (regularized) Lagrangian and the full Green's functions. Consider the generating functional of (connected) Green's functions The Action Principle states (in particular) that where λ is a any parameter in the Lagrangian L. The action principle works for DR Green functions [60] (modulo axial anomalies). An example: the (renormalized) QCD Lagrangian with n l massless quarks and a massive (top) one is customarily written as With properly chosen renormalization constants Z i (4.13) should produce finite Green's functions. If one differentiates W (J) w.r.t. the (renormalized) top quark mass m t the functional should also be finite. However, let us consdider eq. (4.14) at J = 0. It corresponds, obviously, to the VEV of the operator Z m Z 2t t(x), which is not finite already at order α 0 s (a couple of typical diagrams contributing to (4.14) at J = 0 are shown on figure 4). This means that our QCD Lagrangian (4.13) is not full : the term responsible for the renormalization of the vacuum energy is missing. The full QCD Lagrangian reads [61]

JHEP06(2016)175
here E 0 (µ) is the (renormalized) vacuum energy and Z 0 is the corresponding renormalization constant. Now, a four-fold differentiation of the generating functional (4.11) (with L = L full QCD ) w.r.t. m t immediately leads us to the conclusion that the combination should be finite. As a result, we arrive at the following identity valid in all orders in α s : where the dots stand for terms which have a different dependence on the SM coupling constants than y 4 t α n s and 6 is the anomalous dimension of the vaccuum energy The vacuum anomalous dimension plays an important role in the description of the renormalization mixing of all three scalar gauge-invariant operators with (mass) dimension four: It was proven in [61,62] that the matrix of anomalous dimensions in (4.21) reads:

JHEP06(2016)175
In addition, the two-point scalar correlator (4.3) at q = 0 is obviously related (via the action principle) to the four-point correlator (4.9), as the latter can be obtained from the former by a double differentiation w.r.t. m t . This, in turn, leads to the following remarkable relation (again valid in all orders in α s ): Thus, one could compute γ 0 in a few different ways.
1. Direct renormalization of the vacuum energy diagrams. This was done for two and three loops in the papers [61] and [63] respectively. At four loops it was first found in this way for a space-time dimension D = 3 [64,65] (in the process of computing the free energy in the effective high temperature QCD) and (implicitly, via eq. (4.17)) in the present paper.
2. By renormalizing the 4-loop scalar correlator [66] in the limit of small quark mass. 7 The result was later used in [67] to compute the quartic mass corrections to R had at O(α 3 s ).
3. By computing the lowest moment of the scalar correlator at 4 loops [68].
4. By computing the lowest low-energy moment of the axial-vector correlator (related via a Ward identity to the VEV of the scalar current) [69].
Note, finally, that the last two evaluations dealt with massive tadpole diagrams and that all calculations of γ 0 at the four-loop level described above are in mutual agreement as well as the two results for β (4) λ : the one displayed in (3.4) and the one obtained via relation (4.17).

Numerical analysis
In this section we want to numerically evaluate the results for the β-functions presented in section 3. The couplings g s , g 2 , g 1 , y t and λ in the MS-scheme at some fixed scale µ 0 can be computed by matching the experimentally measured parameters α s (M Z ), G F , M W , M Z , M t and M H to them [16,28], where the uncertainties of G F , M W and M Z have a negligible influence on the MS-couplings as compared to the other three. The Higgs mass is taken to be M H = 125.09 ± 0.24 GeV [70] and the strong coupling is extracted from α s (M Z ) = 0.1181 ± 0.0013 [71].
The top mass is a more difficult subject. At the moment a theoretically well-defined mass is not available to high prescision. The extraction of an MS-mass from cross section measurements leads to an uncertainty of 4 − 5 GeV [71], the extraction of the pole mass from cross section measurements to [71] M t = 174.6 pole ± 1.9 GeV. (5.1)

JHEP06(2016)175
The most precise top mass measurements from LHC and TEVATRON give M MC t ≈ 173.34 ± 0.76 GeV [72] where the Monte Carlo mass parameter M MC t would correspond to the pole mass in a purely perturbative setup. Since M MC t is affected by real emission and the implementation of the parton shower it eludes a theoretical definition based on the Lagrangian but has to be calibrated to fit the experimental data. The pole mass of a heavy quark also suffers from a conceptual problem, namely the renormalon ambiguity [73,74]. This is due to the fact that a quark is not observed as a free particle and non-perturbative effects spoil the convergence of the relations between the pole mass and e. g. an MS-mass. One way to connect the Monte Carlo and the pole mass is to identify the first with a theoretically well-defined short-distance mass at a scale R ∼ 1 − 3 GeV, i. e. the parton shower cut-off of the Monte Carlo simulation [75]. This so-called MSR mass does not suffer from non-perturbative effects such as renormalons due to the IR cutoff. The MSR mass can in turn be connected to the pole mass (by adding the contributions from the IR region) [76]  For the following analysis we take this value which is very close to the Monte Carlo mass, noting however that the discrepancy with the current PDG top pole mass value (5.1) indicates that the uncertainties on this parameter might be even larger and a direct and precise extraction of theoretically well-defined quantities like the MSR mass or the MS top-Yukawa coupling itself from the experimental data of a linear collider will be necessary in order to match the size of the other uncertainties entering the vacuum stability analysis. For M t = 173.39 +1.12 −0.98 GeV [72,77] we get the couplings in the MS-scheme at the scale of the top mass using two-loop matching relations [16]  where the experimental uncertainty (exp) stems from M t , M H and α s (M Z ) and the theoretical one (theo) from the matching of on-shell to MS parameters [16]. Evaluating β y t at the scale M t we find

JHEP06(2016)175
which shows the expected suppression of higher order contributions. The four-loop terms are significantly smaller than the lower loop terms. For β λ however the picture is different. Already in previous works [10,32,78] we found a slow convergence of the perturbation series up to three-loop order at the electroweak scale and this is also found true at four-loop level: As discussed in [32,78] there is a remarkable cancellation between the terms containing only g s and y t at two-loop 8 and three-loop level at the scale of the top mass making the overall contribution much smaller than the size of the individual terms. These cancellations seem to accidental, so we cannot predict there existence at four-loop level. However, if the cancellation also exists at this order the terms ∝ y 6 t g 4 s , ∝ y 8 t g 2 s , ∝ y 10 t , etc could make the result significantly smaller and hence increase the convergence of the perturbative series. At higher scales, where g s and y t become smaller the convergence is of course better also if we only include the term ∝ y 4 t g 6 s . At µ = 10 9 GeV we find where the leading four-loop contribution is a factor 4 − 5 smaller than the three-loop result but still not completely negligible. We will now check what this means for the evolution of the Higgs self-coupling λ up to the Planck scale.

Evolution of λ and vacuum stability
We evolve λ from the scale M t using the initial conditions (5.3) and the full SM β-functions (including g s , y t , g 2 , g 1 , λ) up to three-loop order and at four-loop level β    The difference between the evolution of λ with three-loop β-functions (blue curve) and including the leading four-loop terms (red curve) should give some indication on the uncertainty stemming from the truncation of the perturbative series for the β-functions. In order to see this difference we have to zoom in. We choose to do this at the scale where λ becomes negative, which is shown in figure 6. The error bands are calculated using the partial four-loop results and hence centered around the red curve.
In order to further illustrate the dependence of the vacuum stability problem on the top mass and the importance of the issues related to the top pole mass as an input parameter we show the evolution of λ also for the PDG value extracted from cross section measurements of the top pole mass (5.1) and the corresponding uncertainties in figure 7.
The conclusion for vacuum stability remains the same as in our previous works [10,11,[13][14][15]. It looks as if λ becomes negative at around log 10 µ GeV ∼ 9.6 (or log 10 µ GeV ∼ 8.7 in figure 7) rendering the SM not stable if extended up to scales above, but a definitive answer is pending on a more precise extraction of y t (µ 0 ) from experimental data. It is worth noting, however, that due to the reduction in the top mass uncertainty since the combined LHC and TEVATRON analysis [72] a stable SM up to the Planck scale is strongly disfavoured. In this work we have presented analytical results for the leading four-loop contributions to the β-function for the Higgs self-coupling λ and the top-Yukawa coupling as well as to the anomalous dimension of the Higgs field. These results have been connected to pure QCD RG functions and a relation between the anmalous dimension γ 0 of the vacuum energy and γ SS m was found. We have performed an analysis of the evolution of the Higgs self-coupling updating the analyses presented in previous works [10,11,[13][14][15] and establishing a nice hirarchy between the different sources of uncertainty.
With the computation of the leading four-loop terms to β λ , β yt and β gs the β-function uncertainty to the question of vacuum stability becomes significantly smaller than the matching uncertainty (before the two were comparable) which is in turn significantly smaller than the experimental top mass uncertainty. We expect that a full calculation of four-loop β-functions in the SM will confirm this conclusion rendering the remaining β-function uncertainty almost negligible in comparison to the other sources of uncertainty for a vacuum stability analysis.