Manifestly Soft Gauge Invariant Formulation of vNRQCD

Homogeneous power counting in Non-Relativistic QCD (NRQCD) implies the simultaneous existence of both soft and ultrasoft gluons. In the velocity renormalization group (vNRQCD) formalism we show that operators involving soft fields interacting with heavy potential quarks can be put in a manifestly gauge invariant form by utilizing gluon and quark building blocks which contain Wilson lines, and are analogous to those used in the soft collinear effective theory. This leads to several simplifications, in particular significantly reducing the size of the operator basis, which stream-lines matching and anomalous dimension calculations at subleading order in the velocity expansion. Also, soft ghosts no longer couple via potential like interactions to the heavy quark fields, and hence do not appear until two loops. Furthermore the the color structures that appear at each order in $\alpha_s$ in the static potential are only those which should arise according to non-Abelian exponentiation. We also discuss how zero-bin subtractions clarify the role of $i\varepsilon$ terms in the heavy quark propagator poles carrying soft momenta. Even though the choice of the direction of the soft Wilson lines in our building blocks does not matter, there is still a limited set of consistent choices, and the predicted form of the eikonal poles does influence results at two-loops and beyond.


Introduction
Non Relativistic QCD (NRQCD) is an effective field theory (EFT), proposed in its original form in [1,2], which has been successfully used to describe the dynamics of non-relativistic heavy quark-antiquark bound states. This EFT includes three scales which are the heavy quark mass m, momentum mv and energy mv 2 . What complicates NRQCD compared to more canonical effective field theories is that two infrared scales, i.e energy and momentum, are not independent.
They are instead related by the free dispersion relation, a fact encoded by carrying out the power expansion in the heavy quark relative velocity, and correspondingly the logarithms of p and E are in general not independent. A systematic way to deal with these coupled scales was proposed in [3], see [4,5] for reviews and [6][7][8][9][10][11] for earlier foundational work. In this formulation, usually called "vNRQCD", one carries out the matching in a single stage, and simultaneously sum logs of both the energy and momentum by using what is called the "velocity renormalization group", as has been carried out in Refs. [12][13][14][15][16][17][18][19][20]. The framework of pNRQCD provides an alternative formulation of the EFT where the soft and ultrasoft scales are treated independently-for more of the setup, see Refs. [9,[21][22][23][24][25][26][27][28].
While the present formulation of vNRQCD in [3] is quite useful for summing logs, and yields correct results, the operators involving the interaction of soft gauge fields and heavy quarks are not manifestly soft gauge invariant. This shortcoming of the formalism makes computing the soft contribution to the running of the subleading potentials cumbersome, as can be seen in Refs. [12,13]. Although the final results for these potentials are gauge invariant, this lack of manifest symmetry proliferates the number of integrals that need to be calculated. It should be expected that if we instead formulate the theory in a manifestly gauge invariant way (in both the soft and ultrasoft sectors), then simplifications will follow.
In this paper, following ideas put forward in [29] in the context of forward scattering of energetic particles in the soft collinear effective theory [30][31][32][33], we formulate NRQCD so as to maintain soft gauge invariance at all stages of the calculation. This is accomplished by including soft Wilson lines in operators that connect soft and potential fields in the action through gauge invariant soft gluon and quark building block fields. We demonstrate the utility of this formalism by showing how it simplifies the calculation of the anomalous dimensions of the heavy quark potential to order O(v 2 ). We also show how our formalism, in conjunction with the vNRQCD zero-bin subtractions [34], clarifies the issue regarding the known importance of the iε prescription in obtaining the proper result for the two-loop Coulomb potential [35,36].

Review of NRQCD with Single Stage Matching
In vNRQCD we match onto the EFT at m, the mass scale of the heavy quarks. In the low energy theory there are two relevant scales mv (soft) and mv 2 (ultrasoft), where v is a scaling parameter that is of order the relative velocity between the heavy particles. The heavy quarks and antiquarks whose dynamics we wish to describe have kinetic energy p 0 ∼ mv 2 and momenta of p ∼ mv. At the p µ ∼ mv scale there are corresponding soft gluon and massless quark fields, A µ q and ϕ q , and at the p µ ∼ mv 2 scale there are ultrasoft gluon and massless quark fields, A µ us and ϕ us . We use dimensional regularization with d = 4 − 2 to regulate ultraviolet divergences.
At the matching scale m we integrate out all of the hard modes with momenta of order m as well as other offshell modes which include potential gluons with energy p 0 ∼ mv 2 and 3-momentum p ∼ mv. As a result of integrating out offshell modes with momenta ∼ mv we generate a set of potentials between the heavy quarks, which are given by spatially non-local four fermion operators. In this one-stage approach there is no further matching at the scale mv since the mv 2 and mv scales are treated simultaneously.
We also choose to integrate out soft heavy quarks with p 0 ∼ mv, which avoids simultaneously having both soft and potential heavy quarks in the EFT. 1 Thus when a soft gluon, with p µ ∼ mv scatters off a potential heavy quark with p µ ∼ (mv 2 , m v) (relative to its rest mass), the resulting quark is thrown off-shell with p 0 ∼ mv. Therefore, Compton scattering graphs between potential heavy quarks and soft gluons are matched onto non-local (at scales shorter than 1/(mv)) operators. These operators essentially act like additional potential interactions between the potential heavy quark and soft gluon and soft light quark modes. We will see that a crucial distinction between NRQCD and the SCET with forward scattering interactions, is that this potential-soft interaction gets corrected beyond tree level in NRQCD, which is not the case for the analogous Glauber exchange operators in SCET [29]. We will in fact see that these higher order corrections to this gluon-quark potential are what leads to the new color structure starting at three loops (i.e. the structure that is not "maximally non-Abelian" ∝ C i C k A , where C i is the tree-level Casmir). To distinguish the scales mv and mv 2 , in the EFT we split the external 3-momenta (P full ) of the heavy quark fields into a large label piece (p ∼ mv) and a small residual momentum (k ∼ mv 2 ), so that P full = p + k. To maintain homogeneous power counting, we use the label formalism developed in [3], where momentum or energy of ∼ mv are represented by field labels in momentum space. At tree level the relation between the heavy quark fields in QCD (ψ) and 1 This is not strictly necessary. If we instead keep soft heavy quarks ψs together with the potential heavy quarks ψp, then the structure of the EFT Lagrangian will be somewhat different. There will be interactions between the soft and potential heavy quarks and soft gluons. The purely soft heavy quark Lagrangian will be identical to HQET, except that we must include 0-bin subtractions to avoid double counting with the potential region. The one-stage matching formalism obtained from this construction would be closer to the formalism of pNRQCD. vNRQCD (ψ p ) is (2.1) An analogous relation holds for the heavy antiquark fields χ(x) and χ p (x). The full theory derivative operator is then decomposed as where the label operator P [32] only operates on the momentum label space, In this way all derivative operators ∂ µ acting on fields ψ q (x) will scale as v 2 .
The effective vNRQCD Lagrangian can be split into soft, ultrasoft, and potential components, 3) The ultrasoft Lagrangian involves only ultrasoft fields where G µν u is the ultrasoft field strength and the ellipses are terms that are higher order in v. The covariant derivative here is D µ = ∂ µ + iµ u ι /2 g u (µ u )A µ u , which only contains the ultrasoft gluon field and scales as D µ ∼ v 2 , depends on a coupling at the ultrasoft scale µ u = mν 2 , where ν is the subtraction velocity scale. Here ι = e γ E /(4π) appears because the strong coupling is being defined in the MS scheme. The field strength scales as G µν u ∼ v 4 . For convenience we will often suppress the renormalization Z factors that relate bare and renormalized quantities.
The potential Lagrangian involves both potential heavy quarks and ultrasoft fields and has terms [3,8,12] where L pu contains ultrasoft gluon couplings to potential operators and ultrasoft kinetic corrections to the potentials [13], and the ellipses denote terms of higher order in the v expansion. The potential heavy quark fields scale as ψ p ∼ v 3/2 and χ p ∼ v 3/2 , and B i u = ijk G jk u ∼ v 4 . In terms of the velocity subtraction scale the soft scale is µ s = mν, and spin and color indices in V and the fermion fields have been suppressed. Matching perturbatively at m and integrating out the potential gluons generates the terms where k = p − p, T A andT A are color generators in the 3 and3 representations respectively, Here S 1 and S 2 are the quark and anti-quark spin operators respectively. The relation between the bare and renormalized coefficients is given by V bare These Wilson coefficients were calculated to one-loop in Ref. [13] and the analogous potentials for unequal quark and anti-quark masses were computed at one-loop in Ref. [37].
The soft Lagrangian consists of pure soft field terms, plus interaction terms involving both soft and potential fields, We will use the notation iD s for the soft covariant derivative in position space. In Eq. (2.8) the soft covariant derivative is written in momentum space iD q,q = δ q,q i∂ − g s (µ s )µ s ι /2 δ q,q +p A p , and the gauge fixing has been performed in a general covariant gauge with parameter ξ s , and hence includes soft ghosts c q . The soft fields A µ q , ϕ q , and c q have incoming momentum q, whereas ϕ q andc q have outgoing momentum q. Notice that the gauge fixing of the ultrasoft and soft gluon fields can be performed independently. The ultrasoft fields have momentum components which are parametrically smaller then their soft counter-parts and therefore these fields do not directly couple to each other. The heavy (potential) fermions do not transform under soft gauge transformation, as such a transformation has energies of order mv and would throw them off-shell.
In the formalism of Ref. [3] the soft-potential interaction terms are contained in The notation is such that σ denotes a corrections which is down by v σ . The leading contribution to this soft Lagrangian arises from the Compton graphs shown in Fig.(1), plus the analogous diagrams with soft quarks and ghosts. Taking q µ ∼ mv and p µ ∼ (mv 2 , mv), and keeping only the leading order result in v (which is σ = 0), the result in Feynman gauge is [3] U (0) Here the ghost field, c s , is the same one used in Eq. (2.8). The normalization for the generators in the fundamental representation have been taken to be Tr(T a T b ) = 1 2 δ ab . The terms that are higher order in the velocity with σ = 1, 2 have also been calculated in Refs. [12,38]. The softpotential interactions are responsible for part of the running of the Wilson coefficients potential, V (p, p ) in Eq. (2.5), and they are solely responsible for the running of the leading order Coulomb potential operator until 3-loops [3,16,17].
Note that in the Abelian case the sum of the two Compton scattering graphs cancel, and hence W The delta function forces the gluon to have vanishing energy, i.e. it becomes a potential gluon.
Thus in order not to double count we should not include this twice. The rigorous procedure which enforces this are the vNRQCD "zero-bin" subtractions formulated in [34], which avoids double counting between soft and potential contributions. This will be discussed further below.
We emphasize that the individual results in Eq. (2.10) for σ = 0 and their analogs for σ > 0 are gauge dependent. As previously mentioned, there is nothing wrong in working this way, as the S-matrix elements will always be gauge invariant. However, by not keeping manifest soft gauge invariance we lose a useful organizing principle and there is a corresponding loss of simplicity. The main goal of this work is to setup a formalism which will make the coefficients in the Lagrangian L int s manifestly soft gauge invariant.

Reformulation of the Soft NRQCD Lagrangian
The gauge non-invariance of the Lagrangian L int s arises because even though we have used the onshell conditions k 2 = 0 for the external soft gluons and analogs for the heavy fermions, we have not used our complete freedom to exploit the on-shell transversality condition ∂ · A = 0. In this section we will construct soft gauge invariant operators that match the full theory amplitudes on-shell when exploiting the full freedom of the equations of motion. This same organization was used in Ref. [29] in constructing gauge invariant Glauber exchange operators in SCET.
First we introduce the soft gluon and soft quark gauge invariant building blocks, which in where S v is a soft Wilson line in the v µ direction, which parallel transports the derivative to time-like −∞.
Here v = (1, 0) in the heavy QQ rest frame (we use a roman v to make it clear that this is not the relative velocity scaling parameter v). In particular the fundamental soft Wilson line is and we note that . The equations of motion v · D s S v = 0 obeyed by the soft Wilson line imply that v · B = 0. The iε prescription in the denominators is determined by the limits of integration, and kills the contribution from the endpoint at infinity. Note that we used the same Wilson line on both sides of the covariant derivative in the definition of B µ . This ensures that B µ (x) is a pure octet. In particular using standard relations between fundamental and octet Wilson lines and Eq. (3.1) gives This shows that B a µ is a soft field strength attached to an adjoint Wilson line S v = S v (x, −∞). The building blocks have scaling B µ ∼ v and Ξ(x) ∼ v 3/2 . They both start off linear in the soft gluon or quark fields, respectively, for which they are soft gauge invariant extensions. In momentum (label) space using Eq. (3.2) leads to where here k 1 + k 2 = k, and all denominators have k 0 i = k 0 i + iε with our convention for the soft Wilson lines. Note that the operators B µ and χ s transform under a global gauge transformation (equivalent to its transformation at −∞), but this transformation will be canceled by other fields in the interaction Lagrangian we will construct below.
The convention for the direction of the Wilson lines in Eq. (3.1), or equivalently the sign of the iε in Eq. (3.4), should be irrelevant since the pole only matters in the potential region as described in Eq. (2.11), and hence is not a part of the soft sector of the theory 3 . Therefore, the direction of the Wilson line must be irrelevant if we perform the proper potential zero-bin subtraction [29].
To demonstrate this we consider an alternate definition for the soft gluon building block operator This definition still leads to an octet B µ field, but it now gives a field strength attached to an . Unless otherwise specified we will by default work with the definition in Eq. (3.1).
At lowest order in the v power counting we now consider building operators for the interaction Lagrangian L int s with two potential heavy-quark fields and two soft gluon building blocks B µ (which is the minimum number allowed that is consistent with soft momentum conservation), This product of fields is O(v 5 ), and according to the power counting theorem of Ref. [3] it must be multiplied by a coefficient C(q, q , p, p ) ∼ v −1 to give a leading power contribution. (Although this will make it leading order in v, the interaction also starts at O(α s ), and hence is perturbative if µ s Λ QCD .) Imposing the fact that the leading power coefficient must also be m independent, and can not be more singular than linear or quadratic in a momentum, due to the locality of QCD, 5 we can eliminate the possibility of contracting the indices µ and ν 3 Similar reasoning in the EFT of black hole binaries [39] explains apparent IR divergences in the calculation of gravitational potentials at O(v 8 ) [40]. 4 One may wonder why there is no such operator in HQET given that in the one particle sector there should be no distinction between NRQCD and HQET. However, we must keep in mind that these operators will be inserted into time ordered products and thus the relevant dispersion relation for the quarks will always be that of NRQCD.
In HQET the Compton scattering is reproduced by a time ordered product whereas in NRQCD we use these higher dimensional, leading order in v, operators with the soft heavy fermion integrated out. 5 Here we are restricting ourselves to two to two scattering of point particles.
with momenta, and must therefore have a g µν at leading power. There are three independent ways to construct a color singlet out of 3,3, and two 8s, given by the following tensors if abc T c αβ , d abc T c αβ , δ ab δ αβ . This gives are functions of (p, p , q) with q = p − p − q, and g s = g s (µ s ) as before. Note that only spatial B i s appear since v · B = 0. It is convenient to include the factor of 1/2 since the Feynman rule will include terms where the contraction with

Leading Power Soft-Potential Interactions
Let us now match the leading order operator to obtain U (0) ij . Its value is fixed so as to reproduce the Compton scattering tree level graphs in Fig.1. Consider first matching only the spatial polarizations for both gluons, for which the only a contribution stems from the diagram with a 3-gluon vertex. We find which agrees with the result for U (0) ij from Eq. (2.10). Naively it would appear that with this fixed value of the coefficient, we will not match the time-like polarizations. However, gauge invariance implies that on-shell the results for the time-like polarizations must also match. To check this explicitly, we compute the tree level full QCD amplitude involving time-like polarizations from the graphs in Fig. 1, and write it back in the form of soft gluon fields µA obtain as the structure appearing in ig 2 ψ † p M full ψ p when expanding the amplitude to leading order in v. To derive this result we have made use of momentum conservation and the equations of motion To obtain the analog of M full in the effective theory we simply replace each of B q and B q in Eq. (3.7) by the terms with one gluon field from Eq. (3.4) to obtain the terms involving at least one time-like polarization, again substituting µA Using the onshell conditions q 2 = q 2 = 0 and momentum conservation implies that q · q = 1 2 (p − p) 2 − (q 0 ) 2 , and hence that Eq. (3.9) matches exactly with Eq. (3.8) onshell. From tree level matching the remaining leading power coefficients are given by Note that since we are matching on-shell we no longer have any analog of the operator in Eq. (2.9) that involves an interaction of soft ghosts with potential heavy quarks. Eq. (3.6) alone reproduces the correct gauge invariant anomalous dimension results as discussed further in the next section.
While at tree level W ij vanish, these coefficients will be generated by the matching at one loop. In the full theory we have the three one loop diagrams shown in Fig. 2 which can contribute to the d abc and δ ab color structures, while there are no one-loop vNRQCD diagrams with these color structures to subtract. These diagrams generate non-zero results W In the language of the threshold expansion [10], these one loop matching results come from the potential region of the full theory loop integrals. In vNRQCD these contributions generate potentials for the soft gluons as in Eq. (3.6). There will also be one loop corrections to U  one-loop vertices yields three loop corrections to the static potential which contribute to the color structure [41][42][43][44].
We may also investigate operators that have more building block fields. For example, we may consider operators with three B's, ψ † p B i q B j q B k q ψ p with all possible color structures. It is straightforward to show that the matching onto these operators vanishes at tree level. To match we can consider full theory graphs which are leading power and have three spatial polarizations for the gluons. Since the attachment to the heavy quark has a time-like polarization in the full theory at leading power, the only such diagram is the one shown in Fig. 3(a). The result for this graph is exactly reproduced on-shell by the vNRQCD graph in Fig. 3(b) with all spatial polarizations for the gluons. Thus the tree level matching onto the three B operator leads to a vanishing Wilson coefficient. Examining operators with more spatial gluons does not change this tree level correspondence between QCD and vNRQCD graphs when we have a single gluon attached to the heavy quark, generalizing the graphs in Fig. 3. This remains true even when there is more than one gluon attached to the heavy quark. Thus all operators of the form ψ † p B i 1 q 1 · · · B in qn ψ p with n ≥ 3 have vanishing Wilson coefficients at tree level. When matching at higher loop orders the coefficients of these operators may no longer vanish.

Power Corrections to Soft-Potential Interactions
The above analysis can be extended to use the gauge invariant soft building blocks to include terms at subleading orders in the v expansion in Eq.(3.6). The Wilson lines in the B µ q and Ξ q contain all the coupling to the time-like gauge boson polarizations, including terms like U (σ) 00 and U (σ) 0i in Eq. (2.10). Therefore we may extract the required subleading power matching results from [12,38] by keeping only the results with spatial indices, which to 1/m 2 read This provides a much simpler set of matching coefficients than the full list of terms including direct A 0 couplings. Here c F (µ S ), c S (µ S ), and c D (µ S ) are Wilson coefficients in HQET Lagrangian whose renormalized low scale values can be found in [7] and whose leading logarithmic RGE equations can be found in [45].

Running of the Leading Power Potential
To test some aspects of our gauge invariant formulation of L int s we consider in this section the soft one-loop diagrams that are responsible for the running of the Coulomb potential. We also discuss the independence of results to the direction chosen for the Wilson lines in B µ .
The one loop diagram shown in Fig. 4 arises from time ordered product of two U (0) ij operators and is responsible for leading log running of the Coulomb potential coefficient V c . Since the softfermion loop calculation is identical to that considered in earlier papers (see for example Ref. [3]), we will only consider the contributions from soft gluons. This calculation can be simplified by considering the contraction of the full B µ field, rather than the soft A s fields it contains (where A s is still the quantized field). Since B 0 = 0, we only need the time ordered product of two spatial B i fields. To leading order in g this is given by the gauge invariant result Using dimensional regularization with d = 4 − 2 for the graph in Fig. 4 we find where we have included a symmetry factor of 1/2, and theT A is for the antiquark which is in the 3 representation. Here k = p − p is the momentum transfer and the loop integrand is To go from the first to second line of Eq. (4.3) we have written all terms in the numerator in terms of either propagator denominators, or q 2 , (q 0 ) 2 , or k 2 , and dropped power law divergent scaleless integrals. The first term in the integrand in Eq. (4.3) gives a standard integral, For the second term in the integrand in Eq. (4.3) we have a pinch singularity from the q 0 ± iε denominators, as can be verified by working out terms in the residue theorem for the q 0 contour integral. However this pinch occurs from the potential region q 0 ≈ 0 where the approximations used to derive the soft Lagrangian are not valid. In the vNRQCD theory the contribution is removed from the soft integrand by a zero-bin subtraction [34] that is intrinsic to the definition of the soft gluon propagators. Defining T s2 (q 2 0 , q, k) = q 4 2 +q 2 k 2 + k 4 4 / [q 2 + iε][(q + k) 2 + iε] the second term in the integrand with the potential zero-bin subtraction included gives the integral Combining the results for I 1 and I 2 then gives the final result which is precisely the expected result for both the 1/ uv and constant terms from the soft region at this order.
Combining Eq. (4.6) with the contribution from n f soft fermions (which using Ξ q is still the same calculation as in Ref. [3]), we get the full beta function of QCD, which runs the Coulomb potential at one-loop, Note that we reproduce the correct running with the factor of 11C A /3 despite the absence of ghosts in Eq. (3.6). 6 This occurs because the B two point function itself is gauge invariant.
Beyond one loop, the ghosts in Eq. (2.8) will contribute to internal loops, and hence will appear at intermediate steps when determining the O(g 2 s ) correction to Eq. (4.1). But once all contributions are added up the result for the B two point function will again be gauge invariant. It should also be noted that all of the 1/ poles in Eqs. (4.4) and (4.5) are UV in origin as indicated by writing uv . Naively one may have thought there could be IR poles in the integral I 2 , but these cancel out between the terms in the numerator. In general if there are 1/ IR poles in a soft loop graph they are converted to 1/ uv divergences by ultrasoft 0-bin subtractions [34] (which field theoretically implements the pull-up mechanism [16]). However there are no ultrasoft zero-bin subtractions in these diagrams. This is natural because the ultrasoft gluons decouple from the potentials at leading order. The easiest way to see this is to note that we can make a field redefinition to decouple ultrasoft gluons at leading power from the D 0 term in Eq. (2.4), which is analogous to the BPS field redefinition [33] in SCET,

Simplifying the Running of the O(v 2 ) potential
Next we look at the running of O(v 2 ) potential in NRQCD at one loop. In vNRQCD the running of these subleading operators was performed in Refs. [12,17,38] in the non-gauge invariant basis given in (2.9), which involved many more operators due to the need to keep track of the timelike polarizations and their associated Feynman rules. Here we will repeat the calculation of time-ordered products involving U The contribution from the time ordered product with two U (1) ij vertices in Fig. 4 is Since W (0) ij = 0 we only have a contribution from two W (1) insertions which give Adding up all the contributions we find the known result [38] which contributes a soft-loop contribution to the running of the potential coefficients V here. Ultrasoft modes cannot contribute in the one body sector [38], instead they renormalize the time ordered product of two or more soft vertices with quarks and antiquarks as in Fig. 5(a) through graphs like those in Fig. 5(b) and Fig. 5(c). These products appear local as far as ultrasoft gluons are concerned and it is only these products which affect observables. For a complete discussion of this subject we refer the reader to [17,26]. Here we are only interested in showing that results obtained in [17] follow in our manifestly soft gauge invariant formulation. To setup the ultrasoft renormalization from graphs involving the time ordered product of two σ = 0 vertices from L int s in Eq. (3.6) we define the operators where the Γ's are functions of the momenta, for example Γ B,abcd (p, p , q). Note that we do not have an operator in Eq. (4.13) with soft ghost fields, unlike Ref. [17]. These structures are determined by the form of the time-ordered products in Fig. 5(a), so for example the structure in the first operator is: 14)   Fig.5(b,c), are . Since ultrasoft modes do not couple to soft modes, an ultrasoft loop will generate the same Wilson coefficient for the operators with B µ fields as it would if we used A µ fields instead of B µ fields.

Higher Order Loop Corrections to the QCD Potential
In the full theory the static potential V (R) is often defined by the time like Wilson loop In vNRQCD the leading power Coulomb-like potential is obtained from a combination of the V (1,T ) c potential matching coefficients in Eq. (2.6) and time ordered products involving soft fields.
The soft time ordered products also induce the running of the strong coupling in the potential.
In the Wilson loop calculation of V (R) using Eq. (5.1), at three loops one encounters an infrared divergence known as the ADM singularity [46]. However, this divergence is not part of the potential [23], as it arises from the ultrasoft region of integration and thus must be subtracted.
In the vNRQCD formalism this divergence is removed by an ultrasoft zero-bin subtraction to the soft loop computation [34]. The two loop static potential has been calculated in Refs. [35,36], and the three loop soft potential has been calculated in Refs. [41][42][43][44].
An interesting aspect of the two loop static potential, pointed out in [36], is that the iε's in the propagators of soft heavy quarks must be accounted for to obtain the correct result. In our one loop discussion we remarked that the sign of these iε's only contributes a term where the energy of a soft gluon goes to zero, i.e a potential exchange that is removed by the zero-bin subtractions. At two loops, the calculation in [36] indicates that this is no longer the case, and we will explore how this relates to our formalism with B operators in this section.
In particular, it should be noted that we do not have complete freedom to pick the iε's that appear in Feynman rules associated with B µ , because the allowed structure of Wilson lines is constrained by requiring that B µ = B µa T a is purely an octet field, as in Eq. (3.3), which in turn requires the use of the identity S † v S v = 1 and, for example, would not be the case if we tried to Figure 7. Two loop EFT diagram for running of the leading order potential that arises when using the gauge potential as the canonical variable and this absent when quantizing using B.
This is consistent with the fact that not all iε's from the soft Wilson lines in B µ lead to iterated potentials. As we will see below, for some of these terms there is no zero-bin subtraction and therefore the soft loop integrals have non-trivial contributions where these iε's are important.
Working with the B µ fields in the EFT introduces some apparent technical complications that we will now show are straightforward to handle. In particular, when using Eq. (5.1) any integral that is ill defined due a pinch in an energy contour integral, need not be calculated since non-abelian exponentiation [47,48] implies that they are iterations of lower order potentials terms from the matrix element in Eq. (5.1) that are proportional to T k≥2 , where T is defined in Eq. (5.1) . In contrast, in vNRQCD the pinch singularities in soft loop calculations are removed by potential zero-bin subtractions [34]. Also, the contributions which need not be calculated in the Wilson loop calculation, due to non-abelian exponentiation, do not even appear from the vNRQCD Feynman diagrams generating purely soft loops. However, in our formalism using B's the non-trivial soft loop contributions can themselves come with pinch singularities, unlike results obtained from static potential calculations starting with Eq. (5.1). Thus in the EFT we must deal with pinched integrals to calculate the potential, and the corresponding zero-bin subtractions.
As an example of how to properly perform the zero-bin subtractions, consider the two loop integral which arises from the diagram in Fig. 7. We will neglect the terms in the integrand that do not have the maximal number of time-like k 0 propagator factors since this allows us to focus on the most complicated term, which has overlapping subtractions. After expanding out the B fields in terms of soft gluon fields to higher order using Eq. (3.4), and performing the contractions, the relevant integral is is the net momentum exchange between heavy quarks and We have also defined , (5.4) which is antisymmetric under k 0 i → −k 0 i . Once again we find that this same result is obtained whether we start with Eq. Let us simplify this term in the integrand by putting it into a form which makes all the pinch singularities manifest. We will use the notation in [36], introducing and defining Making the change of variables such that k 3 → k − k 1 , and using Eq. (5.4) we can simplify the first term in the integrand in Eq. (5.2) as follows: The last equality follows from the symmetry of the integral under v → −v and under {k 1 → k − k 2 , k 2 → k − k 1 }. The last term in Eq. (5.8) has a double pinch singularity, and the second to last term has a single pinch singularity, whereas the first two terms have no pinch singularities.
All four terms will have zero-bin subtractions, but the zero-bins of the first two terms will vanish upon contour integration. The iε factors in the S 1122 and S 1122 are nevertheless important for evaluating these integrals. In contrast, the last two terms in have Eq. (5.8) have subtractions that must be included, which we can check by examining their power counting in the potential limit.
Consider the scaling of the term S 1122 in Eq. (5.8) which has a double pinch singularity, and call its full integralJ prior to including any subtractions. HereJ ∼ 1/v 2 is the scaling of this integral when both loop momenta are soft. The relevant zero-bin limits to consider and resulting scaling for the corresponding zero-bin loop integrals J From this we conclude that all three of these zero-bin subtractions must be included, since For example, to construct the first integral we take k 0 1 ∼ v 2 , k 0 2 ∼ v and hence expand the integrand ofJ with k 0 1 k 0 2 , k 1 , k 2 . This gives the zero-bin integral . This discussion follows the standard zero-bin approach for multiloop integrals with overlapping subtractions. These terms give the same integral as the third k 0 1 ∼ k 0 2 ∼ v 2 subtraction listed in Eq. (5.9). This integral is .
Thus the net effect of including J subtraction. For the total integral combining these terms then gives Altogether this leads to a combined result that has well defined integrals in both k 0 1 and k 0 2 . Although here we are discussing graphs that contribute to the C 2 A (T e ⊗T e ) color structure, the nature of the zero-bin subtractions for the S 1122 contribution is actually identical to the two loop double box graph in the HQET based loop computation of the potential.
The S 1122 term in Eq. (5.8) is even simpler to handle since it only has a single non-zero zerobin subtraction for k 0 2 ∼ v 2 . The analysis for S 1122 is identical to the two loop graph that is a cross box with an extra gluon rung added. Thus all pinch singularities in Eq. (5.8) are converted to well defined integrals once the zero-bin subtraction terms are included.
Note that the first two terms inside the square bracket in Eq. (5.2) are equal by the symmetry of the original integral under k 1 ↔ k 2 , so our discussion so far also suffices to cover the second term in this expression.
where k 5 was defined in Eq. (5.6). In the final expression here the first three terms do not have pinch singularities, and can be directly computed using standard two-loop techniques. In general the iε factors in the eikonal propagators are needed to obtain correct results for these terms. In contrast, the last three terms in Eq. (5.13) have pinch singularities. This is immediately apparent for the fourth and fifth terms which involve a S 11 . For the sixth term the pinch singularity is revealed after performing one of the energy integrals, for example integrating in k 0 2 yields an answer with a pinch singularity in k 0 1 . Since the nature of the zero-bin subtractions for the last three terms in Eq. (5.13) is quite similar, we will choose as an example to go through the analysis for S 1125 with only the k 4 term kept in the numerator and again pulled outside the integrand.
We call the corresponding full integralL prior to including any subtractions.
Once again we first enumerate the necessary zero-bin subtractions by considering the scaling of theL integrand and measure in potential limits: 14) Note that here because of the presence of a S 5 term that we also considered a k 0 5 ∼ v 2 potential limit. The non-zero zero-bin integrals are .
Again there is an additional subtraction made on L (k 1 ) 0 to ensure its k µ 2 momentum remains soft (or equivalently, that this integral does not double count the result from L ). This term is given by .

(5.16)
Note that in this case L . The final purely soft result is then obtained by To evaluate this result we first note that after combining integrands the combination −L no longer has a pinch singularity in k 0 2 , and has a k 0 1 integral that yields zero. This leaves L =L − L (k 1 k 2 ) 0 . Putting these two integrands over a common denominator gives terms in the numerator involving energies as {(k 0 1 ) 2 , (k 0 2 ) 2 , k 0 1 k 0 2 } or with higher powers. These all yield well defined integrals. Thus the zero-bin again removes the pinch singularities, yielding a well defined result for L. The analysis for the S 1125 and S 1125 terms in Eq. (5.13) is similar, with the subtractions again yielding well defined final results.
In the above analysis, only the pinched poles are removed by zero-bin subtractions, while the other ±iε terms in the eikonal propagators are important for obtaining the final results for integrals. These poles can lead to additional π 2 terms, which are related to the contribution discussed by [36]. We thus find that the precise path for the Wilson lines used to define B is again not relevant at two-loops (and presumably to all orders). However, this does not imply that the signs of the iε obtained from expanding the soft Wilson lines are not needed to obtain the correct results. Indeed, it is important that B is defined as an octet field. With this constraint satisfied, changing the path of the Wilson lines flips the signs of all eikonal iε's in a correlated manner, and these correlations matter.
Finally we note that working in the effective theory can illuminate higher order corrections to the potential. First of all, there are fewer diagrams since there is not a term with a single soft gluon coupling to the potential quarks. The color structure is also clarified since in the EFT only terms which actually contribute to the potential arise in matching. For instance, at two loops the soft contribution of all the integrals proportional to C 3 F as well as those proportional to C 2 F C A must sum to zero if we use the Wilson loop definition of the potential (or using the HQET implementation of the soft sector), whereas in the treatment advocated here no such color structures in soft integrals arise.

Conclusions and Outlook
We formulated the interaction terms between potential heavy quarks and soft fields in vNRQCD in a way which is manifestly soft and ultrasoft gauge invariant. This improves upon the construction of Ref. [3] as it greatly reduces the size of the operator basis for this Lagrangian, and leads to simpler calculations for soft loop diagrams at both leading and subleading power. Most importantly perhaps this formulation in terms of the soft gauge invariant building blocks B and Ξ is likely to lead to simplifications for higher order matching and anomalous dimension calculations in the soft sector.
Our formalism also helps to elucidate the role of the iε prescription in soft heavy quark propagators. These appear from soft Wilson lines in the formalism used here, and give additional pinch singularities that are then removed by zero-bin subtractions from the potential region of momentum space. Interestingly, the zero-bin subtractions in this analysis do not have a 1-to-1 correspondence with EFT diagrams with potential loop momenta, unlike the typical situation that is encountered for zero-bin subtractions. This is easily seen since the zero-bin subtractions Given the analogy between the use of our B for vNRQCD, and the use of the soft and collinear building block fields for Glauber operators in SCET [29], it seems quite likely that a similar conclusion will apply there as well.
An interesting possibility raised by the formalism developed here is to consider quantizing the soft B µ and Ξ fields themselves rather than the soft gauge field and soft quark field. An analogous method of quantization fields that involve collinear Wilson lines was considered for SCET in Ref. [49], but has not yet found wide use there. 7 If we consider this approach in 7 For the same reasons discussed here, it may however prove to be useful for calculations involving Glauber vNRQCD, then it leads to even simpler Feynman diagrams since Eq. (3.6) now yields four point vertex Feynman rules, but no higher point vertices. With soft fields quantized in this manner, the two point function for B is given by Eq. (4.1) dropping terms at and beyond O(g 2 ) since the B fields are now fundamental. Furthermore, the three and four point interactions for Bs are identical to those for the triple and quadruple gluon vertices in QCD [49]. We find that considering this approach in vNRQCD is like working in a v · A s = 0 gauge for the soft gluon fields, utilizing zero-bin subtractions to handle the extra singularities in the gauge propagator that otherwise make this gauge complicated (see eg. [50]), while systematically avoiding power counting violating corrections that are present for HQET in this gauge [51]. 8 It would therefore be interesting to further explore the direct use of a quantized B and Ξ for carrying out soft NRQCD calculations.
Another interesting idea for further exploration is using the field redefinition in Eq. (4.8) to obtain a simpler set of ultrasoft couplings to potential heavy quarks. Just like in SCET this will lead to the appearance of ultrasoft gauge invariant gluon building block fields through the identity Y † v iD µ us Y v = i∂ µ us − gB µ us , where v · B us = 0. Here B us is defined in the same manner as Eq. (3.1) but with the S v Wilson lines replaced by the Y v Wilson lines.