Gauge dependence of the perturbative QCD predictions under the momentum space subtraction scheme

The momentum space subtraction (MOM) scheme is one of the most frequently used renormalization schemes in perturbative QCD (pQCD) theory. In the paper, we make a detailed discussion on the gauge dependence of the pQCD prediction under the MOM scheme. Conventionally, there is renormalization scale ambiguity for the fixed-order pQCD predictions, which assigns an arbitrary range and an arbitrary error for the fixed-order pQCD prediction. The principle of maximum conformality (PMC) adopts the renormalization group equation to determine the magnitude of the coupling constant and hence determines an effective momentum flow of the process, which is independent to the choice of renormalization scale. There is thus no renormalization scale ambiguity in PMC predictions. To concentrate our attention on the MOM gauge dependence, we first apply the PMC to deal with the pQCD series. We adopt the Higgs boson decay width, $\Gamma(H\to gg)$, up to five-loop QCD contributions as an example to show how the gauge dependence behaves before and after applying the PMC. It is found that the Higgs decay width $\Gamma (H\to gg)$ depends very weakly on the choices of the MOM schemes, being consistent with the renormalization group invariance. It is found that the gauge dependence of $\Gamma(H\to gg)$ under the $\rm{MOMgg}$ scheme is less than $\pm1\%$, which is the smallest gauge dependence among all the mentioned MOM schemes.


I. INTRODUCTION
Quantum chromodynamics (QCD) is believed to be the field theory of hadronic strong interactions. Due to its asymptotic freedom property [1,2], the QCD strong coupling constant becomes numerically small at short distances, allowing perturbative calculations for the highenergy processes. The QCD theory in a covariant gauge with massless quarks has three fundamental propagators which are for the gluon, the ghost and the quark fields, respeectively, and four fundamental vertices, namely the triple-gluon, the four-gluon, the ghost-gluon and the quark-gluon vertices. In the literature, various renormalization schemes have been adopted to regularize and remove the ultraviolet divergences emerged at higher perturbative orders. Among them, the momentum space subtraction (MOM) scheme [3][4][5][6][7][8] has also been frequently used in addition to the conventional minimum substraction scheme [9], which carries considerable information of various quark and gluon interaction vertices at spe-cific momentums and leads to a better convergence for some cases. Initially, the MOM scheme is defined via renormalizing the three-point vertices of the QCD Lagrangian at the completely symmetric point [3,4], i.e. the squared momentum of each external momentum of the vertex is equal. Lately, the asymmetric point with one of the external momentum vanishes for the threepoint vertex has been suggested [7,10,11], which has the property of avoiding the infrared divergence in massless QCD theory. More explicitly, the minimal MOM (mMOM) scheme [7] which subtracts at the asymmetric point where one external momentum vanishes has been suggested as an alternation of the original symmetric MOM scheme. It is an extension of the MOM scheme on the ghost-gluon vertex and allows the strong coupling to be fixed solely through a determination of the gluon and ghost propagators. Further more, there are other four kinds of asymmetric MOM schemes, e.g. the one with vanishing momentum for the incoming ghost in the ghost-gluon vertex, the one with vanishing momentum for the incoming quark in the quark-gluon vertex, and the two schemes in dealing with the case of vanishing momentum for the incoming gluon in the triplegluon vertex, respectively. Following the same notions as those of Ref. [11], we label the first two MOM schemes as MOMh and MOMq schemes, and the other two schemes as MOMg and MOMgg schemes [12,13], respectively. Even though the MOM schemes have been successfully applied in various high-energy processes, in different to the minimum substraction scheme, it has been found that the MOM scheme breaks down the gauge invariance. It is interesting to show whether the gauge dependence exists for all (typical) kinds of MOM schemes, or to find a MOM scheme with minimum gauge dependence.
The strong coupling is the most important component of the pQCD theory, we need to know its exact magnitude at any scale so as to derive an accurate pQCD prediction. The scale running behavior of the strong coupling is controlled by the renormalization group equation (RGE), or the β-function. The RGE for the MOM scheme can be related to the one under the modified minimal subtraction scheme (e.g. the MS scheme [14]) via proper relations. At the present, the explicit expressions of the {β i }-functions under the MS scheme have been known up to five-loop level in Refs. [15][16][17][18][19][20][21][22][23][24][25]. Thus the five-loop {β i }-functions for the MOM schemes (mMOM, MOMh, MOMq, MOMg and MOMgg) can be determined with the help of the known five-loop relations [7,11,26,27] to the MS scheme. Another way of deriving the running behavior of the MOM strong coupling up to five-loop level can be found in Ref. [28]. A key component for solving the β-function is the QCD asymptotic scale Λ. The asymptotic scale in MS-scheme can be fixed by using the PDG world average of the strong coupling constant at the scale of Z 0 boson mass, α MS s (M Z ) = 0.1181 ± 0.0011, which leads to Λ n f =5 MS = 0.210 ± 0.014 GeV [29]. The asymptotic scale for a MOM scheme can be derived by using the Celmaster-Gonsalves relation [3][4][5][6][7], e.g.
The MOM scheme could be a useful alternative to the MS scheme for studying the behavior and truncation uncertainty of the perturbation series. Many MOM applications have been done in the literature, e.g. two typical MOM applications for the Higgs-boson decays to gluons and the R-ratio for the electron-positron annihilation can be found in Refs. [6,[30][31][32][33]. Moreover, the processes involving three-gluon or four-gluon vertex provides an important platform for studying the renormalization scale setting problem. For the three-gluon vertex, it has already been pointed out that the typical momentum flow which appears in the three-gluon vertex should be a function of the virtuality of three external gluons [34]. As an example, because of the improved convergence, a more accurate and reliable pQCD prediction for Pomeron intercept can be achieved under the MOM scheme other than the MS scheme [35][36][37][38]. The MOM scheme can also be helpful to avoid the small scale problem emerged in MS scheme [39,40] The Higgs boson is a crucially important component of the Standard Model (SM), its various decay channels are important components for Higgs phenomenology. Among those decay channels, the decay width of H → gg have been calculated up to five-loop level under the MS scheme [41][42][43][44][45][46][47][48][49][50]. Using the relations among the strong coupling constants under various renormalization schemes, one can obtain the corresponding fiveloop MOM expression for the Higgs boson decay width Γ(H → gg) from the known MS expression. A way to transform the pQCD predictions from one renormalization scheme to another renormalization scheme has been explained in detail in Ref. [51]. In the paper, we shall adopt the decay width Γ(H → gg) up to five-loop QCD contributions as an explicit example to show how the gauge dependence of the MOM prediction behaves with increasing known perturbative orders.
Following the standard renormalization group invariance (RGI), a physical observable (corresponding to an infinite order pQCD prediction) should be independent to the choices of renormalization scale and renormalization scheme. For a fixed-order pQCD prediction, conventionally, people uses guessed renormalization scale together with an arbitrary range to estimate its uncertainty, which leads to the mismatch of strong coupling constant with its coefficient at each order and then results as conventional renormalization scheme-and-scale ambiguities. Many scale setting approaches have been suggested to solve the renormalization scale ambiguity. Among them, the principle of maximum conformality (PMC) [52][53][54][55][56] has been suggested to eliminate the conventional renormalization scheme-and-scale ambiguities simultaneously. In different to other scale-setting approaches such as the RG-improved effective coupling method [57,58] and the Principle of Minimum Sensitivity [59][60][61][62][63] and the sequential BLM [64,65] or its alternated version Modified se-BLM [66], the purpose of PMC is not to find an optimal renormalization scale but to fix the running behavior of the strong coupling constant with the help of the RGE, whose argument is called as the PMC scale. The PMC scale is physical in the sense that its value reflects the "correct" typical momentum flow of the process, which is independent to the choice of renormalization scale. After applying the PMC, the convergence of the pQCD series can be greatly improved due to the elimination of divergent renormalon terms. The PMC has a solid theoretical foundation, it satisfies the standard RGI and all the self-consistency conditions of the RGE [67]. Detailed discussions and many applications of the PMC can be found in the reviews [68][69][70][71]. In the paper, we shall first adopt the PMC to eliminate the renormalization scale ambiguity and then discuss the gauge dependence of the MOM predictions on the decay width Γ(H → gg).
The remaining parts of the paper are organized as follows. In Sec.II, we give the basic components and the formulas for transforming the strong coupling constant from various MOM schemes to MS scheme, which are important to transform the known MS pQCD series to MOM one. In Sec.III, we give a brief review on the PMC single-scale approach, which shall be adopted to do our present PMC analysis. In Sec.IV, we discuss the gauge dependence of the decay width Γ(H → gg) under the above mentioned five asymmetric MOM schemes. Sec.V is reserved for a summary. Some detailed formulas are given in the Appendix.

II. THE MOMENTUM SPACE SUBTRACTION SCHEMES
The scale dependence of the strong coupling is controlled by the following β-function, β(a(µ)) = µ 2 ∂a(µ) where µ is the renormalization scale, a(µ) ≡ α s (µ)/(4π). The {β i }-functions are scheme dependent, and their expressions up to five-loop level under the MS-scheme are available in Refs. [15][16][17][18][19][20][21][22][23][24][25]. For short, when there is no confusion, we set a = a(µ) in the following discussions. For an arbitrary renormalization scheme R, the respective renormalization of the gluon, quark and ghost fields are of the form where Z R 3 , Z R 2 andZ R 3 are the renormalization constants of the gluon field A, the quark field ψ, and the ghost field c, respectively. The superscripts 'B' and 'R' denote the bare and the renormalized fields, respectively. The superscript 'b' is the color index for the adjoint representation of the gauge group.
By using the usually adopted dimensional regularization [72] (we work in D = 4 − 2ǫ dimension), the renormalized strong coupling a and the gauge parameter ξ can be written as follows: where we have used the fact that the gauge parameter is also renormalized by the gluon field renormalization constant. The bare strong coupling is scale invariant, and the D-dimensional β-function for the renormalized strong coupling can be derived by doing the derivative over both sides of Eq. (7): Then, we obtain The renormalization of the gluon, ghost and quark selfenergies can be performed as follows and the renormalization of the triple-gluon, the ghostgluon and the quark-gluon vertexes can be performed as follows: where the vertex renormalization constants are related to the field and coupling renormalization constants via the Ward-Slavnov-Taylor identities (i.e. the generalized Ward-Takahashi identities [76][77][78][79]) by Under the minimal subtraction scheme (MS) [9] in which the ultraviolet divergence (1/ǫ-terms) in pQCD series are directly subtracted, the renormalized parameters Z MS k can be written as where the coefficients b MS m,n are free of µ-dependence [73]. The renormalized constant Z MS a is gauge independent, which takes the following form Here the {β i }-functions are for the MS scheme, which are the same for all the other dimensional-like renormalization schemes. This is due to the fact that the strong coupling among the dimensional-like schemes can be simply related via a scale shift [55], e.g. the MS scheme differs from the MS scheme by an additional absorbtion of ln 4π − γ E , which corresponds to redefining the MS scale µ MS as µ 2 MS = µ 2 MS exp (ln 4π − γ E ). Gross and Wilczek found that the LO {β i }-functions under the dimensionallike renormalization schemes are gauge independent [74], and lately, Caswell and Wilczek gave a proof of such gauge independence up to all orders [75] 2 .
Using Eq.(8), one obtains the following relations for the strong coupling and gauge parameter between the MOM and MS schemes: It has been found that the MOM scheme is gauge dependent. In MOM scheme [7,11,26], the gluon, ghost and quark self-energies are absorbed into the field renormalization constants at the subtraction point q 2 = −µ 2 : Using Eq.(8), we obtain the following relationship of the gauge parameters under the MS scheme and MOM scheme: In the following subsections, we make a simple introduction of five asymmetric MOM schemes, giving the relations of the strong couplings under those schemes with the one under the conventional MS scheme, and their gauge-dependent basic components, which are done by renormalizing the three-point vertices, such as the ghostgluon, the gluon-quark and the triple-gluon ones, at the asymmetric point with one of the external momentum of the vertex vanishes, respectively. The gluon, the quark and the ghost propagators, as shown by Fig.1, take the following form where a and b are color indices, i and j denote quark flavors. The gauge parameter ξ = 0 is the Landau gauge, ξ = 1 is the Feynman gauge, and etc. The self-energies Π(q 2 ), Σ V (q 2 ) andΠ(q 2 ) can be extracted from the corresponding one-particle irreducible diagrams by applying proper projection operators [11] (the same holds for the vertex functions discussed below).

B. The ghost-gluon vertex
The tree-level ghost-gluon vertex is −ig s f abc q µ , where q µ is the outgoing ghost momentum. There are two possibilities to set one of the external momenta to zero for the ghost-gluon vertex. One is to set the gluon momentum to zero, whose diagram is shown by Fig.2(a), and the renormalized vertex can be written as and the other is to set one of the incoming ghost momentum to zero, whose diagram is shown by Fig.2(b), and the renormalized vertex can be written as Γ abc µ (q; −q, 0) = −ig s f abc q µΓh (q 2 ).
The ghost-gluon vertex: (a)Γ abc µ (0;-q, q) for the case of the incoming gluon has zero momentum and (b)Γ abc µ (q;q, 0) for the case of one incoming ghost has zero momentum.
Here,Γ g (q 2 ) orΓ h (q 2 ), is the Lorentz invariant function with vanishing gluon or ghost momentum, respectively. At the tree-level, we havẽ The MOMh scheme is defined by renormalizing the ghost-gluon vertex ( Fig.2 h (q 2 = −µ 2 ) = 1,(33) Using Eqs. (15,16,17,18,22,23,24,25), we can connect the strong coupling in the MOMh scheme to the one in the MS scheme through the following equation: In addition, motivated by the non-renormalization of the ghost-gluon vertex in the Landau gauge [79], the vertex renormalization constant for this vertex is chosen as the same as that in MS, i.e.
which is equal to 1 in the Landau gauge. We can then derive the following relation for the coupling constants in those two schemes, We put the derivation in Appendix A. It shows that the MOMh scheme is equivalent to mMOM scheme for the Landau gauge (ξ mMOM = ξ MOMh = 0).

C. The quark-gluon vertex
There are two non-trivial cases with vanishing incoming external momentum for the quark-gluon vertex, the case of a vanishing incoming gluon momentum as shown by Fig.3(a)  as shown by Fig.3(b). It is clear that nullifying the incoming quark momentum is equal to the result of nullifying the outgoing quark momentum, therefore, the vertex Fig.3(a) can be written as and the vertex Fig.3(b) can be written as The subscript 'g' in Eq. (37) and 'q' in Eq. (38) indicate the functions with vanishing gluon momentum and incoming quark momentum, respectively. T a ij is the SU(3) color group generator for the quark. At the tree-level, we have The MOMq scheme is defined by renormalizing the quark-gluon vertex with vanishing incoming quark momentum, e.g., Therefore, the relation of the coupling constants in the MOMq scheme and the MS scheme is

D. The triple-gluon vertex
The triple-gluon vertex is symmetric under the exchange of any two of the gluons. As shown in Fig.4, one can set the momentum of the right-hand gluon to zero without loss of generality. Under this condition, the triple-gluon vertex generally takes the following form where f abc is the structure constant of the SU(3) color group. T 1 (q 2 ) corresponds to tree-level vertex, i.e.
is always absent at tree-level but arises from radiative corrections. T 3 (q 2 ) vanishes due to the Ward-Slavnov-Taylor identity for the triple gluon vertex [11,26]. The MOMg scheme is defined by renormalizing the above triple-gluon vertex with vanishing incoming gluon momentum, e.g., Therefore, the relation of the coupling constants in the MOMg scheme and the MS scheme is Another MOM scheme, which is also based on the triple-gluon vertex is MOMgg scheme, is defined by the following renormalization condition [12,13]: This gives the following coupling relations between MOMgg scheme and the MS scheme: At the present, the gluon self-energies Π MS A (−µ 2 ), the ghost self-energiesΠ MS c (−µ 2 ), the quark self-energies Σ MS V (−µ 2 ), the quark-gluon vertex with vanishing incoming quark momentum Λ MS q (−µ 2 ), the ghost-gluon vertex with vanishing incoming ghost moment, and the two functionsΓ MS h (−µ 2 ) and T MS 1 (−µ 2 ) defined in Eq.(45) and the function T MS 2 (−µ 2 ) defined in Eq. (47) have been calculated up to four-loop QCD corrections under the MS scheme, c.f. Refs. [11,26]. Using the formulas given in those two references and using the relations (34,36,42,45,47) and the equation (26), we can obtain the expressions for the strong couplings and the gauge parameters under various MOM schemes. For convenience, we put those relations in the Appendix B.
Those relations are helpful to transform the conventional MS series to the one under a specific MOM scheme. They are also important to get the MOM scheme βfunction and hence determine the correct α s -running behavior in MOM scheme. The MOM β-function is explicitly gauge dependent, which can be written as The anomalous dimension of gauge parameter γ MS is the gluon field anomalous dimension. Therefore, the MOM-scheme β-function takes the form: It has been stated that the property of the gauge invariance of the renormalization schemes is a sufficient but not a necessary property for the factorization of the QCD β-function [80]. Thus a reliable α s -behavior can be determined and hence a reliable pQCD prediction for various MOM schemes by applying a proper scale-setting approach to deal with the {β i }-terms of the process.

III. GENERAL PMC ANALYSIS OVER THE PERTURBATIVE SERIES
Conventionally, a pQCD approximant, δ(Q), of a physical observable takes the form where Q represents the scale at which the observable is measured, the index p indicates the α s -order of the leading-order (LO) prediction. Here the perturbative coefficients C i are usually in n f power series, where n f is the number of light flavors involved in the process. Using the degeneracy relations among different orders [55,56,81], the pQCD series can be written as the following β i series: where r i,0 (i = 1, 2, 3 . . .) are conformal coefficients which are generally free from renormalization scale dependence, and r i,j (1 ≤ j < i) are non-conformal coefficients, r m,n = r m,n | µ=Q . As a subtle point, any n f -terms that are irrelevant to determine the α s -running behavior should be kept as a conformal coefficient and cannot be transfomred into {β i }-terms [68].
For the standard PMC multi-scale approach described in Refs. [53,55], one needs to absorb the same type of {β i }-terms at various orders into the strong coupling constant via an order-by-order manner. Different types of {β i }-terms as determined from the RGE lead to different running behaviors of the strong coupling constant, and hence, determine the distinct PMC scales at each order. Because the precision of the PMC scale for highorder terms decreases at higher orders due to the less known {β i }-terms in its higher-order terms. Due to the unknown perturbative terms, the PMC prediction has residual scale dependence [37], which is however quite different from the arbitrary conventional renormalization scale dependence. The PMC scale, reflecting the correct momentum flow of the process, is independent to the choice of renormalization scale, and its resultant residual scale dependence is generally small due to both the exponential suppression and the α s suppression [71]. As an alteration, the PMC single-scale approach has been suggested to suppress the residual scale dependence [82]. It effectively replaces the individual PMC scales derived under the multi-scale approach by a single scale in the sense of a mean value theorem. The PMC single scale can be regarded as the overall effective momentum flow of the process; it shows stability and convergence with increasing order in pQCD via the pQCD approximates. The prediction of the PMC single-scale approach is schemeindependent up to any fixed order [83], thus its value satisfies the standard RGI. The examples collected in Ref. [83] show that the residual scale dependence emerged in PMC multi-scale approach can indeed be greatly suppressed. In the present paper, we adopt the PMC singlescale approach to do our discussions.
Using the standard procedures for the PMC singlescale approach, we can eliminate all the non-conformal {β i }-terms and rewrite Eq.(50) as the following conformal series: where the PMC scale Q is fixed by requiring all the nonconformal {β i }-terms vanish. The perturbative series of ln Q 2 /Q 2 over a(Q) up to next-to-next-to-next-to-leading log (N 3 LL) accuracy takes the following form: For convenience, we put the perturbative coefficients λ i (i = 0, 1, 2, 3) in the Appendix C. One may observe that both the resultant PMC conformal series (52) and the scale Q are free of renormalization scale (µ), and thus, the conventional renormalization scale dependence has been eliminated. There is residual dependence for δ(Q) due to the unknown terms (e.g. the unknown N 4 LLterms and higher) in the perturbative series (53).

IV. GAUGE DEPENDENCE OF THE TOTAL DECAY WIDTH Γ(H → gg)
In the present section, we adopt the total decay width Γ(H → gg) under various MOM schemes as an explicit example to show how the gauge dependence behaves with increasing known perturbative orders before and after applying the PMC.
Up to α 6 s -order level, the decay width of H → gg takes the following form where a(µ) = α s (µ)/(4π), µ is the renormalization scale, M H is the Higgs boson mass, and G F = 1.16638 × 10 −5 GeV −2 is the Fermi coupling constant. The coefficients C i∈[0,4] (M H ) under the MS scheme can be read from Refs. [41][42][43][44][45][46][47][48][49][50]. Those coefficients are usually given in n f -power series, and before applying the PMC, the perturbative series (54)   The effective scalesQ (Left column) and their corresponding coupling constants α s (Q) (Right column) versus the gauge parameter are shown in Fig.5. At the two-loop level, the determined scaleQ is gauge independent; However at the three-loop level or even higher, Fig.5(a1-a5) show that Q is gauge dependent under the five mentioned MOM schemes. This is due to the fact that all the r i,j terms of MOM schemes are gauge dependent except for those of i − j = 1 terms. For a specific gauge, we observe that the difference between the two nearby values of Q becomes smaller when more loop terms have been included, indicating the precision ofQ is improved by knowing more loop terms, which agrees with the perturbative nature ofQ. At the present,Q can be fixed up to N 3 LL-accuracy, which is of high accuracy, e.g. the N 3 LL-term shall only shift the N 2 LL-accurateQ by ∼ +1 GeV for mMOM, MOMh  parameter |ξ MOM | for all those MOM schemes. Fig.5(b1-b5) show that the effect coupling (α s (Q)) is also gauge dependent for all the five MOM schemes, the only exception is the MOMgg scheme whose gauge dependence is small and is even zero at the two-loop level. Numerically, the scaleQ and the effect coupling α s (Q) are almost the same for the three MOM schemes, e.g. the MOMh, the MOMq and the MOMg schemes; and if the magnitude of the gauge parameter |ξ MOM | for those three schemes are less than 1, the differences between the two nearby values of α s (Q) at different orders are almost unchanged for a fixed gauge parameter, indicating the effective couplings for those three schemes quickly achieves its accurate value at lower orders. It is interesting to find that under the MOMgg scheme, the effective coupling α s (Q) also can achieve its accurate value at lower orders, whose value is almost gauge independent for |ξ MOMgg | ≤ 1, indicating the gauge dependence ofQ is well compensated by the gauge dependence of Λ MOMgg QCD . More explicitly, the asymptotic scale for various MOM schemes can be derived from Λ MS by using the Celmaster-Gonsalves relation [3][4][5][6][7]. We present the ratios of Λ MOM /Λ MS for n f = 5 in Fig. 6, where the ratios for the following mentioned three symmetric MOM schemes are also presented. It is found that the asymptotic scales of MOMg, MOMq, and MOMh schemes are the same, together with the close values of Q, one can explain the close behaviors of α s (Q) and then close Γ(H → gg) under those schemes. Almost all of the ratios show explicit gauge dependence; the only exception is the ratio of MOMgg scheme, whose value is free of ξ MOMgg and is fixed to be e 50/69 ≈ 2.06 for n f = 5.

B. The gauge dependence of Γ(H → gg) for the five asymmetric MOM schemes
We present the total decay width Γ(H → gg) up to two-loop, three-loop, four-loop and five-loop level before and after applying the PMC in Figs. (7,8,9,10,11). Agreeing with conventional wisdom, the conventional renormalization scale dependence, estimated by varying µ ∈ [M H /4, 4M H ], becomes smaller when more loop terms have been included, e.g. the shaded bands becomes narrower when more loop terms have been included. At the same time, one may observe that the PMC predictions under various MOM schemes are scale independent at any fixed order, but becomes more accurate when more loop terms have been included. Thus the conventional renormalization scale uncertainty is eliminated by applying the PMC, which is consistent with previous PMC examples done in the literature. As mentioned in the Introduction, the scale independence of the PMC prediction is reasonable, since the determined scaleQ re- flects the overall typical momentum flow of the process which should be independent to the induced parameters.
Figs. (7,8,9,10,11) show that the gauge dependence cannot be eliminated by including more and more higherorder terms for both the conventional and the PMC scale-setting approaches. More explicitly, the total decay width of H → gg up to five-loop level under the mMOM, MOMh, MOMq, MOMg, and MOMgg schemes are presented in Tables I and II inate it by including higher order terms. As discussed in Introduction, the MOM schemes have some advantages in dealing with the perturbative series; And as will be shown below, the scheme independence of variance MOM predictions can be greatly suppressed at higher orders. It is thus interesting to find in which region of ξ MOM , the MOM prediction is more reliable. Being consistent with the suggestion of Ref. [10], Figs. (7,8,9,10,11) show that the gauge dependence is relatively weaker within the region of ξ MOM ∈ [−1, 1]. As an exception, Fig.11 indicates that at enough higher orders, the MOMgg prediction could be unchanged for a wide range of |ξ MOMgg |.Thus among all the MOM schemes, we prefer the MOMgg scheme. More explicitly, the coefficient functions C i≥2 in Eq.(50) and r i≥2,0 in Eq. (51) are series in powers of ξ MOM , thus a small magnitude such as |ξ MOM | ≤ 1 could lead to a better convergence and a more steady prediction over the change of ξ MOM ; this is why the Landau gauge with ξ MOM = 0 is usually adopted in the literature [4]. We also agree that it is better to choose a smaller value of |ξ MOM | for various MOM schemes. And in the following discussion, we shall adopt ξ MOM ∈ [−1, 1] to do our discussion.

Conv.
is free of gauge dependence at the two level, which shall be changed by about 1% for n = 3, 4, 5. Thus the gauge dependence of the MOMgg scheme is the smallest. And by taking µ ∈ [M H /4, 4M H ], the total decay width for mMOM, MOMh, MOMq and MOMg schemes behave closely, which shall be changed by about 45%, 10%, 5% and 3% for n = 2, 3, 4, 5, respectively; and the total decay width Γ(H → gg)| MOMgg
Moreover, by using the guessed scale, the convergence of the conventional pQCD series shall generally change greatly under different choices of the renormalization scale, due to the mismatching of the perturbative coefficient with the α s -value at the same order. For example, there is quite large scale uncertainty for each term of the pQCD series of Γ(H → gg) [30,31]; thus, even if by choosing a proper scale, a better convergence can be achieved 4 , one cannot decide whether such a choice leads to the correct pQCD prediction. On the other hand, after applying the PMC, the scale-independent coupling α s (Q) can be determined, and together with the scale-invariant conformal coefficients, one can thus obtain the intrinsic perturbative nature of the pQCD series. By defining a K factor, K = Γ(H → gg)/Γ(H → gg)| Born , one can obtain the relative importance of the high-order terms to the leading-order terms. More explicitly, under the Landau gauge with ξ MOM = 0, we obtain  4 In the literature, the renormalization scale is usually chosen as the one so as to eliminate large logs with the purpose of improving the pQCD convergence; And some scale-setting approaches have been invented to find an optimal scale with the purpose of improve the pQCD convergence but not to solve the renormalization scale ambiguity.
Those results show satisfactorily convergent behavior for Γ(H → gg), especially for the mMOM, MOMh, MOMq and MOMg schemes. Similar to the condition of total decay width, the K factor is also gauge dependent. By varying ξ MOM ∈ [−1, 1], we obtain K mMOM = 1.07 +0.03 −0.10 , K MOMh = 1.07 +0. 12 −0.19 , K MOMq = 1.07 +0.12 −0.19 , K MOMg = 1.06 +0. 12 −0.19 and K MOMgg = 1.45 +0.00 −0.02 . As a final remark, one usually wants to know the magnitude of the "unknown" high-order pQCD corrections. The conventional error estimate obtained by varying the scale over a certain range is usually treated as such an estimation, which is however unreliable, since it only partly estimates the non-conformal contribution but not the conformal one. In contrast, after applying the PMC, the correct momentum flow of the process and hence the correct α s -value is fixed by the RGE and cannot be varied; otherwise, one will explicitly break the RGI, leading to an unreliable prediction. As a conservative estimation of the magnitude of the unknown perturbative contributions for the PMC series, it is helpful to use the magnitude of the last known term as the contribution of the unknown perturbative term [69]; As for the present case, we adopt ±|r 5,0 a 6 s (Q)| as the estimation of the unknown O(α 6 s ) contribution, which is ±3.1 KeV for the mMOM, MOMh, MOMq and MOMg schemes, and ±9.3 KeV for MOMgg scheme.

C. A simple discussion on the symmetric MOM schemes
In addition to the asymmetric MOM schemes, several symmetric MOM schemes have also been suggested in the literature. In the original symmetric MOM scheme [4], the triple-gluon vertex function Γ abc µνρ (k, p, l) is defined to be the value at the symmetric point k 2 = p 2 = l 2 = −µ 2 , i.e. the Feynman diagram of the vertex is the same as Fig.4 but should replace the three external momentum (q, −q, 0) there by the present (k, p, l). Similarly, the ghost-gluon vertex functionΓ abc µ (k, p, l) and the quarkgluon vertex function Λ a µ,ij (k, p, l) can also be defined at the symmetric point k 2 = p 2 = l 2 = −µ 2 [32,88,89]. For simplicity, we label those symmetric MOM schemes which are defined at the triple-gluon vertex, the ghostgluon vertex and the quark-gluon vertex as the MOMggg, the MOMh and the MOMq schemes, respectively. Using the relations (21,22), together with the known threeloop vertex functions under the MS scheme collected in Ref. [32], one can obtain the expressions for the strong couplings and the gauge parameters of those symmetric MOM schemes up to three-loop level [89]  We present the PMC effective scalesQ and their corresponding coupling constants α s (Q) of Γ(H → gg) under three symmetric MOM schemes in Fig.(12). The total decay width Γ(H → gg) versus the gauge parameter for the cases of three symmetric MOM schemes are presented in Fig.(13). Similar to the case of asymmetric MOM scheme, the gauge dependence cannot be suppressed by eliminating the scale dependence. More explicit, we present the total decay width Γ(H → gg) at several typical gauge parameters in Table. III.

V. SUMMARY
In the paper, we have made a detailed discussion on the gauge dependence of the total decay width Γ(H → gg) up to five-loop level under various MOM schemes. Our main results are: der the MOMgg scheme is the smallest and is less than ±1%. In this sense the MOMgg scheme could be treated as the best type of MOM scheme.
• By applying the PMC, the scale independent effective momentum flow (Q) of the process can be fixed by using the RGE, and as shown by Figs. (5,12), it differs for various MOM schemes, which ensures the scheme independence of pQCD predictions [40].
The MOMgg decay width has the smallest net error due to the small gauge dependence. It is found that the Higgs decay width Γ(H → gg) varies weakly on the choice of the MOM schemes, being consistent with the renormalization group invariance. Such small differences (less than ∼ 1%) among different schemes could be attributed to the unknown higher-order terms, e.g. the unknown N 4 LL and higher-order terms in the PMC scaleQ's perturbative series [71]. For example, with the help of Eq.(53), if treating ±|λ 3 a 3 (Q)| as an estimation of the contribution of the unknown N 4 LL-term ofQ, the change of the total decay width ∆Γ(H → gg) shall well explain the gaps of the total decay widths of different schemes, e.g. ∆Γ(H → gg)| PMC ≃ ±2.8 KeV for all the MOM schemes.
• The pQCD convergence of the conventional series varies greatly under different choices of the renormalization scale, due to the mismatching of the perturbative coefficient with the α s -value at the same order; Thus it is improper to use the conventional series to predict the unknown terms. On the other hand, after applying the PMC, the scaleindependent coupling α s (Q) is determined, and together with the scale-invariant conformal coefficients, one can achieve the intrinsic perturbative nature of the pQCD series and give more reliable prediction of unknown terms. Using the known five-loop prediction of Γ(H → gg) as an explicit example, if choosing ±|r 5,0 a 6 s (Q)| as a conservative estimation of the unknown six-loop prediction, we obtain the O(α 6 s )-order contribution is ±3. Eq. (18) shows that for any scheme R, we have and then we obtain On the other hand, Eq.(12) leads to In this Appendix, we give the perturbative transformations of the strong couplings and gauge parameters among a specific MOM scheme and the MS scheme. For convenience, we use the short notations (a, ξ) and (a, ξ) to represent (a MS , ξ MS ) and (a MOM , ξ MOM ), respectively. The transformations have firstly been considered in Ref. [4] and then been improved in Ref. [80]. Here for self-consistence and for our present needs, we give a more detailed derivation and one-order higher transformations than those of Ref. [80].
Generally, one can expand the strong couplings and gauge parameters under the MOM scheme over the ones under the MS scheme. Up to the present known order and to suit the needs of our present discussions, we have a = a + χ n (ξ)a n + O(a 4 ) .
Those two expansions (B3) and (B4) are important to transform the known MS perturbative series of a physical observable to the one under a certain MOM scheme. Our task is to derive the coefficients b i (ξ) and χ i (ξ) from the known ones φ i (ξ) and ψ i (ξ). For the purpose, we first do the following Taylor expansions: where ξ − ξ = ξ χ 1 (ξ)a + χ 2 (ξ)a 2 + χ 3 (ξ)a 3 + O(a 4 ) .