Matching at one loop for the four-quark operators in NRQCD

The matching coefficients for the four-quark operators in NRQCD (NRQED) are calculated at one loop using dimensional regularization for ultraviolet and infrared divergences. The matching for the electromagnetic current follows easily from our results. Both the unequal and equal mass cases are considered. The role played by the Coulomb infrared singularities is explained in detail.


Introduction
Effective field theories (EFTs) have become increasingly popular in describing processes where several scales are involved. In particular, two EFTs, namely Heavy Quark Effective Theory (HQET) and Non-Relativistic QCD (NRQCD) have been used for systems with heavy quarks. These EFTs take advantage of the fact that the masses of the heavy quarks (charm and bottom) are much larger than the remaining dynamical scales in the problem.
HQET was designed to study systems with one heavy quark [1,2,3] and has become a standard tool during the last years. Apart from the mass of the heavy quark (m) the remaining dynamical scales in heavy-light systems reduce to a single one Λ QCD . The HQET Lagrangian can be organized in a power series of the inverse pole mass of the heavy quark. Each term in this series consists of a gauge invariant operator. Only two kinds of terms turn out to be important for heavy-light systems: (i) terms containing light degrees of freedom (gluons and light quarks) only (which are irrelevant in most of the phenomenological applications), and (ii) terms containing a bilinear in the heavy quark fields. The size of each term is easily estimated by assigning the scale Λ QCD to whatever is not a heavy mass in the Lagrangian.
NRQCD was designed to study systems with a heavy quark and a heavy antiquark [4,5] and, although it is older than HQET, it has not received much attention until recently [6]- [11]. In this case, apart from the heavy quark mass, there are at least two dynamical scales. Namely the typical relative momentum in the bound state p and the typical binding energy E. Because of the existence of these two scales, the power counting rules are different from the HQET case and the size of each term in the NRQCD Lagrangian is not unique. Nevertheless, counting rules have been given to estimate the leading size of each term (see [5]). Independently of the relative size of each term, the NRQCD Lagrangian also consists of a power series of the inverse pole mass of the heavy quark. Here, though, there are important terms of three kinds: the two first kinds correspond exactly to (i) and (ii) in HQET, where we include in (ii) terms containing a bilinear in the antiquark fields as well. The third kind (iii) corresponds to operators bilinear in both heavy quark and heavy antiquark fields (four fermion terms).
A crucial step in building an EFT for heavy quarks is the so called matching. In the process of matching we enforce the effective theory to reproduce suitable S-matrix elements of the full theory. In this way we fix the parameters (Wilson coefficients) of the effective theory. Through the matching process the high energy contributions are encoded in Wilson coefficients multiplying the operators in the Lagrangian (and in the currents) of the effective theory. The determination of some of these Wilson coefficients of the NRQCD Lagrangian is the main topic of this paper.
The question arises whether the Wilson coefficients of the terms (i) and (ii) in HQET and NRQCD are the same. We shall support below the claim in [6] that this is indeed the case. Therefore, since the mass of the heavy quarks is (by definition) much larger than Λ QCD , the matching may be done order by order in 1/m and α s .
The matching for NRQCD has been known at tree level since long. This can be obtained by enforcing the tree level of S-matrix elements to be equal to those of QCD (QED) as mentioned above. For terms bilinear in the quark (antiquark) fields, this is equivalent to performing a Foldy-Wouthuysen transformation in the QCD Lagrangian. Although for HQET the matching at tree level for the bilinear terms can be carried out exactly as above [12], in most of the works it has been done somewhat differently: either by imposing the off-shell Green functions be equal to those of QCD [2] (see also [13]) or by integrating out the 'antiparticle' degrees of freedom [14]. The Lagrangian obtained in this way is in fact different from the NRQCD Lagrangian. However both Lagrangians are related by local field redefinitions or by using the equations of motion [6]. Results for the matching at one loop have also been known in the HQET for some time [3]. Nevertheless, attempts to perform the matching beyond tree level in NRQCD have not begun until recently. The main obstacle was that in NRQCD, unlike in HQET, the kinetic term was thought to be a necessary ingredient in the quark propagator for a matching calculation, If a hard cut-off is used (µ << m), it can easily be seen that the matching can be performed just like in HQET since k 0 >> k 2 /m in the ultraviolet. However, if dimensional regularization is used, the high energy modes (k > m) are not explicitly suppressed and they give non-vanishing contributions. This can be seen because the behavior of the NRQCD propagator changes at energies larger than the mass. In spite of this, one would like to use dimensional regularization because it keeps all the symmetries of QCD and, moreover, the calculations are technically simpler. Several authors have addressed this problem [7] and recently an appealing solution has been proposed [6]. There, it is claimed that the matching in NRQCD using dimensional regularization should be performed just like in HQET, namely the kinetic term must be treated as a perturbation. Let us make some remarks which support this approach. The key point is that in order to carry out the matching it is not so important to know the power counting of each term in the effective theory as to know that the dynamical scales of the effective theory are much lower than the mass. The power counting tells us the relative importance between different operators but this does not change the value of the matching coefficients. That is, we only need no matter what the relation between |p|, E and Λ QCD is. The above becomes clear if one thinks of the matching as a procedure to integrate out high energy degrees of freedom a la Wilson: the effective Lagrangian that we obtain after integrating out energies and momenta until a scale µ, m >> µ >> |p|, E, Λ QCD does not depend on the relative weight of the lower scales.
In addition, in ref. [6] dimensional regularization was used to regulate both the ultraviolet (UV) and the infrared (IR) divergences in the full and the effective theory [15]. The latter arise when the S-matrix elements are expanded about the residual momentum. In fact, it is not so important to know the way the UV divergences of the full theory are regulated since the comparison is done between S-matrix elements which are UV finite (after renormalization). Nevertheless, it is essential to regulate in the same way the IR divergences in both the full and effective theory in order for them to cancel out. This will always happen since by construction both theories have the same IR behavior. It is also important, from a practical point of view, to regulate the UV divergences of the effective theory using dimensional regularization. In this way, the calculation in the effective theory becomes trivial since there is no dimensionfull parameter in any integral. In the ref. [6], the matching was performed at one loop until O(1/m 2 ) for operators bilinear in the quark fields. It is the aim of this paper to perform the matching at one loop until O(1/m 2 ) for four-quark operators and hence to complete the matching at one loop O(1/m 2 ).
We are thus faced with the computation of S-matrix elements of four heavy quarks in QCD and HQET. The computation of these matrix elements in HQET is unusual, although some related calculations already exists in the literature [16,17]. Indeed, for heavy-light systems four fermion operators are relevant only when two of the fermions are light. For heavy quarkonium systems instead all four quark fields are heavy. In fact it is in these S-matrix elements where we can see the peculiar IR behavior of heavy-heavy systems, which eventually gives rise to the Coulomb pole and hence to the standard non-relativistic weak coupling bound states.
This IR behavior should appear in both the full and the effective theory. If we expand about the residual momentum the matrix elements of the dimensionally regulated QCD, we may expect an IR singularity reflecting the Coulomb pole. However this singularity corresponds to an odd power-like IR divergence and hence it is put to zero in dimensional regularization. This is not a problem. Indeed, since the effective theory has the same IR behavior, it also has an IR divergence reflecting the Coulomb pole which is consistently put to zero by dimensional regularization. The important thing when doing the matching is to take into account all the non-analytical behavior in the heavy quark masses which can not be obtained in the effective theory. Proceeding in this way we are certainly taking into account all the non-analytical behavior in the masses coming from high momenta (QCD logs). The remaining non-analytical behavior (Coulomb pole) is encoded in the effective theory.
Although we have been talking about QCD and NRQCD, the results for QED and NRQED follow trivially from our calculations.
Let us finally mention some of the possible applications of this work. The unequal mass case may be important for the B c system (this system has been studied in refs. [18]) which is expected to be seen in the future. This case is also important in QED for the muonium or Hydrogen-like atoms. For the equal mass case, our results fixes the scale of the α s running constant for annihilation contributions to the four quark interaction. This is important since in QCD, at the scales of Bottomonium and Charmonium, α s strongly depends on the scale. Moreover, since there are many scales in the game (m, p, E) it is not a priori clear which one should be used in order to fix the value of α s in the perturbative [19,20] and non-perturbative potentials [21]. In fact, depending on where the contribution comes from, this value may be different. Recently, the spectrum of Υ(1s) and J/ψ has been obtained from perturbative QCD at O(mα 4 s ) [20]. Next improvement, namely O(mα 5 s ) receives contributions from the matching of four quark operators, and hence our calculation becomes relevant. It should also be taken into account in parameterizations of the non-perturbative heavy quark potential along the lines of reference [21].
We distribute the paper as follows. In sec. 2 we define our four quark operators and their Wilson coefficients. In sec. 3 we calculate the Wilson coefficients for the unequal mass case. In sec. 4 we calculate the Wilson coefficients for equal mass case. In sec. 5 we discuss a few relevant issues in our calculation. The last section is devoted to the conclusions. A few technical points concerning the Coulomb singularity are relegated to an Appendix.

Setting the matching
The piece of the NRQCD Lagrangian containing four quark operators at O(1/m 2 ) reads where ψ is the Pauli spinor field that annihilates a heavy quark and χ is the Pauli spinor field that creates a heavy anti-quark. The subindices 1,2 denotes the possibility of working with different particles (different masses). We will omit these indices when the particleantiparticle case is treated.
There is another possibility of writing down these terms by using Fiertz transformations. It reads The relation between the two bases is Of course, one can always use a redundant bases with the eight operators (2.1) and (2.2). The Lagrangian (2.2) is more convenient, as far as the matching calculation is concerned, when one is dealing with annihilation processes in the equal mass case. Nevertheless (2.1) is a better option when one addresses a bound state calculation. We shall use (2.1) for the unequal mass case and the redundant basis for the equal mass one in order to ease comparison with the actual calculations.
In the QED case we have Now, the relation between the two bases is We shall expand the dimensionally regulated matrix elements about zero residual momentum. Since there are no derivative terms in (2.1) and (2.2), the zeroth order in the expansion will be enough. Namely we only have to calculate the matrix element for the four quarks at rest. This means that the amputated legs in a digram only have to be multiplied either by p + (projector on the particle subspace) or p − (projector on the antiparticle subspace), and the kinematic factor m/E relating relativistic and nonrelativistic normalizations can be put to one. We shall use MS subtraction scheme for both UV and IR divergences and work with the Feynman gauge. The matching coefficients should be gauge independent, but they depend on the subtraction scheme. It is worth emphasizing that we do not work in the on-shell renormalization scheme for the wave function (of course our masses always correspond to the pole mass), but just MS. In this scheme (also in MS or similars) the matching can be carried out straightforwardly. If the on-shell scheme is used for the full or effective theory, one must identify the UV divergences which correspond to a wave function renormalization and subtract them accordingly (not just minimally). This is obviously more tedious than using MS throughout. The little price to be paid for this simplicity is that our fields are not properly normalized. This must be taken into account by including the proper Z factors when calculating on-shell (2.7) Notice finally that at the order we are working at the Wilson coefficients in (2.1) and (2.2) are invariant under the local field redefinitions discussed in [6].

Unequal mass case
In this case annihilation diagrams are forbidden and we are only left with the two QCD diagrams of fig. 1. In these diagrams the Coulomb singularity can be identified and the above mentioned mechanism by which it disappears uncovered. We show this in detail in the Appendix. In short, things go as follow. In order to perform some integrals we have to move to dimensions high enough in order to regulate the IR Coulomb singularity. When coming back to four dimensions we can trace back the IR Coulomb singularity as a pole in higher dimensions, which does not appear in four dimensions since dimensional regularization loses power-like divergences. The point is that we have not provided a suitable dimensionfull parameter (the relative momentum) and hence dimensional regularization has no way to reproduce the Coulomb pole. This fact was pointed out some time ago in ref. [17].
We obtain the following matching coefficients In the case of QED our results reduce to The spin dependent piece, which is subtraction point independent, agrees with the result obtained by Caswell and Lepage in ref. [4]. The scalar piece is new.
Since we will need the equal mass results in the next section let us display them here. For QCD they read For QED we have Recently, the scalar piece for the equal mass case in QED was calculated in ref.
[8] using a cut-off regularization. This result agrees with ours except for a finite piece, which may be due to a different renormalization scheme for the four fermion operators in the effective theory. This contribution is relevant for the full calculation of the positronium energy levels at order O(mα 5 ). Work in this direction is under way [22].

Equal mass case
For equal particles annihilation processes are allowed and they should be taken into account ( fig. 2). From fig. 2a we obtain the well known result [5]  This is the lower order (tree level) contribution. Let us consider first the one loop contributions arising from the gluon self-energy. Each heavy quark loop ( fig. 2b) gives a contribution where T R = 1/2 for QCD and T R = 1 for QED. The QED result had been already obtained in refs. [8,9]. Light quarks (n f ) and gluons give a contribution ( fig. 2c) Let us next consider the vertex corrections (figs. 2e and 2f). Fig. 2e is quite interesting. A singularity associated with the Coulomb pole should appear, but again it does not show up when doing the computation in dimensional regularization for quarks at rest. This is totally analogous to what happened with diagrams in fig. 1 in the previous section. No signal of infinity imaginary anomalous dimension appears either [16]. We refer the reader to ref. [17] for a detailed explanation on what is going on. We obtain For QED fig. 2e has been already computed [8,9].
The diagrams of fig. 2f do not appear in QED. They lead to iπ . (4.6) Finally we consider the contributions from the diagrams of fig. 2g. These diagrams also exist in QED and their contributions in this theory have already been calculated in ref. [8]. We obtain from the diagrams of fig. 2g d c,2g and d c,2g sv is zero. For QED we reproduce the results in [8,9]. Summarizing all the contributions from annihilation diagrams we obtain Recall that we have to add to the annihilation contributions above the contributions (3.7)-(3.10).

Discussion
Let us first address the important question on how the matching calculation helps to fix the scale of α s . In the previous section we have not paid any attention to the flavor dependence of α s . For simplicity let us focus on the case of a single heavy flavor. Suppose that in QCD we have N f flavors. Then in NRQCD we have n f = N f −1 relativistic flavors. Consequently, the MS running coupling constant in NRQCD is expected to run according to N f − 1 flavors. However, this is not obvious from the matching calculation. Notice that the α s s in the NRQCD Lagrangian (both explicit and in the Wilson coefficients) are those inherited from QCD and hence one may be tempted to make them run with N f flavors. In order to clarify this issue consider first the pure gluonic part of the NRQCD Lagrangian [6].
Notice that the kinetic term does not have the standard normalization anymore. This can be recovered by a simple redefinition of the gluon field. Since the remaining gluon fields in the NRQCD Lagrangian are multiplied by g, this is equivalent to make the change in all the gs which multiply the gluon fields. At one loop this is nothing but changing the running coupling constant of N f flavors for the running coupling constant of N f −1 flavors which is a desired result. However there are additional dependences on α s in the NRQCD Lagrangian (which are not multiplying gluon fields) in the Wilson coefficients. Notice however that the difference between α  (4.13). Indeed, the Wilson coefficient in (4.13) is linear in α s and hence it may be sensitive on whether this α s corresponds to N f or N f − 1 flavors. Since this α s is inherited from QCD it corresponds to N f flavors. However the ν dependence enters in such a way that and hence the scale of α s is naturally fixed to m. The Wilson coefficient d c vv in (4.13) should better be written like Nevertheless, if there is no dynamical scale between m 1 and m 2 and we are not interested in any renormalization group improvement of the Wilson coefficients, carrying out the matching in one step or in two steps must give exactly the same result to any fixed order in perturbation theory. Then, we expect our results to be useful for the B c meson in QCD and for the muonium and hydrogen-like atoms in QED. Recall that the Wilson coefficients in (5.2) trivially change into Now rescaling the gluon field to its usual normalization (d 1 = 1) is equivalent to in the coupling constants multiplying the gluon fields, and hence these α s run with N f − 2 flavors. Notice also that the following equality holds, α For QED an analogous discussion implies that (4.15) is given in terms of the QED running coupling constant. This expression in terms of the low energy α (α ∼ 1/137) reads If we add (3.11) and (3.12) to (4.14) and (5.8) respectively, we obtain the results presented in [10].
Imaginary parts appear in the Wilson coefficients at several instances. In order to obtain them from our expressions beware that we have located the cut at the negative real axes of the m 2 complex plain. These imaginary parts have to do with inelastic cross sections which cannot be obtained within the non-relativistic theory alone. They are also related to the decay width of heavy quarkonium states into light hadrons and had been calculated before [4]. Our results agree with this previous calculation.
A word of caution is required when dealing with the Pauli matrices in D dimensions. The Pauli matrices arising in NRQCD have in fact very different origin, as we comment next. For the non-annihilation diagrams the Pauli matrices originate from While the first equality can be understood as a definition, for the second one we have used the following prescriptions (with the proper limit when D → 4) The finite part of d vv depends on these prescriptions. For the annihilation diagrams the Pauli matrices originate from 1 When carrying out calculations in dimensionally regulated NRQCD the same prescriptions have to be used and it may be eventually important to keep in mind the different origin of the various σ k ⊗ σ k and 1 ⊗ 1. For NRQED there are no ambiguities at this order since the spin dependent terms are finite.
Let us finally mention that the matching at one loop for the electromagnetic current at leading order in 1/m arises trivially from the calculation of the diagrams in figs. 2b and 2c. Schematically, the full current is approximated by where δ r encodes the one loop correction due to hard gluons. Now, one only has to realize that the relevant computation (the matching procedure follows analogously to the one for the four-fermion operators) is the one we performed for the diagrams above but taking into account the different color factors. We obtain which agrees with the well known result. We consider this procedure by far the simplest and most efficient method to obtain δ r (one can also trivially obtain the result for QED, δ r = −4α/π). No problem with the Coulomb pole appears through the calculation. Notice also that no anomalous dimension appears either. This could be traced back to the fact that in both QCD and HQET for one quark and one antiquark (the effective theory to which we are matching to from a practical point of view) have symmetries which protect this current. For the effective theory this symmetry is U(4) [23].

Conclusions
We have calculated the matching coefficients of the four quark operators of NRQCD at one loop and O(1/m 2 ). We have considered both the unequal and equal mass cases. We have shown explicitly how some matching coefficients in the NRQCD Lagrangian conspire in such a way that all α s appearing in them must be considered at the scale m whereas the α s multiplying the gluon fields must be considered as running with the number of remaining relativistic flavors only.
The binding energies of Υ(1s) and J/ψ have been recently obtained from perturbative QCD at O(mα 4 s ) [20]. Next improvement, namely O(mα 5 s ) in the NRQCD framework requires the knowledge of the matching coefficient of the four quark operators at one loop calculated here. In the framework of NRQED, these are also necessary to obtain the positronium binding energy at O(mα 5 ).
The unequal mass case in NRQCD may have eventual applications to the B c meson. For NRQED it may be relevant for precision calculations (involving recoil corrections) in muonium and hydrogen-like atoms. In particular it would be relevant for the spectrum of an hydrogen atom at O(mα 5 ) where the electron has been substituted by a τ particle.
Acknowledgments A.P. acknowledges a fellowship from Generalitat de Catalunya. Financial support from CICYT, contract AEN95-0590 and financial support from CIRIT, contract GRQ93-1047 is also acknowledged.

A Coulomb singularity
In this appendix we show how the Coulomb singularity is reflected in our calculation. Upon integration over q 0 we obtain IR singularities from the poles of the quark propagators and from the poles of the gluon propagators.

Consider the following integrals
The poles of the quark propagators produce the Coulomb singularity where Λ → 0 is an IR cut-off. At D=4 this integral has odd power like singularities which are ignored by dimensional regularization. However we expect these singularities to shown up as poles in an odd number of dimensions.
The poles in the gluon propagators also give rise to IR singularities. These read For n = 1 and n = 2 we expect a pole in D = 4 and D = 6 respectively (an extra pole at D = 4 cannot be ruled out a priori for n = 2 but it will not turn up).
The explicit result for I n below fulfills the expectations above where D = 4 + 2ǫ. Notice that the ǫ = n singularities are of UV origin.