The Higgs boson inclusive decay channels $H \to b\bar{b}$ and $H \to gg$ up to four-loop level

The principle of maximum conformality (PMC) has been suggested to eliminate the renormalization scheme and renormalization scale uncertainties, which are unavoidable for the conventional scale setting and are usually important errors for theoretical estimations. In this paper, by applying PMC scale setting, we analyze two important inclusive Standard Model Higgs decay channels, $H\rightarrow b\bar{b}$ and $H\rightarrow gg$, up to four-loop and three-loop levels accordingly. After PMC scale setting, it is found that the conventional scale uncertainty for these two channels can be eliminated to a high degree. There is small residual initial scale dependence for the Higgs decay widths due to unknown higher-order $\{\beta_i\}$-terms. Up to four-loop level, we obtain $\Gamma(H\rightarrow b\bar{b}) = 2.389\pm0.073 \pm0.041$ MeV and up to three-loop level, we obtain $\Gamma(H\rightarrow gg) = 0.373\pm0.030$ MeV, where the first error is caused by varying $M_H=126\pm4$ GeV and the second error for $H\to b\bar{b}$ is caused by varying the $\overline{\rm MS}$-running mass $m_b(m_b)=4.18\pm0.03$ GeV. Taking $H\to b\bar{b}$ as an example, we present a comparison of three BLM-based scale setting approaches, e.g. the PMC-I approach based on the PMC-BLM correspondence, the $R_\delta$-scheme and the seBLM approach, all of which are designed to provide effective ways to identify non-conformal $\{\beta_i\}$-series at each perturbative order. At four-loop level, all those approaches lead to good pQCD convergence, they have almost the same pQCD series, and their predictions are almost independent on the initial renormalization scale. In this sense, those approaches are equivalent to each other.


I. INTRODUCTION
In 2012 a new boson has been discovered by CMS and ATLAS experiments at the Large Hadron Collider (LHC) [1,2], whose properties are remarkably similar to the Standard Model (SM) Higgs [3][4][5][6][7][8]. For example, its mass is 125.5 ± 0.2 +0.5 −0.6 GeV by ATLAS collaboration [8] or 125.7 ± 0.3 ± 0.3 GeV by CMS collaboration [4], where the first error stands for the statistic error and the second one stands for the systematic error.
If the SM Higgs has a mass around 126 GeV, then its decay width shall be dominated by H → bb [9,10]. Theoretically, many efforts have been made on studying the Higgs decays into a bottom pair [11][12][13][14][15][16][17][18][19][20][21][22][23][24]. As a reference, the pure QCD corrections at the two-loop or the three-loop level have been reported in Refs. [25][26][27], and the explicit expressions up to O(α 5 s ) have been given by Ref. [28]. In addition, the Higgs decay channel H → gg also plays a crucial role in studying the properties of Higgs boson. The coupling of Higgs to a pair of gluons, which is mediated at one loop by virtual massive quarks (essentially generated by the top quark alone) and becomes independent of the top-quark mass m t in the limit M H ≪ 2m t . An important feature of H → gg is its affinity to the gluon-gluon fusion mechanism for Higgs production. The effective coupling ggH eventually may provide a way to count the number of heavy quarks beyond SM [29]. Its next-to-leading order * email:wuxg@cqu.edu.cn (NLO) QCD corrections are quite large and amount to about 70% in comparison to the leading order contribution [30][31][32][33][34][35]. Later on, the QCD corrections for H → gg up to three-loop level have been calculated in the limit of an infinitely heavy top quark mass [36][37][38][39]. Those great improvements on loop calculations provide us chances for deriving more accurate estimation on Higgs properties.
The physical predictions of the theory, calculated up to all orders, are surely independent of the choice of renormalization scale and renormalization scheme due to the renormalization group invariance [40][41][42][43][44]. Whereas, as is well-known, there is renormalization scale and renormalization scheme ambiguities at any finite order. It is helpful to find an optimal scale setting so that there is no (or greatly suppressed) scale or scheme ambiguity at any fixed order and can achieve the most accurate estimation based on the known perturbative calculations.
As an estimation of the physical observable, one can first take an arbitrary initial renormalization scale µ r = µ init r and apply some scale setting method to improve the pQCD estimation. For the conventional scale setting, once the renormalization scale is setting to be an initial value, it will always be fixed during the whole analysis. That is, for the processes involving Higgs, in the literature, one usually take µ r ≡ µ init r = M H as the central value, which eliminates the large logs in a form as ln(µ r /M H ). Because there is no strong reasons for such a choice, as a compensation, one will vary the scale within a certain region, e.g. µ init r ∈ [M H /2, 2M H ], to ascertain the scale uncertainty. It is often argued that by setting and varying the scale in such a way, one can estimate contributions from higher order terms; i.e. changing in scale will affect how much of a result comes from Feynman diagrams without loops, and how much it comes from the leftover finite parts of loop diagrams. And because of its perturbative nature, it is a common belief that those scale uncertainties can be reduced after finishing a higher-and-higher order calculation. However, this ad hoc assignment of scale and its range usually constitutes an important systematic error in theoretical and experimental analysis. More explicitly, the conventional scale setting can not answer the questions: why it is M H and not M H /2 or others that provides the central estimation; when there are several typical energy scales for the process, then which one provides the central value ?
Several scale setting methods have been suggested, e.g. the renormalization group improved effective coupling method (FAC) [45], the principle of minimum sensitivity (PMS) [46], the Brodsky-Lepage-Mackenzie method (BLM) [47] and its underlying principle of maximum conformality (PMC) [48][49][50][51][52][53][54][55]. The FAC is to improve the perturbative series by requiring all higher-order terms vanish and the PMS is to force the fixed-order series to satisfy the renormalization group invariance at the renormalization point. The BLM improves the perturbative series by requiring the n f -terms of at each perturbative order vanish. The PMC provides the principle underlying the BLM, and it suggests a principle to set the optimal renormalization scales up to all orders, they are equivalent to each other through the PMC -BLM correspondence principle [49]. Those methods, being designed to eliminate the scale ambiguity, have quite different consequences and may or may not achieve their goals. A detailed introduction and comparison of these methods can be found in a recent review [55]. In the present paper, we adopt PMC for analyzing the Higgs decays.
The main idea of PMC lies in that the PMC scales at each perturbative order are formed by absorbing all nonconformal terms that governs the running behavior of the coupling constant into the coupling constant. At each perturbative order, new types of {β i }-terms will occur, so the PMC scale for each perturbative order is generally different. Even though, one can choose any value to be µ init r , the optimal PMC scales and the resulting finiteorder PMC prediction are both to high accuracy independent of such arbitrariness. After PMC scale setting, the divergent renormalon series does not appear and the convergence of the pQCD series can be greatly improved. Because of these advantages, the PMC method can be widely applied to high energy physics processes, some examples of which can be found in Refs. [49-52, 54, 56, 57].
Because the PMC provides the underlying principle for BLM, the previous features or properties derived by using BLM can also be understood by using PMC. Before applying PMC or BLM to high-energy processes, one needs to use the expression with full initial renormalization scale dependence. That is, those terms that have been eliminated by setting the renormalization scale to be equal to the factorization scale or by setting the initial renormalization to be the typical momentum transfer should be retrieved back. So, the previous scale dependence analysis or conclusions drawn under the BLM should be adopted with great care, since there is misuse of BLM in the literature. It is interesting to show whether the PMC can work well for the inclusive Higgs decays and whether the accuracy of the estimations can be improved. In present paper, we show the newly suggested PMC procedure, the so called R δ -scheme [54,58], with much more detail. A comparison with some other suggestions to extend the BLM scale setting up to any perturbative orders shall also be presented.
The remaining parts of the paper are organized as follows. In Sec.II, we present the calculation technology for applying PMC to Higgs decay processes H → bb and H → gg up to four-loop level. In Sec.III, we present the numerical results and discussions. The final section is reserved for a summary.

II. CALCULATION TECHNOLOGY
In this section, we present an improved analysis for the Higgs decay channels H → bb and H → gg by using the PMC R δ -scheme. For the purpose, • We shall first rearrange the four-loop expressions [28,38] that have been derived under the conventional scale setting in a more general form. That is, the n f -terms in those expressions that are coming from the light-quark loops and are responsible for controlling the running behavior of the coupling constant shall be transformed into {β i }series. Those {β i }-series via the renormalization group equation rightly control the running behavior of the strong coupling constant. Every process has its own {β i }-series and its own optimal (PMC) scales. Thus, after absorbing all those {β i }-series into the strong coupling constant via an order-byorder way, we can obtain the optimal running coupling constant for the specific process.
• As stated in the Introduction, before applying the PMC scale setting, we need to obtain the expressions with full (initial) renormalization scale dependence. For the purpose, at present, we need to transform the four-loop results of H → bb or the three-loop results of H → gg derived under the conventional assumption of µ init r ≡ M H to a more general form that explicitly contains the initial renormalization scale µ init r , which may or may not equal to M H . This can be achieved by using the strong coupling constant's scale transformation equation up to four-loop level, i.e.
where a s = α s /4π, Q * and Q are two arbitrary renormalization scales.
• We shall set the PMC scales in an order-by-order manner according to R δ -scheme. By doing the loop calculations, the b-quark mass is treated as massless, cf. the review on Higgs properties [59], there is only an overall m 2 b -factor in the decay width. Since its value is irrelevant to the PMC procedures and should be kept separate during the PMC scalesetting, either the choice of MS-running mass or the pole mass is reasonable. In the formulae, we fix its value to be m b (M H ) within the MS-scheme. A detailed discussion on this point is in preparation, which shows by applying PMC properly, ei-ther the choice of pole mass or MS-running mass can get consistent estimation. Recently, a discussion on H → γγ presents such an example [60].
In the following, we sequently present the results before and after PMC scale setting for the two channels H → bb and H → gg. By taking the initial renormalization scale µ init r = M H , the analytic decay width with explicit n f dependence for the channel H → bb can be formally written as The R δ -scheme [54,58] not only illuminates the {β i }pattern of the process but also exposes a special degeneracy of the coefficients at different perturbative orders. Such degeneracy is necessary, which, similar to the PMC -BLM correspondence principle [49], ensures the one-toone correspondence between the n f -series and the {β i }series at each perturbative order.
Applying R δ -scheme [54,58] to Eq.(2), one can derive the general form of H → bb for µ init r = M H with the help of Eq.(1), which can be written as, Here for later convenience, we have transformed the n fseries into the required {β i }-series. Explicit expressions for β 0 , β 1 and β 2 in n f -series can be found in Refs. [61][62][63].

B. Results for H → bb after PMC scale setting
Following the standard procedure of the R δ -scheme, by absorbing all non-conformal {β i }-series that control the running behavior of the coupling constant into the coupling constant, we can reduce Eq.(5) to the following conformal series, where in the second line, we present the values for the conformal coefficients over the α s expansion, which show the relative importance of the perturbative series. Here Q i with i = (1, · · · , 4) are PMC scales, which can be obtained through the following formulas where β = −a 2 s ∞ i=0 β i a i s being the conventional QCD renormalization group β-function. We note that the last scale leaves some ambiguity in PMC scale setting, since there is no {β i }-terms that can determine its optimal value, we set Q 4 = µ init r . Because the PMC scales (Q 1 , Q 2 and Q 3 ) themselves are in perturbative series, the residual scale dependence due to unknown higher-order {β i }-terms shall be highly suppressed.
C. The decay width of the H → gg By taking the initial renormalization scale µ init r = M H , the analytic expression for the decay width of H → gg can be written as For convenience, in the third line, we present the values for the coefficients over the α s -expansion by setting n f = 5, which explicitly show the relative importance of the perturbative series. The coefficients c i,j are [38] From Eq. (20), the general form for the Higgs decay process H → gg with µ init r = M H can be written as the following form, Following the same procedures of R δ -scheme, the conformal or non-conformal coefficients r i,j (µ init r ) can be written as It is noted that, as required, the conformal coefficients r i,0 are independent of µ init r and we can omit the argument (µ init r ) in r i,0 . By absorbing all non-conformal {β i }-series that control the running behavior of the coupling constant into the coupling constant, we can reduce Eq.(23) into the following conformal series, where in the second line, we present the values for the conformal coefficients over the α s expansion. Here r i,0 (µ init r ) are conformal coefficients. The PMC scales Q i with i = (1, 2, 3) can be obtained from the same way as that of Eqs. (17,18,19), Q 4 is also the undetermined scale due to the unknown higher order {β i } terms and we also fix it to be µ init r .

III. NUMERICAL RESULTS AND DISCUSSIONS
To do numerical calculation, we take G F = 1.16638 ×  [64], and by using the quark mass anomalous dimension expressions listed in Ref. [65,66], we obtain m b (M H ) = 2.78 ± 0.02 GeV.

A. Basic results for H → bb
The decay widths of H → bb before and after PMC scale setting are presented in Table I, where Γ i stands for the decay width at each perturbative order with i = LO, NLO and etc., and Γ tot stands for the total decay width. We set the renormalization scale µ r ≡ µ init r = M H for the conventional scale setting, and we take the initial renormalization scale µ init r = M H to initialize the PMC scale setting. Table I shows that, either before or after PMC scale setting, about 80% contribution comes from the LO order terms, which is exact and free from the strong interactions. The total decay width for H → bb remains almost unchanged, Γ(H → bb) ≃ 2.39 MeV. This shows the choice of µ r ≡ M H is a lucky guess for the conven-tional scale setting. Because of the elimination of the renormalon terms, one could expect a better pQCD convergence after PMC scale setting. This is clearly shown in Table I, e.g. the four-loop terms only give ∼ 0.1% contributions to the total decay width. For the H → bb process at O(α 4 s ), we need to introduce four PMC scales, i.e. the LO PMC scale Q 1 , the NLO PMC scale Q 2 , the N 2 LO PMC scale Q 3 and the N 3 LO PMC scale Q 4 . As stated in the last section, since there is no β-terms that can determine its optimal value, we set Q 4 ≡ µ init r 1 . Using the formulas (17,18,19), we 1 This corresponds to the second type of residual scale dependence after PMC scale setting [57], which, as is the present case, can also be highly suppressed when the pQCD convergence is under well control. show how the PMC scales depend on the initial renormalization scale, which are presented in FIG.(1). FIG.(1) shows that the PMC scales Q 1,2,3 are highly independent on the choice of initial renormalization scale. This indicates that the PMC scale setting do provide a principle for setting the optimal (solitary) renormalization scale of high energy processes. For example, setting µ init r = M H , we obtain Q 1 = 62.3 GeV, Q 2 = 40.5 GeV, Q 3 = 53.1 GeV. (35) These PMC scales are smaller than M H to a certain degree due to the exponential suppressions from the absorbtion of higher order {β i }-terms. These PMC scales are different, which shows that they are controlled by different {β i }-series at each perturbative order.
As a further comparison, we show the total decay width Γ(H → bb) versus the initial renormalization scale µ init r before and after PMC scale setting in FIGs. (2,3 process or a value which minimizes the contributions of the loop diagrams. For the present process, µ r ≡ µ init r = M H . The total decay width Γ(H → bb) shows a relatively strong dependence on the value of µ init r at the NLO level, e.g. We would like to stress that even if a proper choice of µ init r may lead to a value close to the experimental data by using conventional scale setting, this may not be the correct answer for a fixed-order estimation. Especially, if a process does not converge enough, one has to finish a more and more complex loop calculations so as to achieve the precision goal.
• FIG.(3) shows that after PMC scale setting, the total decay width of H → bb are almost flat versus the choice of renormalization scale even at the NLO level. This is due to the fact that the PMC scales Q 1 , Q 2 and Q 3 themselves are highly independent on the choice of µ init r , as shown by FIG.(1). The residual scale dependence due to unknown higherorder {β i }-series has been highly and exponentially suppressed.
It is noted that there is slight difference for the decay widths at different perturbative orders, e.g.
Such difference shows that even though by absorbing the non-conformal terms into the coupling constant, one can greatly improve the pQCD convergence and simultaneously eliminate the scale dependence at lower perturbative orders, one may still need to know higher-order conformal contributions if one wants to achieve even higher precision. For examples, the N 2 LO improves NLO estimation by about 2% and the N 3 LO improves N 2 LO estimation by about 1%. More over, the unknown higher-order non-conformal contributions can be roughly estimated by varying the final undetermined PMC scale as Q 4 via the conventional way, e.g. [Q 4 /2, 2Q 4 ].

B. Basic results for H → gg
We can estimate the properties of H → gg in a similar way as that of H → bb.
We put the decay width for H → gg before and after PMC scale setting in Table II, where Γ i stands for the decay width at each perturbative order with i = LO, NLO and etc., and Γ tot stands for the total decay width.  (38) and after PMC scale setting, we have This shows that for the decay of H → gg, only after a three-loop correction, one can obtain a desirable convergent estimation. Note the pQCD convergence after PMC is weaker than the case for the conventional scale setting for H → gg, since the N 2 LO part becomes (−20%). This could mean that we need to know an more accurate {β i } series so as to determine the N 2 LO PMC scale. For the H → gg process at O(α 5 s ), we need to introduce four PMC scales, i.e. the LO PMC scale Q 1 , the NLO PMC scale Q 2 , the N 2 LO PMC scale Q 3 and the N 3 LO PMC scale Q 4 . As stated in the last section, since there is no β-terms that can determine its optimal value, we set Q 4 ≡ µ init r . FIG.(4) shows that the PMC scales Q 1,2,3 are highly independent of the choice of initial renormalization scale. This indicates that the PMC scale setting do provide a principle for setting the optimal (solitary) renormalization scale of high energy processes. By setting µ init   These PMC scales are also smaller than M H to a certain degree due to exponential suppressions from the absorbtion of higher order {β i }-terms. We show the total decay width Γ(H → gg) versus the initial renormalization scale µ init r before and after PMC scale setting in FIGs.(5,6). In these two figures, the dash-dot, dashed, dotted and solid lines are for LO, NLO, N 2 LO and N 3 LO estimations, respectively.
These results show that for the decay channel of H → gg, we can obtain similar conclusions as those of H → bb. More explicitly, • FIG.(5) indicates that, under the conventional scale setting, the total decay width Γ(H → gg) shows a strong dependence on µ init r at the LO level, e.g. N 2 LO and N 3 LO estimations, respectively. This shows that as one includes higher-and-higher orders, the guessed scale will lead to a better estimation. For the H → gg decay, only after a three-loop correction, one can obtain a desirable small about several percent scale error.
• FIG.(6) shows that, after PMC scale setting, the total decay widths of H → gg up to the mentioned perturbative orders are almost flat versus the choice of initial renormalization scale µ init r . The residual scale dependence due to unknown higher-order {β i }-series has been highly suppressed. Similar to the case of H → bb, this is due to the fact that the PMC scales Q 1 , Q 2 and Q 3 themselves are highly independent on the choice of µ init r , which are shown clearly by FIG.(4).

C. A comparison of the approaches underlying BLM scale setting
The BLM scale setting is designed to improve the pQCD predictions by absorbing the n f -terms via a proper way into the coupling constant [47]. Since its in-vention by Brodsky-Lepage-Mackenzie in 1983, the BLM has been widely accepted in the literature for dealing with high energy processes, such as the e + e − → hadrons, the deep inelastic scattering, the heavy meson or baryon productions or decays, the exclusive processes such as the pion-photon transition form factors, the QCD lattice perturbative theory, and etc.. Encouraged by its great successes, several approaches have been tried to extend BLM to any perturbative orders or put it in a more solid background, such as the PMC-I approach (first approach to achieve the goal of PMC via the PMC-BLM correspondence principle) [48,49], the R δ -scheme (second approach to achieve the goal of PMC) [54,58] and the seBLM approach [67,68].
It is noted that the role of the running coupling in any gauge theory is to absorb the physics of the β function, which governs its running behavior via the renormalization group equations. Any approach that properly identifies the {β i }-series for a physical observable will lead to equivalently the same result. Practically, one usually calculate the n f -terms by considering the vacuum polarization contributions. However, different ways of identifying n f -series to {β i }-series may lead to: I) different effective {β i }-series at each known perturbative order; II) different residual {β i }-dependence because of unknown perturbative orders; III) different conformal terms leaving at each perturbative order; IV) and different pQCD convergence. If one can do the perturbative corrections to a higher enough perturbative order, different effective schemes may result in consistent physical predictions. The equivalence of the PMC-I approach and R δ -scheme have already been shown in Refs. [54,55,58]. In the following, we take the Higgs decay channel H → bb as an explicit example to show that the seBLM approach is also consistent with the PMC approaches. One subtle point of such a comparison (an also any applications of those approaches) lies in that, we should first transform the estimations with full renormalization scale dependence with the help of the transformation formulae (1).
In PMC-I approach, by introducing an PMC-BLM correspondence principle in which the {β i }-series for a physical observable has the same parton of the running coupling itself [49], the number of the effective independent {β i }-terms exactly corresponds to the number of n f -terms at each perturbative order. In R δ -scheme, by introducing the "degeneracy" properties of the {β i }series observed by a generalization of the conventional dimensional regularization scheme to any dimensional-like ones [54], one can also obtain a one-to-one correspondence between the {β i }-series and the n f -series. The calculation technologies for those two self-consistent approaches can be found in the corresponding references, the interesting readers may turn to those references or very recent review [55] for detail. In Sec.II, we have presented our analysis under the R δ -scheme.
While the seBLM scheme provides quite a different way from those two PMC approaches, in which a general {β i }series at each perturbative order have been introduced, and in order to get an one-to-one correspondence between the n f -series with the {β i }-series, some extra approximations (or equivalently some extra degrees of freedom) have to be introduced [67]. More explicitly, the seBLM scheme transforms the standard power series a n s (µ init r ) to the series of the products n i=1 a s (Q i ). After applying the seBLM scheme to Eq.(2), the decay width of the process H → bb can be expressed as follows: where r ′ i stands for the conformal coefficients of seBLM leaving at each perturbative order, A i = β 0 a s (Q i ) stand for the redefined coupling constant. Because we have no higher-order {β i }-terms to determine the scale for A 4 , we set Q 4 = Q 3 as suggested by seBLM. Three effective seBLM scales are ln where the explicit expressions of the coefficients ∆ i,j can be found in Ref. [67] 2 . Then, we obtain Q 1 = 62.5 GeV, Q 2 = 29.1 GeV, Q 3 = 127.0 GeV (45) for µ init r = M H = 126 GeV. It is noted that the expressions for Q 1 in Eq. (17) and Eq.(42) are equal to each other at the LO and NLO level, we obtain almost the same value for Q 1 under the R δ -scheme and the seBLM approach.
We present a comparison of those three approaches in Table III, which shows that the pQCD convergence of the perturbative series clearly. Here as a comparison, the estimations for conventional scale setting and the results for the BKM scheme [69] and the fractional analytic perturbation theory (FAPT) scheme [70] are also presented.
At the present considered estimation up to four-loop levels all of those schemes including the conventional scale setting lead to good pQCD convergence. Especially, the PMC-I scheme, the R δ -scheme and the seBLM scheme have almost the same pQCD series. This is reasonable, since those three scale-setting schemes are designed to deal with the {β i }-series of the process. In addition, one will observe that after eliminating the non-conformal {β i }-series, one may also derive the (initial) renormalization scale independence for a fixed order prediction for those approaches. The initial scale dependence for the PMC-I, R δ -scheme and seBLM approaches on the NLO decay widths Γ NLO are presented in Table  IV, where there typical initial scales µ init r = M H /2, M H and 2M H are adopted. It shows clearly that the value of Γ NLO are almost unchanged with µ init r . The higher-order terms have similar properties.

IV. SUMMARY
The conventional scale setting procedure assigns an arbitrary range and an arbitrary systematic error to fixedorder pQCD predictions. And its error analysis can only get a rough estimation of the β-dependent nonconformal terms, not the entire perturbative series. As a possible solution, the PMC provides a systematic way to set the optimized renormalization scales for high energy processes. In principle, the PMC needs an initial value to initialize renormalization scale and renormalization procedures. It is found that the estimates of PMC are to high accuracy independent of the initial renormalization scale; even the PMC scales themselves are in effect independent of the initial renormalization scale and are 'physical' at any fixed order. This is because the PMC scale itself is a perturbative series and those unknown higherorder {β i }-terms will be absorbed into the higher-order terms of the PMC scale, which is strongly exponentially suppressed. Since the renormalization scale and scheme ambiguities are removed, the PMC can improve the precision of tests of the Standard Model and enhances the sensitivity to new phenomena. It is noted that the PMC applies the known non-conformal {β i }-terms in a strict and systematic way to determine the behavior the coupling constant at each perturbative order. It provides an accurate estimation for the known perturbative series, and one may still need higher order calculations so as to known even higher-order conformal contributions, especially, when the perturbative series is not converge enough.
The PMC can be applied to a wide-variety of perturbatively-calculable collider and other processes. In addition to previous examples done in the literature, following its standard R δ -scheme procedures, we have done a through analysis of these two processes up to four-loop and three-loop levels. A comparison of the estimations under three approaches, i.e. the PMC-I approach, the R δ -scheme and the seBLM approach, have also been presented. We observe, • Under conventional scale setting, it is often argued that by finishing a higher enough perturbative calculation, one can get desirable convergent and scale-invariant estimations. For the present considered channels, when considering up to three-loop level or four-loop level, as shown by FIGs.(2,5), the decay width becomes almost invariant within the region of µ init However, even if a proper choice of µ init r may lead to a value close to the experimental data by using conventional scale setting, this is a guess work and may not be the correct answer for a fixed-order estimation. Especially, if a process does not converge enough, one has to finish a more and more complex loop calculations so as to achieve the precision goal.
• As shown by FIGs.(3,6), after PMC scale setting, the total decay widths of H → bb and H → gg show a fast trend of stabilization versus the change of initial renormalization scale, which are almost flat even at the NLO level. The residual scale dependence due to unknown higher-order {β i }-series has been greatly suppressed. This indicates that the PMC scale setting do provide a principle for setting the optimal renormalization scale of high energy processes.
• In comparison to the previous estimations, e.g.
Refs. [71,72], a more accurate predications for those two Higgs decay processes can be obtained. After PMC scale setting, we obtain the total decays widths for those two channels:  .(7).
• The β function governs the running behavior of the coupling constant via the renormalization group equations, thus any approach that can properly identify the {β i }-series for a physical observable will surely lead to equivalently the same predictions. Practically, one usually calculate the n fterms by considering the vacuum polarization contributions. After a certain scale setting, different ways of identifying n f -series to {β i }-series may lead to different effective {β i }-series at each known perturbative order, different residual {β i }-dependence because of unknown perturbative orders, or different pQCD convergence. A comparison of H → bb for the PMC-I approach, the R δ -scheme and the seBLM approach has been presented in Table III. At the four-loop level all those approaches lead to good pQCD convergence, they have almost the same pQCD expansion series, and all of them are almost independent on the wide choice of the initial renormalization scale. This shows these three approaches are equivalent to each other. The residual differences of these approaches are caused by the unknown {β i }-terms that could be suppressed to a required accuracy by finishing a more higher-order calculation.
• As one subtle point, one may meet the problem of quite small (or near the fixed point) effective scales for a specific scale setting method. For example, we have noted that for the case of R(e + e − → hadron) at the measured scale Q, we can obtained a convergent and precise conformal series up to fourloop level by applying the R δ -scheme, whose LO, NLO and NNLO PMC scales are [54]: Q 1 = 1.3Q, Q 2 = 1.2Q and Q 3 = 5.3Q, respectively. In contrast, by using the seBLM scheme, we shall obtain ln(Q 2 /Q 2 2 ) ∼ 167 [67], which leads to quite small Q 2 out of pQCD domain. If, as suggested by PMC, we only deal with the n f series that rightly controls the running behavior of the coupling constant into the coupling constant via the standard way of se-BLM, then we shall obtain more moderate seBLM scales, Q 1 = 1.3Q, Q 2 = 1.1Q and Q 3 = 228.9Q. And, similar to the present Higgs decays, we can obtain consistent results for R(e + e − ) under both the seBLM and the PMC scale settings. Moreover, it is noted that such situation could be softened to a certain degree by using the commensurate scale relation [73], or one may solve it by using proper running behavior of the coupling constant in lower scale region [74].
• As another subtle point, even if one can eliminate the scale dependence at lower perturbative order as NLO, it may necessary to know the higher-order conformal contributions if we want to achieve even higher precision. Taking the case of H → bb as an example, its N 2 LO terms improves the NLO estimation by about 2% and the N 3 LO terms improves the N 2 LO estimation by about 1%. The unknown higher-order non-conformal contributions can be roughly estimated by varying the final undetermined PMC scale as Q 4 via the conventional way, e.g. [Q 4 /2, 2Q 4 ]. If after PMC scale setting the final terms at a certain fixed order give negligible contribution, then we shall obtain quite accurate estimations at such fixed order.
As shown by Table IV, it is noted that by setting µ init r = M H /2, the H → bb NLO decay width Γ NLO under the conventional scale setting is close to the PMC estimations. In this sense, a choice of µ r ≡ M H /2 is better than the choice of µ r ≡ M H for the conventional scale setting. In fact, under such choice, one can also obtain a more convergent pQCD series for the conventional scale setting.
• As a minor point, taking H → bb as an example, we point out a wrong way of estimating the conventional scale error. The correct way is to set an initial scale µ init r and get the full µ init rdependent expression (5), e.g. those terms proportional to ln(µ init r /M H ) are kept, and then by varying µ init r ∈ [M H /2, 2M H ] to discuss its scale error. In this way, we have shown that the conventional scale error up to four-loop level is almost eliminated, as shown in FIG.(2). The wrong way is to adopt the expression (2) as the starting point, and directly varies the scale of coupling constant from M H to M H /2 or 2M H to discuss the scale error. In this wrong treatment, the log terms involving µ init r = M H disappear, which however may have sizable contributions. In fact, such a naive treatment shows the conventional scale error is still about ±2% [71,72] for varying the scale within the region of [M H /2, 2M H ] even at the four-loop level.