Three-loop QCD matching of the flavor-changing scalar current involving the heavy charm and bottom quark

We compute the matching coefficient between the quantum chromodynamics (QCD) and the non-relativistic QCD (NRQCD) for the flavor-changing scalar current involving the heavy charm and bottom quark, up to the three-loop order within the NRQCD factorization. For the first time, we obtain the analytical expressions for the three-loop renormalization constant Z~s(x,Rf)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{Z}_s(x,R_f)$$\end{document} and the corresponding anomalous dimension γ~s(x,Rf)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\gamma }_s(x,R_f)$$\end{document} for the NRQCD scalar current with the two heavy bottom and charm quark. We present the precise numerical results for those relevant coefficients (CFF(x0),…,CFBB(x0))\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(C_{FF}(x_0), \ldots , C_{FBB}(x_0))$$\end{document} with an accuracy of about thirty digits. The three-loop QCD correction turns out to be significantly large. The obtained matching coefficient Cs(μf,μ,mb,mc)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_s(\mu _f,\mu ,m_b,m_c)$$\end{document} is helpful to analyze the threshold behaviours when two different heavy quarks are close to each other and form the double heavy Bc\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$B_c$$\end{document} meson family.


I. INTRODUCTION
Heavy quark system provides a unique window to understand the perturbative and nonperturbative nature of Quantum Chromodynamics (QCD) theory.The key quantity of the heavy quark system is the heavy quark mass with m Q Λ QCD , which is naturally employed to distinguish the perturbative and nonperturbative interactions.For the two heavy quark system such as heavy quarkonium and the threshold production of top quark pair, the nonrelativistic QCD (NRQCD) effective theory is a very successful approach to separate the long-distance nonperturbative interactions and the short-distance perturbative interactions [1].In this effective theory, the QCD observable can be further expanded into the NRQCD effective operator matrix elements with the corresponding matching coefficients.The matching coefficients can be calculated order by order and can be in series of two small parameters, i.e. the strong coupling constant α s and the quark relative velocity v.
Up to now, various matching coefficients for the two heavy quark currents, involving the b b, bc and cc systems, have been computed to higher-order accuracy within the NRQCD effective theory.Taking the double-heavy B + c = bc meson as an example, the next-to-leading order (NLO) matching coefficient for the axial-vector current was first obtained in 1995 [2].Then, the NLO matching coefficient for the vector current was calculated in 1999 [3].Later on, the NLO QCD corrections combined with the higher-order relativistic corrections for the axial-vector current and the vector current are investigated in Ref. [4].The approximate results of the next-to-next-to-leading order (NNLO) QCD correction for the axial-vector current matching coefficient was first obtained in 2003 [5].However, the complete analytical expression of the two-loop matching coefficient for the axial-vector current became available in 2015 [6].This year, the next-to-next-to-next-to-leading order (N 3 LO) QCD correction for the matching coefficient of the axial-vector current was first numerically calculated in Ref. [7].Very recently, the two-loop matching coefficient of the vector current was evaluated by the authors of this work [8].After that, the three-loop matching coefficient of the vector current was also achieved in Ref. [9].
However, the matching coefficient for the flavor-changing scalar current involving the heavy bottom and charm quark has not yet been found in previous works.Considering the state-of-the-art theoretical accuracy of the matching coefficients for other heavy quark currents, we will compute the matching coefficient for the flavor-changing heavy quark scalar current of the B c meson up to the three-loop QCD corrections within the NRQCD factorization frame.The novel calculation is helpful to analyze the threshold behaviours when two different heavy quarks such as the bottom and the charm quark are close to each other.In addition, the three-loop QCD matching procedure for the heavy flavor-changing scalar current of the B c meson provides a window to check the convergence of the perturbative series of the NRQCD effective theory.
The rest of the paper is arranged as follows.In Sec.II, we introduce the matching formula between the full QCD theory and the NRQCD effective theory.We present the analytical results of the three-loop scalar current renormalization constant and the corresponding anomalous dimension in the NRQCD effective theory.We also discuss the matching and running equations for the strong coupling constant α s .In Sec.III, we give our calculation procedure and the projection for the heavy flavor-changing scalar current.In Sec.IV, we present our numeric results of the matching coefficient up to the three-loop order accuracy.Finally, we summarize in Sec.V.

II. MATCHING FORMULA
The heavy flavor-changing scalar current involving b and c quark in the full QCD can be written as j s (y) = Ψb (y)Ψ c (y), which can be further expanded into the NRQCD effective scalar current js with the reduced heavy quark mass m red = m b m c /(m b + m c ) and the quark relative momentum p at the leading order of quark relative velocity, according to the definition as given in Ref. [10].ψ is the two-component Pauli spinor field that annihilates a heavy quark, while χ is the two-component Pauli spinor field that creates a heavy anitiquark.
The matching coefficient can be determined through the conventional perturbative matching procedure.Namely, one performs renormalization for the on-shell vertex functions in both the perturbative QCD and the NRQCD sides, then solves the matching coefficient order by order in α s .
The matching formula with renormalization procedure reads where the left hand part of the equation represents the renormalization of the full QCD current while the right hand part represents the renormalization of the NRQCD current.The term O(v 2 ) denotes the higher order relativistic corrections in powers of the heavy quark relative velocity v between the bottom quark b and the charm quark c.Γ 0 s ( Γ 0 s ) denotes the on-shell unrenormalized heavy flavor-changing current vertex function in the QCD (NRQCD) theory, and the leading-order (LO) matching coefficient is normalized into 1.In this paper, we will consider higher-order QCD corrections up to O(α 3 s ) but at the lowest order in heavy quark relative velocity v [19,29].Z s is the QCD current renormalization constant in on-shell (OS) scheme, i.e., Z s = . And Zs is the renormalization constant of the NRQCD effective current in the modified-minimal-subtraction (MS) scheme.Z 2 and Z m are QCD on-shell quark field and mass renormalization constants, respectively.The three-loop analytical results of the on-shell quark field and mass renormalization constants allowing for two different non-zero quark masses can be found in literature [30][31][32], which can be evaluated to high numerical precision with the package PolyLogTools [33].The QCD coupling MS renormalization constant can be found in literature [34][35][36].The NRQCD on-shell quark field renormalization constants Z2,b = Z2,c = 1, since all light particles in NRQCD are massless.The matching coefficient C s (µ f , µ, m b , m c ) depends on the NRQCD factorization scale µ f and the QCD renormalization scale µ in a finite order QCD correction calculation.
After implementing the quark field, quark mass and the QCD coupling constant renormalization, the QCD vertex function gets rid of ultra-violet(UV) poles, while still contains uncancelled infra-red(IR) poles starting from order α 2 s .The remaining IR poles in QCD should be exactly cancelled by the UV divergences of Z s in NRQCD, which renders the matching coefficient finite.With the aid of the obtained high-precision numerical results, combined with the features of the NRQCD current renormalization constants investigated in other known literature [7,9,12,29], we have successfully reconstructed the exact analytical expression of the NRQCD renormalization constant for the flavor-changing heavy quark scalar current through numerical fitting recipes [33,37,38].Here we directly present the final result as following where the coefficients Z s (x) and Z (3) s x, are of the following form 34 + 72 ln 2 − 9 ln x + 18 ln(x + 1) + 9 ln where The corresponding anomalous dimension γs for the NRQCD scalar current is related to Zs by [22,[39][40][41][42][43] where Z [1] s denotes the coefficient of the 1 pole in Zs , and the NRQCD factorization scale µ f is used because both Zs and γs are defined in the NRQCD effective theory.The explicit expressions of the coefficients γ(2) s (x) and γ(3) in Eq. ( 6) are of the form of From above expressions, one can see that both of Zs and γs explicitly depend on the NRQCD factorization scale µ f , which is a feature found in other NRQCD currents [7,9,12,19,29].Note that the above three-loop expressions of Z s and γs for the flavor-changing scalar current involving the heavy charm and bottom quark are known for the first time.The obtained Z s and γs have been checked with several different values of m b and m c .To verify the correction of our results, on the one hand, one can check that above Z s (γ s ) is symmetric under the combined exchange of m b ↔ m c and n b ↔ n c .On the other hand, in the equal quark mass case of x = 1, our Z s and γs are in full agreement with the known results as given in Refs.[10][11][12].
In our calculation, we include the contributions from the loops of charm quark and bottom quark in the full QCD, which however are decoupled in the NRQCD.To match the QCD with the NRQCD, one need apply the decoupling relation [13,[43][44][45][46] of α s , i.e., the coupling constants α (n l +1) s (µ) in QCD (with n l + 1 flavours) and α (n l ) s (µ) in the NRQCD (with n l light flavours) are related by [13] α where L = ln(µ 2 /m 2 Q ) and mQ is the on-shell mass of the decoupled heavy quark.
Besides, we can evolve the strong coupling from the scale µ f to the scale µ with renormalization group running equation [47] in D = 4 − 2 dimensions as following To calculate the values of the strong coupling, we also use the renormalization group running equation [48] in D = 4 dimensions as where . At the one-, two-and three-loop level, the coefficients of the QCD β function are of the form of In our numerical evaluations, n b = nc = 1 and n l = 3 are fixed through the decoupling region from µ = 1 GeV to µ = 6.25 GeV, and the typical QCD scale Λ (n l =3) QCD = 0.3344GeV is determined using the three-loop formula with the aid of the package RunDec [48][49][50][51] by inputting the initial value α (n f =5) s (mZ = 91.1876GeV)= 0.1179.

III. CALCULATION PROCEDURE
Our high-order calculation consists of the following steps.First, we use FeynCalc [52] to obtain Feynman diagrams and corresponding Feynman amplitudes.By $Apart [53], we decompose every Feyman amplitude into several Feynman integral families.Second, we use FIRE [54]/Kira [55]/ FiniteFlow [56] based on Integration by Parts (IBP) [57] to reduce every Feynman integral family to master integral family.Third, based on symmetry among different integral families and using Kira+FIRE+Mathematica code, we can realize integral reduction among different integral families, and further on, the reduction from all of master integral families to the minimal set [58] of master integral families.Last, we use AMFlow [59], which is a proof-of-concept implementation of the auxiliary mass flow method [60], equipped with Kira [55]/FiniteFlow [56] to calculate the minimal set of master integral families.
In order to obtain the finite results of the high-order QCD corrections, one has to perform the conventional renormalization procedure [6,10,61,62].Equivalently, we can also use diagrammatic renormalization method [63] with the aid of the package FeynCalc [52], which at N 3 LO sums contributions from three-loop diagrams and four kinds of counter-term diagrams, i.e., tree diagram inserted with one α 3 s -order counter-term vertex, one-loop diagram inserted with one α 2 s -order counter-term vertex, one-loop diagram inserted with two αs-order counter-term vertexes, two-loop diagrams inserted with one αs-order counter-term vertex.Our final finite results by these two renormalization methods are in agreement with each other.
We want to mention that all contributions up to NNLO have been evaluated for general gauge parameter ξ and the NNLO results for the scalar current matching coefficient are all independent of ξ, which constitutes an important check on our calculation.At N 3 LO, we work in Feynman gauge.By FeynCalc, there are 1, 1, 13, 268 bare Feynman diagrams for the QCD vertex function with the flavor-changing heavy quark scalar current at tree, one-loop, two-loop, three-loop orders in αs, respectively.Some representative Feynman diagrams up to three loops are displayed in Fig. 1 and Fig. 2. In the calculation of multi-loop diagrams, we have allowed for n b bottom quarks with mass m b , nc charm quarks with mass mc and n l massless quarks appearing in the quark loop.Physically, n b = nc = 1 and n l = 3.To facilitate our calculation, we take full advantage of computing numerically.Namely, before generating amplitudes, m b and mc are chosen to be particular rational number values [64][65][66].Following the literature [10], we employ the projector constructed for the flavor-changing heavy quark scalar current to obtain intended QCD amplitudes, which means one need extend the scalar current projector with equal heavy quark masses in Eq. ( 8) of Ref. [10] to the different heavy quark masses case.Adopting the same notation of Ref. [10], we choose q1 = mc m b +mc q + p and q2 = m b m b +mc q − p denoting the on-shell charm and bottom momentum, respectively, and present the projector for the flavor-changing heavy quark scalar current as where the small momentum p refers to relative movement between the bottom and charm, q represents the total momentum of the bottom and charm, q To match with the NRQCD, one need extract the contribution from the hard region in the full QCD amplitudes for the scalar current, which means one need first introduce the small relative momentum p to momenta in the amplitudes as above and then series expand propagator denominators with respect to p up to O(p) in the hard region of loop momenta [10].As a result, the number and powers of propagators in Feynman integrals constituting the amplitudes for the scalar current will remarkably increase compared with the vector current case, the zeroth component of the axial-vector current and the pseudoscalar current case.In our practice, the total number of propagators in a three-loop Feynman integral family is 12.In our calculation, the most difficult thing is the reduction from Feynman integrals with rank 5, dot 4, and 12 propagators to the master integrals.By trial and error, we find it is more appropriate for Fire6 [54] to deal with this problem than Kira [55] or FiniteFlow [56].After using Kira+FIRE+Mathematica code to achieve the minimal set of the master integral families based on symmetry among different integral families, the number of three-loop master integral families is reduced from 829 to 26, meanwhile the number of three-loop master integrals is reduced from 13251 to 300.
b c FIG. 1.Typical Feynman diagrams for the QCD vertex function with the bc system up to two-loop order.The cross " " implies the insertion of the flavor-changing heavy quark scalar current.

FIG. 2.
Typical Feynman diagrams labelled with the corresponding color factor for the QCD vertex function with the flavor-changing heavy quark scalar current at three-loop order.The cross " " implies the insertion of the scalar current.The thickest solid closed circle represents the bottom quark loop, and the other solid closed circle represents the charm quark loop.The dotted closed circle represents the ghost loop.

IV. RESULTS
Following Refs.[7,9,29], the dimensionless matching coefficient Cs for the flavor-changing scalar current involving the heavy bottom and charm quark can be decomposed as: where the coefficients γ(2) s (x) and γ(3) have been defined in Eqs.(6,7,8).The parameters C (i) (x)(i = 1, 2, 3) in above equation are independent of ln µ and ln µ f and are the nontrivial parts of Cs at O α i s .It's also well known that the coefficients Cs and C (i) (x) satisfy the following symmetric replacements [2][3][4][5][6][7]9]: The one-loop QCD correction to Cs, denoted by C (1) (x), can be analytically achieved as: The two-loop and three-loop matching coefficients in ( 14) are C (2) (x) and C (3) (x), which can be decomposed in terms of the different color/flavor structures by following the conventions as being used in Refs.[7,9,12,19,29,67]: Due to limited computing resources, we choose to calculate the matching coefficient Cs at three rational GeV, mc = 475 100 GeV} (i.e., x = 1), respectively.The results Cs obtained at the physical point and the check point verify the symmetric features of C (i) (x) in Eq. ( 16).Our results Cs obtained at the equal mass point x = 1 are consistent with the known matching coefficient results Cs for the scalar current with the equal quark mass case in the previous literatures [10][11][12].To confirm our calculation, we have also applied the same calculation procedure to the evaluation of the three-loop matching coefficients Cv, Cp, C (a,0) for the flavor-changing heavy quark vector current, the flavor-changing heavy quark pseudoscalar current, the zeroth component of the axial-vector current with the flavor-changing heavy quarks, respectively, where our results verify Cp ≡ C (a,0) and our results Cv and C (a,0) agree with the known results in previous literature [7][8][9][10][11][12]19].
In the following, we will present the highly accurate numerical results of C (2) (x) and C (3) (x) at the physical heavy quark mass ratio x = x0 = 1.5/4.75 with about 30-digit precision.The various components of C (2) From the numerical results in Eqs.(20,21), we find the dominant contributions in C (2) (x0) and C (3) (x0) come from the components corresponding to the color structures C 2 F , CF CA, C 2 F CA and CF C 2 A , and the contributions from the bottom and charm quark loops are negligible.
Fixing the renormalization scale µ = m b = 4.75GeV, mc = 1.5GeV, and setting the factorization scale µ f = 1 GeV, Eq. ( 14) then reduces to With the values of α From Fig. 3 and Tab.I, one can find that the higher order QCD corrections have large influence compared with the NLO correction.Especially, the O(α 3 s ) correction looks quite sizable, which confirms the breakdown of the convergence at O(α 3 s ) for other heavy quark currents in previous literatures [7,9,12].Note that, at each truncated perturbative order, the matching coefficient Cs is renormalization-group invariant [7,9], e.g., at N 3 LO, Cs obeys the following renormalization-group running invariance: where C N 3 LO s (µ f , µ, m b , mc) has dropped the O(α 4 s ) terms in Eq. ( 14).Namely the scale-dependence of the N 3 LO matching coefficients rises from O(α 4 s ), which is suppressed compared to lower-order.However, the scale-independence coefficients of α 4 s such as C (3) (x) and ln µ f , which come from the O(α 3 s ) order in Eq. ( 14), are considerably large by aforementioned calculation within the framework of the NRQCD theory.These terms will lead to a significantly larger renormalization scale dependence at N 3 LO.From Fig. 3, we also find the NRQCD factorization scale µ f has a large influence on the matching coefficient.When µ f decreases, both the convergence of αs expansion series and the independence of µ will be improved.The understanding of the observed large NNNLO corrections in this paper and the previous literatures [7,9,12] is another important topic, which may be related to the higher order corrections to NRQCD long-distance nonperturbative matrix elements, higher order relativistic corrections and the resummation techniques.We will leave it in the future work.

V. SUMMARY
In this paper, we have performed the N 3 LO QCD corrections to the matching coefficient for the flavor-changing heavy quark scalar current involving the bottom and charm quark within the framework of the NRQCD effective theory.For the first time, we obtain the analytical expressions of the three-loop renormalization constant and the corresponding three-loop anomalous dimension of the NRQCD scalar current with different heavy quark mass m b and mc.Meanwhile, the three-loop matching coefficient Cs(µ f , µ, m b , mc) has also been obtained with high numerical accuracy, which is helpful to analyze the threshold behaviours when two different heavy quarks are close to each other.The obtained N 3 LO QCD corrections are considerably large compared with lower-order corrections, and exhibit stronger dependence on the QCD renormalization scale and the NRQCD factorization scale at higher order, which suggests that higher order QCD corrections should be combined with other techniques in order to get a reliable prediction within the NRQCD effective theory.
and the parameter x = m c /m b representing the ratio of the heavy charm and bottom quark mass.

(n l FIG. 3 .
FIG. 3. The renormalization scale dependence of the matching coefficient Cs at the LO, NLO, NNLO and N 3 LO accuracy.The central values of the matching coefficient Cs are calculated inputting the physical values with µ f = 1 GeV, m b = 4.75GeV and mc = 1.5GeV.The error bands come from varying µ f from 1.5 to 0.4 GeV.