FeynOnium: Using FeynCalc for automatic calculations in Nonrelativistic Effective Field Theories

We present new results on FeynOnium, an ongoing project to develop a general purpose software toolkit for semi-automatic symbolic calculations in nonrelativistic Effective Field Theories (EFTs). Building upon FeynCalc, an existing Mathematica package for symbolic evaluation of Feynman diagrams, we have created a powerful framework for automatizing calculations in nonrelativistic EFTs (NREFTs) at tree- and 1-loop level. This is achieved by exploiting the novel features of FeynCalc that support manipulations of Cartesian tensors, Pauli matrices and nonstandard loop integrals. Additional operations that are common in nonrelativistic EFT calculations are implemented in a dedicated add-on called FeynOnium. While our current focus is on EFTs for strong interactions of heavy quarks, extensions to other systems that admit a nonrelativistic EFT description are planned for the future. All our codes are open-source and publicly available. Furthermore, we provide several example calculations that demonstrate how FeynOnium can be employed to reproduce known results from the literature.


Introduction
In the last decades we witnessed how Effective Field Theory (EFT) methods [1,2] were successfully applied to describe various phenomena governed by electromagnetic, weak, strong and gravitational interactions. A modern pedagogical introduction to the main ideas and techniques of EFTs can be found e.g. in [3][4][5]. Taking advantage of the hierarchy of widely separated dynamical scales found in many physical systems, we can construct suitable EFTs that precisely capture the behavior of the given system at energies below a certain scale. The resulting EFT, which is based on the underlying symmetries, relevant degrees of freedom and power-counting rules, allows us to describe low energy phenomena in a simple but yet rigorous and systematic way.
From the technical point of view, calculations in EFTs are organized as expansions in small dimensionless parameters (e.g. ratios of energy scales). The power-counting rules of the theory precisely tell us where the expansion should be truncated to achieve the precision we are aiming at. Even though the leading order predictions can be often obtained in a short pen and paper calculation, the usage of EFT methods does not imply that everything becomes trivial. On the contrary, the determination of higher order corrections routinely necessitates the usage of elaborate codes for automatic calculations. The time needed to develop such codes and subsequently run them on a sufficiently powerful computer often becomes a bottleneck in the task of obtaining higher order EFT predictions matching experimental accuracies.
Nonrelativistic Effective Field Theories (NREFTs) constitute a subbranch of EFTs for describing systems, where the relevant velocity scales are typically much smaller than the speed of light. Examples for such systems are nonrelativistic bound states (e.g. positronium, muonium [6], heavy quarkonia [7,8]) or systems made of nonrelativistic atoms [9] and molecules [10]. Even though NREFT methods are most commonly employed in nuclear and atomic physics, nowadays they are becoming increasingly popular for studying possible beyond the Standard Model scenarios such as nonrelativistic dark matter [11][12][13][14][15][16][17][18] or heavy neutrinos [19].
One can roughly distinguish between two ways to approach perturbative calculations in NREFTs. The first method attempts to "hide" the nonrelativistic nature of the theory by rewriting (whenever possible) operators and amplitudes in terms of Lorentz covariant quantities. In return, one hopes to benefit from existing codes for automatic calculations and to avoid dealing with nonrelativistic expressions as much as possible. The other approach is to embrace the loss of manifest Lorentz covariance and to perform the calculations directly, working with nonrelativistic integrals, Pauli matrices and Cartesian tensors.
In our view, the best strategy consists of finding the right balance between the two approaches without sacrificing any physical insight or computational convenience. In particular, we believe that depending on the size of the calculation and the questions one is trying to answer, it may be sometimes more advantageous to work with noncovariant quantities directly, rather than trying to eliminate them altogether.
However, since nowadays most calculations are carried out using computer codes, the choice between the two above mentioned approaches has also a technical dimension.
However, this is not the case once one becomes interested in performing the given calculation is a nonrelativistic fashion. While we do not see any intrinsic difficulties that make the automation of nonrelativistic calculations more challenging than the relativistic ones, one can hardly find any public codes applicable to this scenario. In principle, nothing prevents a programming-savvy user to implement the necessary operations in FORM [27,28], Mathematica, Maple or any other symbolic manipulation system. This is also what most practitioners usually do, when they face the necessity of carrying out a largescale nonrelativistic calculation without being in the possession of suitable in-house codes. Unfortunately, the NREFT community suffers from a visible lack of interest in making such codes publicly available, which effectively means that a lot of people have no other choice than to write their codes from scratch.
In our view, this situation is very unfortunate and deserves to be improved. Our contribution to the solution of this problem and the novelty of this work is, therefore, to provide software packages that are publicly available, well documented, easy to use and most importantly suitable for automatizing nonrelativistic calculations. What is more, we will also explicitly show how these tools can be used to reproduce important results from the literature.
The FeynOnium project started in 2016 [29] and its main focus is still directed towards tree-and 1-loop level calculations in Nonrelativistic QCD (NRQCD) [6,7] and potential Nonrelativistic QCD (pNRQCD) [30,31] as well as the electromagnetic counterparts Nonrelativistic QED (NRQED) and potential Nonrelativistic QED (pNRQED). However, as it will become clear in the course of the paper, most of the provided routines are in no way limited to a particular theory and can be employed in very generic nonrelativistic calculations. Our key deliverables are a new version of the FeynCalc package [32], capable of dealing with nonrelativistic quantities out of the box and a special add-on (also called FeynOnium) for NREFTs. This paper is organized in the following way. In section 2 we introduce NRQCD and pNRQCD, which will appear in many of our example calculations. To set the stage for our tools, we provide a brief overview of the existing codes for EFT calculations in section 3, while our technical implementation is described in section 4. The installation and usage of the packages are explained in section 5. In section 6 we demonstrate how the presented codes can be employed to reproduce some well-known (NR)EFT results from the literature. Our conclusions and possible future directions of this work are summarized in section 7. We provide useful formulas for algebraic manipulations of Pauli matrices in appendix A, and we list the Lorentz and Cartesian tensors in FeynCalc's internal (FCI-notation) and external (FCE-notation) notations in appendix B. Finally, we present derivations of NRQCD and pNRQCD Feynman rules in appendix C.

Nonrelativistic QCD and potential nonrelativistic QCD
NRQCD is an effective field theory of QCD that is appropriate for describing bound states of a heavy quark and a heavy antiquark, like heavy quarkonia. The heavy quark and the heavy antiquark with mass m have a typical velocity v inside the bound state, which is the small expansion parameter of this EFT. The degrees of freedom of NRQCD are the two-component Pauli spinor fields ψ and χ that describe the heavy quark and the heavy antiquark, respectively, which interact with the gluon field A through the Lagrangian, given up to order 1/m by [6,7] q i are massless quark fields with flavor i, and c n are the matching coefficients of NRQCD. The heavy-quark mass m that appears in the NRQCD Lagrangian is the pole mass. At order 1/m 2 and beyond, operators of higher dimensions appear, which include heavy quark bilinears, four-quark operators, and gluonic operators. The velocity expansion in NRQCD is an expansion in powers of momentum divided by the heavy quark mass m, where the momentum may scale like mv, the typical size of 3-momenta of the heavy quark and antiquark in the quarkonium, or mv 2 , the typical size of the binding energy. This is similar to the heavy-quark effective theory (HQET) [33][34][35][36][37], where the expansion parameter is Λ QCD /m. Both cases can be regarded as a formal expansion in powers of 1/m, and the two effective field theories have the same Lagrangian in the two-fermion sector, although they have different power counting rules.
The matching coefficients c n are determined by requiring the EFT to reproduce QCD for processes involving nonrelativistic heavy quarks. That is, where M QCD is a QCD amplitude, O n are NRQCD operators, with corresponding matching coefficients c n . The NRQCD matrix elements B|O n |A are computed with the same initial and final states as the QCD amplitude. The NRQCD matrix elements scale in v, and hence, the sum over n is organized in powers of v. In practice, in order to work with a finite number of NRQCD matrix elements, the sum is truncated at a given order in v.
The matching coefficients c n can be determined perturbatively, by computing the QCD amplitude on the left-hand side of eq. (2.2) in perturbative QCD (pQCD), and the NRQCD matrix elements on the right-hand side of eq. (2.2) in perturbative NRQCD. The Feynman rules of perturbative NRQCD are listed in appendix C.1. Then, the c n are determined by requiring that the right-hand side of eq. (2.2) reproduces the perturbative QCD amplitude on the left-hand side to a desired accuracy in an expansion in powers of the momenta of the heavy quarks, antiquarks, and soft gluons that appear in the perturbative amplitude.
Unlike what is usually done in perturbative QCD, perturbative NRQCD calculations are organized in terms of nonrelativistic quantities like the 3-momenta of quarks and gluons. Two-component Pauli spinors and Pauli matrices handle the heavy-quark spin. On the other hand, amplitudes in perturbative QCD are usually given in terms of relativistically covariant quantities, like 4-momenta, gamma matrices and Dirac spinors. Hence, in order to compute c n from eq. (2.2), it is necessary to rewrite the pQCD amplitude in terms of nonrelativistic quantities so that it can be compared with the NRQCD matrix elements. That is, we need to rewrite 4-momenta in terms of 3-momenta, gamma matrices in terms of Pauli matrices, and Dirac spinors in terms of two-component Pauli spinors. This can involve a considerable amount of nonrelativistic algebra that is best done on a computer. NRQCD involves two dynamical scales mv and mv 2 . When mv Λ QCD , the scale mv can be integrated out perturbatively to obtain a new effective field theory called potential NRQCD (pNRQCD). The degrees of freedom of pNRQCD are a singlet and an octet field, low energy gluons and light quarks. Since for the lowest quarkonium resonances the typical size of the relative coordinate is smaller than the inverse on the confinement scale Λ QCD , one can employ the multipole expansion at the Lagrangian level. The heavy quark sector of the pNRQCD Lagrangian in the weakly coupled case (r Λ −1 QCD ) at next-to-leading order in the multipole expansion and at leading order in the 1/m expansion is given by [30,31] L pNRQCD heavy quark where the S and O are the singlet and octet fields, respectively, that depend on time, the relative coordinate r and the center-of-mass coordinate R. They have the following color indices: (2.4) where N c is the number of colors and T F = 1/2. All gluon fields in eq. (2.3), such as the chromoelectric field E i = G ai0 T a and the covariant derivative are evaluated at the center-of-mass coordinate R. The light part of the pNRQCD Lagrangian is the same as in eq. (2.1). The singlet and octet Hamiltonians h s and h o can be split into a kinetic term (for simplicity we just write the leading one) and a potential h s (r, p, S 1 , where V s and V o are the color-singlet and color-octet potentials, respectively. The functions V A (r) and V B (r), as well as the potentials V s (r, p, S 1 , S 2 ) and V o (r, p, S 1 , S 2 ) are the matching coefficients of pNRQCD. The matching between NRQCD and pNRQCD in the weak coupling regime can be carried out by requiring the Green's functions in NRQCD and pNRQCD to be equal order by order in 1/m, α s and r. In the case of perturbative matching, Green's functions in pNRQCD are computed using the pNRQCD Feynman rules, which are listed in appendix C.2. As it was the case for the matching between QCD and NRQCD, also the matching between NRQCD and pNRQCD involves 3-dimensional vectors, which calls for a computer environment that can deal with nonrelativistic algebra.

Existing approaches to automatic EFT calculations
Many of the publicly available packages for EFT calculations aim at exploring the phenomenology of the Standard Model extended with nonrenormalizable operators, the socalled Standard Model Effective Field Theory (SMEFT) [38,39], at tree-or 1-loop level. In practice one usually integrates out heavy fields from an assumed underlying theory of the physics beyond the Standard Model and constructs the corresponding EFT operators. Other common tasks include deriving operator bases from the given set of symmetries, switching between different operator bases, computing the renormalization group evolution of Wilson coefficients or extracting Feynman rules from the effective Lagrangian.
An explicit evaluation of Feynman amplitudes for a given process usually lies beyond the scope of such packages. This part of the calculation can be accomplished e.g. by exporting the Feynman rules for the relevant part of the effective Lagrangian to some common format and then employing suitable codes for perturbative calculations. For example, a model in the UFO [54] format can be imported into popular tools such as Mad-Graph5_aMC@NLO [55], GoSam [56,57] The authors of the above-mentioned EFT codes often stress that their packages are not limited to SMEFT, but can be also employed for more generic theories. While this is certainly true, such theories are nonetheless expected to be manifestly Lorentz covariant, which is problematic for NREFT calculations. Additional limitations equally apply to EFTs that contain nonstandard propagators such as eikonal propagators in HQET, soft-collinear effective theory (SCET) [65][66][67][68][69][70] or chiral perturbation theory (ChPT) [2,71,72].

Some aspects of dark matter studies in an EFT framework can be automatized with
DirectDM [73], which can match the user-provided relativistic high-energy theory onto a low-energy EFT in which dark matter particles interact with nonrelativistic nucleons [74,75]. The matching is nonperturbative and is done at leading order (LO) in the chiral expansion (i.e. one expands in the momentum transfer instead of the strong coupling constant). Furthermore, the determination of the Wilson coefficients is performed in a fully automatic fashion and does not require any explicit manipulations of nonrelativistic quantities. This is very different from the approach we follow in FeynOnium, where the user is required to carry out the matching calculation explicitly but can do so in a much more flexible way.
As far as EFTs of strong interactions are concerned, the number of useful publicly available codes is rather low. For mesonic ChPT, the package Ampcalculator [76] can be employed to automatically calculate selected processes up to 1-loop. PHI 1 , a FeynCalc add-on developed by F. Orellana to generate and manipulate amplitudes in generic ChPT processes is, unfortunately, not compatible to the current version of FeynCalc. Integrals arising from propagator diagrams in HQET (up to 3-loops) can be automatically evaluated with Grinder [77], a special package available for REDUCE and Axiom computer algebra systems. In the case of SCET, numerical calculation of soft functions at NNLO is possible with SoftSERVE [78]. Regarding NRQCD, tree-level amplitudes for heavy quarkonium production and decay processes can be generated with MadOnia [79] or HELAC-ONIA [80,81]. A library of amplitudes for the heavy quarkonium hadroproduction at NLO that were already evaluated with the private FDC [82] code is available via the FDCHQHP package [83].
When applying EFT methods to strong interactions, many practitioners prefer to rely on their in-house codes, which often combine multiple public and private tools in one framework. The first step usually involves the diagram generation, which can be accomplished with FeynArts or QGRAF [84]. After that, the output can be processed with suitable FORM or Mathematica codes, although one might also want to use other computer algebra systems such as Maple, Reduce or Redberry [85]. The codes can be either completely self-written or based on publicly available tools like FeynCalc and Feyn-CalcFormLink [86]. After having carried out all the necessary algebraic simplifications, one would like to evaluate the resulting loop integrals either symbolically or numerically. The loop integral calculus is an interesting topic on its own and we refer to [87] for a pedagogical introduction to the existing methods. Let us merely remark that the simpler (only few mass scales and legs) 1-loop EFT integrals can be very often calculated analytically via a direct application of the Feynman parametrization. Of course, more complicated cases might still require more elaborate techniques, such as Integration-By-Parts reduction (IBP) [88,89], differential equations [90][91][92][93][94][95], sector decomposition [96][97][98] or Mellin-Barnes representation [99][100][101][102]. If the quantity one wants to calculate depends on the phase-space integration over a squared matrix element (possibly multiplied with other functions), it is common to do the evaluation using numerical methods. Unless the final state involves at most 2 or 3 legs, analytic results are very difficult to obtain, irrespective of whether one calculates the phase-space integrals directly or employs special methods such as reverse unitarity [103,104].
The obvious difference between the existing approaches that rely on private codes and FeynOnium is not only that our codes are public but also that we are explicitly interested in providing the complete scripts required to reproduce a particular result. This should hopefully motivate other members of the EFT community to share their software tools and also make the EFT techniques more accessible to a broader audience, including students and researchers working in different areas of quantum field theory.

FeynCalc 9.3
Turning FeynCalc into a tool that could support both relativistic and nonrelativistic calculations on the same footing was a challenging endeavor, both technically and conceptually. The reason is that FeynCalc was originally created to work with manifestly Lorentz covariant quantities, therefore it was not possible to design everything from scratch, but one had to ensure that the new features nicely fitted into the existing framework. One of the main goals was to preserve backward compatibility and to allow the user to employ already familiar functions such as Contract, ExpandScalarProduct or Dirac-Simplify without worrying whether the input contained nonrelativistic expressions or not. New methods were added only for manipulations that were not available or not required in the previous FeynCalc version, e.g. LorentzToCartesian for breaking manifest Lorentz covariance. This means that Cartesian tensors (just as Lorentz tensors) now belong to the most fundamental quantities that can be manipulated using FeynCalc.
While it is not our scope to give a full account of the implemented modifications, in the following we will describe the main design decisions that were taken to make FeynCalc useful for nonrelativistic calculations. This should hopefully help the reader to gain a better a feeling for the new abilities of the package.

Lorentz and Cartesian indices and vectors
The three fundamental FeynCalc objects used to manipulate 4-vectors, corresponding to scalar products, Levi-Civita tensors and Dirac matrices, are called Pair, Eps and DiracGamma, respectively. Essentially, Pair is a symmetric function with two slots. Each of the slots can accept two types of arguments, which are LorentzIndex (for Lorentz indices) and Momentum (for 4-momenta). Both of them also have two arguments. The first argument of LorentzIndex denotes the name of the corresponding index (e.g. µ, ν, ρ . . . ), while the first argument of Momentum specifies the name of the 4-momentum(e.g. p, q, l . . . ). The second argument of both functions fixes the spacetime dimension, which can be 4, D or D − 4. Moreover, the second argument is optional and when it is missing the spacetime dimension defaults to 4. Depending on the combination of its arguments Pair may represent a Lorentz vector (LorentzIndex and Momentum), a metric tensor (twice LorentzIndex) or a scalar product (twice Momentum). As far as DiracGamma is concerned, its first argument (LorentzIndex or Momentum) specifies whether we have a Dirac matrix with a free Lorentz index γ µ or a Feynman slash / p, while the optional second slot is used for setting the spacetime dimension. The representation of Levi-Civita tensors follows the same pattern, with Eps being a function that has four slots for LorentzIndex-or Momentum-type arguments.
This symbolic representation of Quantum Field Theory (QFT) quantities within Feyn-Calc is called internal or FCI-notation. In addition to that, FeynCalc is also equipped with an external or FCE-notation, which consists of convenient shortcuts that are useful for the manual input or when exporting FeynCalc results to other programs. FeynCalc functions usually output the results in the internal notation but accept the input in both notations. The routines for switching between the two notations are FeynCalcInternal (abbreviated with FCI) and FeynCalcExternal (abbreviated with FCE). For example, to input a D-dimensional vector p µ in the FCI-notation we need to write , while in the FCE-notation the same expression can be entered as FVD [p,µ]. A summary of FeynCalc symbols that represent tensors and matrices in both notations can be found in appendix B.
Automation of nonrelativistic calculations requires support for additional tensors that carry explicit temporal or spatial (Cartesian) indices. In particular, the code must be able to deal not only with manifestly Lorentz covariant quantities (e.g. p µ or l · q) but also with objects like p 0 , p i , l 0 q 0 or l · q. In FeynCalc 9.3 this has been achieved by extending the internal notation with the following symbols: CartesianPair, CartesianMomentum, CartesianIndex, TemporalPair, TemporalMomentum and PauliSigma.
The first three are conceptually similar to the above-mentioned Pair, Momentum and LorentzIndex. For example, CartesianPair is a special pairing that accepts Cartesian-Momentum or CartesianIndex as arguments for its two slots and can be used to represent 3-vectors, Cartesian scalar products or Kronecker deltas. It is important to stress that when CartesianMomentum and CartesianIndex have no second argument, their default dimension is 3. For calculations in dimensional regularization one should write them as D−1 dimensional quantities. This is because in FeynCalc every Cartesian tensor is always understood to be the spatial piece of the corresponding Lorentz tensor. Thus, if such a Lorentz tensor (e.g. 4-vector) lives in D dimensions, its spatial (a 3-vector) component must be a D − 1-dimensional object.
As the name already suggests, the main purpose of introducing TemporalPair and TemporalMomentum is to have a suitable representation for the temporal components of 4-vectors. Notice that we do not require a new dedicated symbol for the 0th index of a tensor, since it can be written using the already existing ExplicitLorentzIndex[0]. Moreover, TemporalMomentum has only one argument that denotes the original 4-vector. This simply reflects the fact that in dimensional regularization the temporal component of a 4-vector still remains a 1-dimensional object, while its spatial components are analytically continued to D − 1 dimensions.
CartesianPair and CartesianIndex are equally valid arguments of DiracGamma. In this case they obviously represent a Dirac matrix contracted to a 3-vector (e.g. γ i p i ) or a Dirac matrix with a Cartesian index (e.g. γ i ). The same also applies for Levi-Civita tensors. For example, an Eps with three distinct CartesianIndex arguments stands for ijk . Symbolic Pauli matrices σ µ and σ i are available via PauliSigma, which (similar to DiracGamma) can represent a Pauli matrix with a Lorentz or a Cartesian index as well as a Pauli matrix contracted to a Lorentz or a Cartesian vector. Here again we would like to refer to appendix B for the list of available quantities and the commands to enter them in FeynCalc.

Upper and lower indices
Many software frameworks for automatic QFT calculations do not explicitly distinguish between upper (contravariant) and lower (covariant) Lorentz indices. This simplification does not introduce any ambiguities, provided that all input expressions obey Einstein's summing convention and are written in a Lorentz covariant fashion. Given that every pair of Lorentz indices appearing in a single term is understood to be contracted with each other, it is clear that one of the indices must appear upstairs, and the other one downstairs. For example, in it is irrelevant whether FV[p,µ] stands for p µ or p µ , as long as µ is understood to be a dummy index. When dealing with free (i.e. uncontracted) indices, it is enough to know that in a manifestly Lorentz covariant expression the position of the index on one side of the equation must match its position on the other side. Consider e.g. the symbolic expression This can be interpreted asū as well asū Notwithstanding the ambiguity of this symbolic notation, it does not lead to inconsistencies so that the calculations will always return sensible results. Moreover, once we mentally fix the positions of the free indices in the input expression, manifest Lorentz covariance guarantees that these positions will not change in the course of all intermediate symbolic manipulations and will be preserved in the final result. While such a "mixing" of covariant and contravariant indices may seem aesthetically unpleasant, this trick greatly helps to improve the performance of symbolic codes, which is especially important when working with very large expressions. Things become more complicated once we want to handle expressions that contain both Lorentz and Cartesian tensors. Depending on the metric signature, moving a Cartesian or a temporal index into an opposite position may introduce a minus sign. For example, for Furthermore, such indices are not constrained to appear in the same position on both sides of an equation, so that expressions like 2. In a contraction of two Lorentz indices it is understood that one of them is upstairs and the other is downstairs.
3. In a contraction of two Cartesian indices, both indices are understood to be upper indices.

A free Lorentz or Cartesian index is always understood to be an upper index.
While the first two rules merely formalize something that was always implicitly assumed in FeynCalc calculations, the last two rules are new. The third rule was never required before, since earlier versions of FeynCalc could not deal with Cartesian tensors. The fourth rule helps to avoid ambiguities when interpreting expressions with free indices. Let us briefly illustrate how, by applying the above rules, we can interpret different FeynCalc expressions in a sensible way Notice that our notation also applies to tensors that carry both Lorentz and Cartesian indices. Such quantities often arise at different stages of nonrelativistic calculations and are, therefore, fully supported in FeynCalc 9.3. For example, where we used the FCI notation to write down a metric tensor with mixed indices, since such an object has no corresponding FCE shortcut.

Nonstandard integrals
FeynCalc is very often employed as a convenient tool for symbolic manipulations of loop integrals, especially at 1-loop. Integrals with only one loop momentum and standard 1/(p 2 − m 2 )-type propagators can be conveniently handled using the Passarino-Veltman reduction technique [105], which is available in FeynCalc since version 1.0. Indeed, tensor reduction and the subsequent analytic or numerical evaluation of the resulting Passarino-Veltman functions is sufficient for a large class of 1-loop calculations in the Standard Model and its extensions. Unfortunately, such methods often turn out to be inadequate when EFTs come into play. For example, eikonal propagators, as they appear in HQET or SCET, cannot be handled by the routines implemented in FeynCalc 9.2. The same is also true for Euclidean and Cartesian integrals as well as integrals involving temporal components of 4-vectors.
On the one hand, it is difficult to find a good strategy for treating such "nonstandard" integrals in FeynCalc in a generic fashion. As far as tensor reduction is concerned, even at 1-loop such integrals often cannot be directly rewritten in terms of scalar integrals with unit numerators. In the lack of a universal basis similar to the Passarino-Veltman functions, 2 the evaluation of the corresponding master integrals often proceeds on a case-by-case basis.
On the other hand, some algebraic manipulations that are needed in the course of the Passarino-Veltman reduction turn out to be applicable to almost all kinds of loop integrals. For example, partial fractioning and minimal tensor reduction to remove loop momenta with uncontracted indices can be straightforwardly applied to Cartesian and eikonal loop integrals such as .
The very first step in making FeynCalc useful for such calculations is to introduce new symbols to represent various nonstandard propagators. In the external notation this can be achieved by adding only 3 new shortcuts, which allow to cover a broad range of nonstandard loop integrals. These are SFAD (StandardFeynAmpDenominator), CFAD (CartesianFeynAmpDenominator) and GFAD (GenericFeynAmpDenominator), which stand for covariant, Cartesian or generic propagators respectively. Of course, for compatibility reasons the original symbol FAD (FeynAmpDenominator) has been kept in FeynCalc as it has been employed there since the very beginning. For the sake of clarity, table 1 summarizes all types of propagators that can be entered using the four above-mentioned shortcuts.

Shortcut in FeynCalc
Meaning Table 1: Implementation of the new propagator types using StandardFeynAmpDenominator, CartesianFeynAmpDenominator and GenericFeynAmpDenominator. Here x can be an almost arbitrary function of loop-momentum dependent scalar products.
While the old FAD covers only a small fraction of propagators that are possible with the new SFAD, the former is still somewhat better integrated into FeynCalc than the latter. This mainly concerns the use of the function FeynAmpDenominatorSimplify for detecting scaleless integrals and finding useful loop momentum shifts. These differences will be gradually eliminated in the future versions of the package, where the treatment of the new integral types will become more mature.
The syntax of SFAD may seem cumbersome at the first sight, but these inconveniences are more than compensated by the great flexibility encoded in this shortcut: Both quadratic and linear propagators are covered and the signs in front of the mass term m 2 and the causality parameter iη can be chosen freely. Furthermore, some common propagator types can be entered faster as follows Notice that in the case of massless eikonal propagators FeynCalc takes special care to preserve the sign of iη. Rewriting of propagators as in where the causality parameter switches its sign, is explicitly avoided. CFAD can be regarded as the Cartesian counterpart of SFAD with the main difference being that the default signs of m 2 and iη are opposite to that of SFAD e.g.
Apart from this characteristic feature, CFAD has virtually the same syntax as SFAD and can be used to enter different types of Cartesian integrals. Last but not least, one should also mention GFAD that acts as a generic placeholder for entering integrals that cannot be represented using SFADs and CFADs alone. For example, the singlet propagator in pNRQCD (cf. appendix C.2) is a quantity that explicitly depends on the temporal component of a 4-momentum flowing through the corresponding line. Hence, we can write its denominator as Since GFAD objects may represent almost arbitrary loop-momentum dependent denominators, FeynCalc will usually abstain from applying any kind of loop momentum shifts to integrals containing such propagators. This means that computations involving GFADs will require significantly more user intervention at intermediate steps than those relying on the simpler but less versatile SFADs and CFADs. It is therefore advisable not to introduce GFADs unless absolutely necessary. Having said that, we would also like to stress that partial fractioning and tensor reduction are nonetheless available also for integrals containing GFAD propagators. An important limitation that the users should be aware of concerns types of integrals that can be manipulated using the built-in functions. Internally, FeynCalc always classifies input integrals into three possible categories: 1. Integrals in which loop momenta appear solely as 4-vectors, meaning that such expressions enjoy manifest Lorentz covariance.
2. Manifestly noncovariant integrals where each integration measure is split into temporal and spatial components e.g. as in (4.17) where the integrand f (k 0 , k) explicitly depends on temporal and spatial components of the loop momentum k.
3. Integrals that are mixtures of covariant and noncovariant quantities e.g. as in where x is some c-number and p is an external 3-momentum.
FeynCalc can readily handle integrals of type 1 or 2, but the "mixed" integrals of type 3 cannot be processed straightforwardly. This is because the underlying code heavily relies on working with linearly independent scalar products involving loop momenta. However, it is hardly possible to guarantee linear independence once 3-momenta and 4-momenta are allowed to appear in the same integral. For example, one might encounter zeros in the form k · q + k · q, with q = (0, q) T , (4.19) which would remain undetected and hence lead to potentially disastrous consequences towards the end of the computation. Owing to the fact that integrals similar to those in eq. (4.18) do arise in many NREFT calculations, it is important to clarify how to handle them properly. Here we propose three different strategies depending on the form of the involved integrals and the expected difficulties in calculating the resulting master integrals.
In some cases it might be possible to recast a mixed integral into a form that is manifestly Lorentz covariant i.e. to convert a type 3 integral into a type 1 integral. As far as numerators are concerned, we can always introduce auxiliary vectors such as n = (1, 0, 0, 0) T and v = (p 0 , 0, 0, 0) T to have for a loop momentum k and an arbitrary external 3-momentum p. A 3-momentum vector with a free index can be written in a covariant fashion at the cost of introducing a metric tensor with Lorentz and Cartesian indices, e.g. as in However, when applied to denominators, these methods may produce inconvenient propagators such as 1 and alike. Furthermore, introducing too many auxiliary vectors will likely make the integral more complicated than it really is. This is why, in general, this approach is not always applicable or even sensible. The other two strategies require us to convert a mixed integral into a type 2 integral first. This procedure is straightforward and unambiguous, since any scalar product of two 4-vectors can be always decomposed into its spatial and temporal components as in (4.23) This means that each integration over a D-dimensional loop momentum k splits into a 1-dimensional integration over the temporal component k 0 and a D − 1-dimensional integration over the spatial comment k.
The second strategy would be to tensor reduce and partial fraction the k-integrals, whereas k 0 will be regarded as an external parameter. The resulting integrals (that still depend both on k and k 0 ) are then declared to be master integrals.
When employing the third strategy we would, on the contrary, integrate over k 0 first, ending up with pure Cartesian k-integrals. Those integrals may explicitly depend not only on the 3-vector k and its scalar products with external momenta, but also on the magnitudes of the 3-vectors, such as |k| and |k − p|. Furthermore, the k 0 -integration requires some care in applying the residue theorem and picking up the correct poles. Nevertheless, this procedure often leads to fewer and simpler master integrals as compared to the case where each master integral must be integrated in k and k 0 .
In principle, FeynCalc can be useful in all 3 scenarios of dealing with mixed integrals, but the level of automation will vary substantially. For example, integrations over the temporal components of loop momenta have to be performed by hand, since such a procedure is too difficult to automatize in full generality.
Another important restriction in the handling of tensor integrals is the requirement that those should not contain vanishing Gram determinants. Although this issue seems to be rarely discussed outside of the context of the Passarino-Veltman functions, the breakdown of naive tensor reduction for integrals with zero Gram determinants can, in principle, occur in all kinds of tensor integrals. For example, the result of the tensor reduction of the 3-point function is proportional to the inverse of the Gram determinant 4((p · q) 2 − p 2 q 2 ). Therefore, a naive attempt to tensor reduce this integral at the special kinematic point p · q = p 2 = q 2 will inevitably fail. In general, it is well known (cf. e.g. [106,107]) that many of such cases can be worked around by considering a larger nonsingular system of linear equations and extracting necessary relations to reduce the original integral. However, as of now, such procedures are not yet implemented in FeynCalc.
More details on practical manipulations of nonstandard integrals in FeynCalc 9.3 can be found in section 5.5.

Installation
The FeynOnium projects consists of two components: the recently released FeynCalc 9.3 [32] that can be used for nonrelativistic calculations and a homonymous add-on that is dedicated to NREFTs. Both FeynCalc and the add-on are open source 3 with the source code hosted on GitHub. 4 FeynCalc requires at least Mathematica 8 or later, while the add-on runs on top of FeynCalc 9.3 or later. The most convenient way to setup the whole framework is to use the automatic online installer. The FeynCalc installer can be invoked by running the following code in a new Mathematica session Although not strictly necessary, it is also recommended to install the FeynHelpers addon [108], which provides convenient and easy-to-use interfaces to other tools for evaluating Passarino-Veltman functions and performing IBP reductions of loop integrals. Last but not least, one should also consider downloading, FeynRules 5 as it is used to create FeynArts model files that are employed in some of the FeynCalc and FeynOnium examples.
Notice that an add-on can be activated only during the loading of FeynCalc. This is why before loading the package the names of the add-ons (as strings) must be specified in a list assigned to the global variable $LoadAddOns. For example, to use FeynOnium and FeynHelpers one should run at the very beginning of a Mathematica session.

Basic nonrelativistic calculations
Most standard FeynCalc routines for amplitude manipulations such as Contract, Uncontract, ExpandScalarProduct, MomentumCombine, ComplexConjugate etc. are directly applicable to expressions containing noncovariant quantities. Therefore, everyone who at least knows how to use FeynCalc for tree-level calculations should have no difficulties to master the new nonrelativistic capabilities of the package. The basic Cartesian tensors required for nonrelativistic studies are 3-vectors (e.g. p i abbreviated as CV[p,i]), Kronecker deltas (e.g. δ ij abbreviated as KD[i,j]), scalar products of two 3-vectors (e.g. p · q abbreviated as CSP[p,q]) and Cartesian Levi-Civita tensors (e.g. ijk abbreviated as CLC[i,j,k]). Notice that all Cartesian vectors are typeset bold, with 3-dimensional vectors having a bar and D − 4-dimensional vectors a hat. The vectors without a bar or a hat live in D − 1 dimensions. This agrees with the existing FeynCalc typesetting of 4-vectors that follows [109]. Explicitly, we have The same notation applies also to Dirac and Pauli matrices. The shortcuts CV, KD, CSP and CLC correspond to 3-dimensional quantities. Their D−1dimensional versions are obtained by attaching a D to the corresponding shortcut, e.g. as in CVD or KDD. Attaching an E yields the respective D − 4-dimensional symbol, e.g. CSE. 6 In this respect the new Cartesian tensors behave in the same way as the existing Lorentz quantities.
The crucial task of manifestly breaking Lorentz covariance of tensors and matrices can be accomplished using LorentzToCartesian. This function rewrites the occurring tensors and matrices with Lorentz indices in terms of their temporal and spatial components as in which is important e.g. when doing an amplitude-level matching between relativistic and nonrelativistic theories. For example, we can write Here the dollar sign indicates an index contraction between the 3-vector p i and the metric tensor with mixed indices g iµ . In the internal representation (currently there are no FCEshortcuts for such objects) we have At this stage one would often like to assign some specific values to the spatial and temporal components of 4-vectors and scalar products. This can can be done via direct assignments as in where we set n 0 = 1 and n · p = 0. As usual, these assignments can be removed via FCClearScalarProducts. One can also exploit additional simplifications by extracting the magnitude (which could be e.g. an expansion parameter in an EFT calculation) of a 3-vector or specifying relations between spatial components of some 4-vectors as in To this end we can assign values directly to a particular CartesianMomentum and let FeynCalc know that some symbols (e.g. |k|) are scalars and hence can be pulled out of expressions that denote vector contractions. The latter is done using the Datatype mechanism, where the corresponding symbols are defined to have datatype FCVariable. For example, the relations given in eq. (5.5) can be implemented via where we also account for the fact that the scalar product of a unit vector with itself is unity. In this case expressions such as (k + p) 2 or ijk k i p j l k can be directly rewritten as

In[7]:= {CSP[k + p], CLC[][k, p, l]} // ExpandScalarProduct
If one would like to differentiate with respect to a 3-vector, the corresponding routine is called ThreeDivergence. It works in exactly the same way as its 4-dimensional analogue FourDivergence e.g.

In[8]:= ThreeDivergence[1/(CSP[p, q] + a) (b + CSP[p]), CV[p, i]]
Let us now show how the introduced machinery can be employed in real-life calculations. To this end we can reproduce the value of the matching coefficient which enters the O(α 0 s v 2 ) differential production cross-section for e + (l 1 )+e − (l 2 ) → χ c 0 (P )+ γ(k) in the NRQCD factorization formalism [110][111][112]. Here s stands for the collision energy in the center of mass frame, while the summation sign implies that we must average over the polarizations of the leptons and sum over the polarizations of the photon. Explicit values of the short-distance coefficients c J=0 1 and c J=0 3 are given in [112]: where r = 4m 2 /s, ε * i (k) denotes the polarization 3-vector of the external photon and v(l 2 )γ i u(l 1 ) stands for the spatial piece of the leptonic current. The kinematics is chosen in such a way that We would like to evaluate the photon polarization sum by introducing an auxiliary vector n µ = (1, 0, 0, 0) T so that only physical degrees of freedom (transverse polarizations) are taken into account. Using FeynCalc, expressions similar to eq. (5.6) can be computed as follows. First of all, we need to specify all kinematic constraints which agrees with eq. (58b) in [112].

Dirac algebra
Since FeynCalc 9.3 all routines related to Dirac algebra support manipulations of Dirac matrices with temporal or spatial indices. This allows the user to evaluate very generic noncovariant expressions involving Dirac matrices such as So far, Euclidean Dirac matrices are not yet supported, but given enough interest from the user side they might be added in the future.

Pauli algebra
Pauli matrices are a completely new class of algebraic objects introduced in FeynCalc 9.3 for the first time. For the sake of consistency and user convenience, their handling was modeled after the existing implementation of the Dirac algebra. Therefore, it should not come as a surprise that FeynCalc is equipped with routines called PauliSimplify, PauliTrace and PauliOrder.

If a chain of 3-dimensional Pauli matrices contains repeated Cartesian indices or contractions with identical 3-vectors as in
If it is necessary to reduce the number of matrices in a chain to at most one by repeatedly applying the relation When doing loop calculations in dimensional regularization it becomes necessary to extend the definition of Pauli matrices to D − 1 dimensions. It is well known that the anticommutator of two Pauli matrices can be consistently generalized to Then, using eq. (5.16) one can derive relations for eliminating pairs of indices and vectors in a chain of Pauli matrices in D − 1 dimensions. The same also applies for traces of an even number of matrices. A collection of such formulas can be found e.g. in [113]. Yet the commutation relation of 3-dimensional Pauli matrices becomes ambiguous in dimensional regularization, as the Levi-Civita tensor ijk is intrinsically a 3-dimensional quantity. Related issues are well known to the practitioners and a valuable discussion of this topic can be found in [113]. In general, traces of odd numbers of Pauli matrices do not naively generalize to D − 1 dimensions. Similar issues arise when trying to reduce products of Pauli matrices (e.g. σ i σ j ⊗ σ i σ j ) to a finite 3-dimensional basis (e.g. 1 ⊗ 1 and σ i ⊗ σ i ). Applying eq. (5.14) or any other projection method will also generate contributions that vanish in the limit D → 4, the so-called evanescent operators. It is worth noting that evanescent operators multiplied by poles in 1/ε produce finite contributions to the final results and in general require dedicated treatments [114,115]. Different prescriptions for dealing with Pauli matrices in D − 1 dimensions can be found in the literature [113,[116][117][118] and it is important to be aware of these issues to avoid inconsistencies.
As far as FeynCalc is concerned, the precise treatment of D − 1-dimensional Pauli matrices can be specified via FCSetPauliSigmaScheme[]. The default value is "None", meaning that only eq. (5.16) is used to simplify chains of Pauli matrices, while σ-odd traces are left unevaluated. In this way the returned results are always unambiguous.
In order to obtain more compact (but also scheme-dependent) expressions, one may want to specify a prescription for evaluating the remaining traces and reducing chains of Pauli matrices to a minimal basis. In the current version of FeynCalc the user can evaluate The occurring products of D − 1-dimensional Levi-Civita tensors are then calculated using where all Kronecker deltas are defined in D − 1 dimensions. In particular, we have We are looking forward to the feedback and suggestions from the NREFT community to implement more useful prescriptions in future iterations of the framework.

Loop calculations
The two main tools provided in FeynCalc for dealing with loop integrals are tensor reduction (via TID and FCMultiLoopTID) and partial fractioning (via ApartFF). Both operations employ FeynAmpDenominatorSimplify (also abbreviated as FDS) for recognizing vanishing integrals and applying suitable loop momentum shifts. As has already been mentioned in section 4.3, FDS mainly relies on heuristics which may not work so well with nonstandard propagators, thus missing some obvious simplifications and not setting scaleless integrals to zero.
Let us discuss tensor reduction. If we are dealing with 1-loop tensor integrals that can be reduced to scalar integrals with unit numerators, it is advantageous to employ TID. Such a reduction is always possible for purely Lorentz or Cartesian integrals with quadratic propagators, but the support for Cartesian integrals is a new feature of FeynCalc 9.3. For simplicity, we can consider a massless Cartesian rank 2 tensor integral with one external momentum Of course, more complicated integrals are also possible, as there are no formal limitations on the tensor rank and the number of external legs that can be processed by TID. The only practical limitation is the degrading performance when handling very complicated tensor integrals.
When confronted with mixed integrals (cf. section 4.3), TID can often automatically perform tensor reduction with respect to the spatial part of the loop momentum. This is certainly true for integrals such as

In[21]:= TID[CFAD[k, k -p] FVD[k, i], k]
where TID essentially applies tricks described in section 4.3. Notice that the result for the integral in eq. (5.22) still contains a scaleless integral d D k/(k 2 − (k 0 ) 2 ). This is an example of the difficulty of enhancing FDS with good heuristics for nonstandard integrals, especially when integrations in the temporal and spatial components of loop momenta must be treated separately.
Tensor reductions of integrals that are expected to contain irreducible denominators should be done using FCMultiLoopTID. Such denominators constitute a common feature of multiloop integrals, but they often arise already at 1-loop once propagators different from quadratic ones come into play. This is also the main reason why TID refuses to handle integrals with eikonal propagators: In such cases there is simply no guarantee that the reduction to integrals with unit numerators can succeed. FCMultiLoopTID is not affected by this problem, because it only considers loop momenta with free indices or those contracted to Dirac or Pauli matrices, Levi-Civita tensors and polarization vectors. For example, FCMultiLoopTID does not regard as a tensor integral and will therefore leave it unchanged. On the contrary, in the case of 2 (5.25) the function will uncontract the scalar product of k and the Pauli matrix σ, producing a rank 1 tensor integral that will be subsequently reduced

In[22]:= FCMultiLoopTID[CSISD[k] CFAD[k, k -p], k]
If we know that uncontracting particular scalar products of loop momenta with other vectors may lead to a simpler result, we can use the option Uncontract to specify those vectors explicitly. One of such examples is we can nonetheless achieve the desired reduction.
The next thing we would like to discuss is partial fractioning. In FeynCalc 9.3 ApartFF has been extended to support the newly introduced nonstandard propagators, thus making it possible to handle many nontrivial cases such as However, we also observed that due to the specifics of ApartFF, some desirable decompositions cannot be obtained automatically. This mainly concerns integrals with propagators that do not form an overdetermined basis. Since ApartFF is applicable only to cases with overdetermined propagator bases (cf. section 3.3 in [21] for more details), it would normally ignore such integrals altogether. The trick to overcome this behavior is to multiply the corresponding integral by unity i.e. by a suitable propagator and its inverse. Given that the product of the original integral and the extra propagator contains an overdetermined basis of propagators, we may freely subject it to partial fractioning. At the end, multiplying back the so-obtained result with the inverse of the auxiliary propagator ensures that the final result is equivalent to the original expression.
For definiteness, let us consider the integral where we would like to trade the numerator k · p for k 2 . We cannot achieve this neither with FCMultiLoopTID nor using the standard mode of ApartFF. This is why we extended the syntax of ApartFF to support the above-mentioned trick. When the second argument of the function is not a list, it is interpreted as the inverse of the auxiliary denominator that has already been added to the integral in the first argument. After having carried out such partial reduction, it is usually advisable to run ApartFF again (this time in the standard mode), to simplify the product of the intermediate result with the inverse denominator. In the case of the integral in eq. (5.28) we obviously need to introduce the unity as k 2 /k 2 = 1. Therefore, we multiply eq. (5.28) by 1/k 2 (written as CFAD[k]) and put k 2 (as CSPD[k]) into the second slot of ApartFF. As far as the nonstandard propagator 1/|k| is concerned, we can write it as a GFAD with √ k 2 . Owing to the abundance of such propagators in nonrelativistic calculations, we deliberately added support for square roots of Cartesian scalar products to FeynCalc. Putting everything together, we have
Using the manipulations described in this section it should be possible to handle a wide range of (NR)EFT calculations, at least at 1-loop level.

FeynOnium
The FeynOnium extension builds upon the new symbols and routines introduced in the previous sections. It provides tools that help to streamline NREFT calculations by reducing the amount of code that needs to be written from scratch. Furthermore, FeynOnium includes a number of worked out examples that explicitly reproduce selected NREFT results from the literature. This should not only help practitioners to quickly master the new framework but also lower the entry barrier for students and researchers from other branches of particle physics who would like to familiarize themselves with NREFT techniques.
Most FeynOnium functions tend to produce rather large output expressions, which are best viewed and processed within Mathematica. Therefore, we prefer not to clutter this section by copying long code samples. Instead, we would like to explain the conceptual ideas behind those routines, making it clear where and why they should be used in practice. For explicit usage examples we refer to the Mathematica notebook accompanying this publication and scripts reproducing physical results that are provided together with the program.
In a matching calculation between a relativistic and a nonrelativistic theory with fermionic degrees of freedom it is often useful to rewrite Dirac spinor chains in terms of Pauli matrices and Pauli spinors. To this end FeynOnium provides two special functions. FMSpinorChainExplicit2 merely rewrites the chains using the Dirac representation of the Dirac matrices without making any additional assumptions about the underlying process. In contrast, FMSpinorChainExplicit is specifically tailored for studying production or decay processes of heavy fermions in the rest or laboratory frame using Jacobi momenta. It implements the threshold expansion method from [119] in 4 dimensions, where the small relative momentum between the two fermions in the rest frame can be used as an expansion parameter. The generalization to 3-body problems first derived in [112] is also implemented. However, prior to applying FMSpinorChainExplicit it is necessary to perform an SPVAT (scalar, pseudoscalar, vector, axial-vector, tensor) decomposition of all Dirac chains using DiracReduce, convert the obtained spinor chains to a special notation via FMToStandardSpinorChains and employ LorentzToCartesian to break the manifest Lorentz covariance. In order to disentangle contributions from different angular momentum components J in an amplitude one may want to explicitly project out the corresponding components of suitable tensors as in The routine FMCartesianTensorDecomposition encodes projections with J = 0, 1 and 2 for 3-dimensional tensors up to rank 5 and can be easily extended to contain more J-values and higher rank tensors.
Another issue that regularly arises in complex nonrelativistic calculations are spurious terms that vanish by the virtue of the 3-dimensional Schouten identity ijk p l − jkl p i + kli p j − lij p k = 0, (5.31) where p is an arbitrary Cartesian vector. In general, it is very difficult to apply this identity in a systematic way, which is why FeynOnium features a tool that facilitates this task. 7 FMCartesianSchoutenBruteForce tries out all possible combinations that can be formed out of the given list of Cartesian vectors and checks if this helps to reduce the number of terms in the expression. Although this approach may seem hopeless at first sight, in practice we observe that it works surprisingly well, eliminating most of the spurious terms after some number of iterations. The use of covariant projectors for heavy nonrelativistic systems introduced in [120] can be automatized via FMInsertCovariantProjector. Production and decay projectors for spin singlet/triplet and color singlet/octet states can be thus applied straightforwardly.
Last but not least, we also implemented Feynman rules for pNRQCD vertices in the weak-coupling regime at order r in the static limit as given in figure 5 of [8]. In the lack of a convenient way to generate pNRQCD Feynman diagrams automatically, 8 our implementation should significantly facilitate the tedious task of entering pNRQCD amplitudes by hand.

Examples
In order to show how various functions of

Euler-Heisenberg Lagrangian
Although a systematic investigation of the EFT approach did not commence before the 70s of the last century, one can find many examples of much earlier applications of these techniques. One of them is the Euler-Heisenberg (EH) Lagrangian [121], an EFT of QED devised to describe photon-photon scattering at energies much below the electron mass m e . The only degrees of freedom in this theory are low-energetic photons that interact with each other via 2n-photon vertices, with n ≥ 2. These vertices arise from considering 2nphoton scattering amplitudes in the full theory (QED) and integrating out the electrons. Since QED, unlike QCD, contains no tree-level gauge boson self-interactions, the matching starts at 1-loop. Effective vertices with an odd number of photons are forbidden by Furry's theorem [122]. At leading order in the 1/m e expansion, the EH Lagrangian reads This theory is often presented in the introductory lectures to EFT (cf. e.g. [5], [123]), but the computation of the matching coefficients is either completely omitted or left as an exercise. More technical details can be found in [124], yet the reader must still work out the missing steps on her or his own. Here we will follow the calculation of [124] and show how the matching coefficients (at leading order) can be determined in a semi-automatic fashion. On the QED side of the matching we need to consider the process with p 2 i = 0. It is sufficient to work with the forward scattering configuration p 1 = p 3 , p 2 = p 4 which leaves us with only two kinematic invariants: p 1 · p 2 and m 2 e . Then we can strip the QED amplitude of the polarization vectors and equate it to the corresponding amplitude in the EH EFT so that where c 1 and c 2 are the unknown matching coefficients. This tensor equation can be converted into a system of two scalar linear equations by contracting it with g µν g ρσ and g µρ g νσ . Then the task of determining c 1 and c 2 is reduced to the calculation of T µνρσ QED g µν g ρσ and T µνρσ QED g µρ g νσ expanded up to the third order in p 1 · p 2 around 0. To obtain the EFT amplitudes automatically, we need to create a FeynRules model of the EH Lagrangian at O(1/m 4 e ) and export it to FeynArts. Since the Lagrangian is Lorentz covariant, this can be done in a straightforward fashion. The corresponding model file EulerHeisenberg.fr is already included in FeynCalc 9.3 and can be converted into a FeynArts model by evaluating the script GenerateModelEulerHeisenberg.m. 9 After having generated the QED and EH EFT amplitudes via FeynArts we need to convert them to the FeynCalc notation which is done with FCFAConvert. Then, Contract and DiracSimplify are employed to perform the contractions of the Lorentz indices and to simplify the Dirac algebra, including the evaluation of the Dirac traces. The quantities T µνρσ QED g µν g ρσ and T µνρσ QED g µρ g νσ are first reduced to scalar 1-loop integral by the means of TID. Then, FIREBurn calls FIRE 5 [125] to eliminate propagators raised to integer powers using the IBP-reduction. The evaluation of the remaining scalar 1-loop integrals, including the expansion in p 1 · p 2 is handled by PaXEvaluate, a frontend to Package-X. As expected, the final result is free of UV and IR divergences, so that we can directly take the limit d → 4. Upon substituting all the contributions back into the system of linear equations, we obtain the known result where α is the fine structure constant. While one certainly could perform this calculation in a more efficient way using FORM and the C++ version of FIRE, the advantage of the presented approach is that it requires almost no familiarity with tools for automatic calculations and can be employed even by undergraduate students. On the other hand, a recent work [126] that explored higher order operators in the EH Lagrangian and its QCD counterpart using FeynCalc and FeynHelpers clearly shows that these tools are very useful also in real research.

Heavy Baryon Effective Theory
FeynCalc's new ability to manipulate loop integrals with eikonal propagators can be handy even in comparably simple cases, such as the 1-loop correction to the heavy nucleon propagator in baryonic χPT [127]. Following [128], we need to evaluate which corresponds to a Feynman diagram with a nucleon emitting and absorbing a pion of mass M . Even though this calculation can be certainly done by pen and paper, with 9 The script is located in Examples/FeynRules/EulerHeisenberg inside the FeynCalc directory.

FeynCalc 9.3 it is effectively a one-liner that consists of applying PauliSimplify and
FCMultiLoopTID to eq. (6.5) and readily yields the two master integrals 1/(k 2 − M 2 ) and 1/[v · (r − k) (k 2 − M 2 )] in agreement with the literature.

Dimension six 4-fermion operators in NRQCD (unequal mass case)
Matching coefficients that multiply NRQCD dimension six 4-fermion operators in the unequal mass case were originally obtained in [117]. The corresponding operators are given by where m 1 (m 2 ) denotes the mass of a heavy quark (antiquark). The matching coefficients are determined by the hard momentum region (i.e. loop momenta of order m 1 , m 2 ) in the QCD box diagrams contributing to It is convenient to rewrite the momenta p i as To extract the values of d ss , d sv , d vs and d vv , it is necessary to expand the QCD amplitudes at 0th order in the small relative momenta q and q . It is convenient to work in the center of mass frame, so that setting q and q to zero is equivalent to setting p 1 = p 2 = p 3 = p 4 = 0. (6.11) In this case all scalar products between the 4-vectors p i can be expressed through polynomials in m 1 and m 2 , e.g.
Having obtained the amplitude from FeynArts we carry out the tensor integral reduction as well as Dirac and color algebra simplifications. Employing the "naive" scheme for dealing with Pauli matrices in D-dimensions we can readily rewrite all Dirac structures in terms of ξ † σ i ξ and η † σ i η. By using Contract with the option EpsContract set to False we actively prevent contractions of products of Levi-Civita tensors and can therefore implement the prescription of [117] via a replacement rule. Explicitly, this amounts to using ijk ijk = (D − 2)δ kk . (6.13) The loop integral structure of the amplitude has already been rewritten in terms of Passarino-Veltman functions which can be directly evaluated via PaXEvaluateUVIRSplit for D = 4 − 2ε and expanded around ε = 0. Using following Fierz identities for color matrices [7] T , we can easily separate color singlet and color octet contributions from each other. The separation into spin singlet and spin triplet pieces is even simpler, the former being proportional to ξ † ξ or η † η and the latter to ξ † σ i ξ or η † σ i η respectively.
Thus we finally recover the known results from the literature given by [117] (confirmed also in [129]) where α s is the strong coupling constant.

J/ψ → 3γ decay in NRQCD
The LO (both in velocity and α s ) NRQCD prediction for the decay J/ψ → 3γ (or Υ(1S) → 3γ) can be extracted by adapting the corresponding calculation for orthopositronium [130]. Nonetheless, it is also instructive to explicitly repeat this calculation by matching the QCD tree-level amplitude Q(p 1 ) +Q(p 2 ) → γ(k 1 ) + γ(k 2 ) + γ(k 3 ) (6.16) to NRQCD. The kinematics is and It is also convenient to parametrize |k 1 | and |k 2 | as where x 2 ranges from 0 to 1 − x 1 , while x 1 will be eventually integrated from 0 to 1. It is sufficient to expand the 6 QCD diagrams to 0th order in |q| and, after switching to Pauli matrices and spinors via LorentzToCartesian and FMSpinorChainExplicit2, we can square the QCD amplitude and sum over the polarizations of the photons. The angular integration of the 3-body phase space can be replaced by the J = 0 projection with respect to the unit vectorsk 1 andk 2 via FMCartesianTensorDecomposition. Upon integrating over x 1 and x 2 (here it can be done analytically using Mathematica's Integrate) and multiplying with the corresponding prefactor we obtain the total decay rate in QCD at LO given by where e Q is the fractional electric charge of the heavy quark Q. Comparing it to the corresponding perturbative NRQCD expression we correctly conclude that [7] Imf em ( 3 S 1 ) = 4 9 (π 2 − 9)α 3 e 6 Q . (6.23)

QQ → γγ decays in NRQCD
Let us consider at tree-level the QCD process with p 1 = 1 2 P + q, p 2 = 1 2 P − q, P · q = 0, (6.25) and p 2 1,2 = m 2 Q , k 2 1,2 = 0, (6.26) where we want to expand the amplitude in the relative momentum of the heavy quark pair, q, up to 4th order. To this end it is necessary to spell out all kinematic dependence on |q|, e.g. to specify q 0 = 0, p 0 1,2 = k 0 1,2 = E q , (6.27) with E q = q 2 + m 2 Q . This nonrelativistic expansion requires us not only to distinguish between spin singlet and spin triplet contributions but also to explicitly project components of the amplitude corresponding to the total angular momentum values J = 0, 1, 2. Of course, the J = 1 contribution must vanish due to the Landau-Yang theorem [131,132].
These steps can be directly performed with the aid of the FeynOnium functions FM-SpinorChainExplicit2, PauliSimplify and FMCartesianTensorDecomposition. Furthermore, by applying FMCartesianSchoutenBruteForce to the J = 2 contribution we can readily remove terms that vanish by Schouten's identity.
The so-obtained QCD amplitudes can be then matched to NRQCD, as it was done in [133], to obtain the matching coefficients relevant for the quarkonium decay processes η Q → γγ, χ Q0 → γγ and χ Q2 → γγ, where the heavy quark flavor Q can be c (for charmonia) or b (for bottomonia). Since we cannot generate NRQCD amplitudes automatically, they must be entered by hand. We write the NRQCD amplitudes up to order |q| 4 . In the next step, we square the NRQCD amplitudes and sum over the polarizations of the photons to arrive to the final heavy quarkonia decay rates in perturbative NRQCD. From there we can read off the matching coefficients of NRQCD decay operators that contribute through the leading heavy quarkonium Fock state |QQ . These are at LO (in agreement with [133]) Imf em ( 1 S 0 ) = α 2 e 4 Q π, (6.28a) Imf em ( 3 P 0 ) = 3α 2 e 4 Q π, (6.28c) Imf em ( 3 P 2 ) = 4 5 α 2 e 4 Q π, (6.28d) The definitions of operators multiplying these coefficients can be found in the appendix of [133].

Inclusive hadronic decays of P -wave quarkonia in NRQCD
The inclusive decay of χ QJ (Q = c or b) into light hadrons (LH) at LO in the velocity in the framework of NRQCD can be written as [7] Γ 29) where χ QJ |O 1 ( 3 P J )|χ QJ and χ QJ |O 8 ( 3 S 1 )|χ QJ denote NRQCD matrix elements. Here we would like to consider virtual next-to-leading order (NLO) loop corrections in the 2 gluon channel to Imf 1 ( 3 P 0,2 ), which were first calculated in [134] (with IR divergences regularized with a gluon mass) and later in [135] using the NRQCD formalism and explicitly distinguishing between UV and IR poles in dimensional regularization.
To extract the Born level contribution we need to consider three 10 tree-level diagrams describing the process Q(p 1 ) +Q(p 2 ) → g(k 1 ) + g(k 2 ) (6.30) in QCD, where all external particles are put on-shell. The momenta of the heavy quarks can be rewritten as where P is the total heavy quarkonium momentum and q is the relative momentum. Since we are not interested in relativistic corrections we can simplify the kinematics by setting Following [135] we project on spin-triplet P -wave states via a suitable spin-triplet color singlet covariant projector where A is the original amplitude, E αβ stands for the polarization of the quarkonium, and I c is a unit matrix in color space. The trace is understood to be taken over spinor and color indices. The trace of A can be implemented with FMInsertCovariantProjector, while the derivative with respect to q β can be obtained using FourDivergence. When squaring the amplitude A S=1,L=1 using ComplexConjugate, we need to sum over the polarizations of the quarkonia with different J-values using Doing so we recover (ε is the dimensional regularization parameter from D = 4 − 2ε) The calculation of the virtual corrections proceeds along the same lines as above, but is technically more challenging. We need to evaluate QCD 1-loop corrections to the process given in eq. (6.30). The results for the J = 0 and J = 2 contributions for each Feynman diagram are available in tables 2 and 3 of [135]. Contrary to the approach chosen in the original publication, we choose to evaluate loop integrals after and not before applying covariant projectors and expanding in q. As it has been observed in [136], in this case we do not encounter any Coulomb singularities. Furthermore, to speed up the calculation, we choose to reduce the number of integrals that need to be evaluated by using IBP reduction. This obscures however the distinction between UV and IR divergences in dimensional regularization, so that we denote both kind of poles with ε. Despite using solely Mathematica we can obtain the final result within half an hour on a modern laptop. Having assembled the virtual color singlet contributions to the decay of 3 P J quarkonia into two gluons from the evaluated 1-loop diagrams and the tree-level amplitude we find The result from [135] agrees with our expression upon setting ε UV = ε IR and dropping the Coulomb singularity.

Relativistic corrections to quarkonium light-cone distribution amplitudes
FeynOnium's capability of handling nonrelativistic objects can be combined with the built-in 4-dimensional algebra of FeynCalc to perform nonrelativistic expansion of QCD hadronic matrix elements that involve lightlike collinear momenta. In this section, we show how FeynOnium can be used to compute light-cone distribution amplitudes (LCDAs) of J P C = 1 −− heavy quarkonia in NRQCD to relative order v 4 accuracy at leading order in α s , which was first done in [137]. Following [137], we compute the LCDA for the J/ψ state with momentum P , which is defined by the matrix element 42) where ε is the polarization 4-vector of the J/ψ, and Q(x) is the QCD heavy-quark field, which is a four-component Dirac spinor field. The Wilson line (6.43) where P is the path ordering operator, ensures the gauge invariance of the nonlocal operator Q α (x). The light-cone vectors n andn are lightlike vectors that satisfy n ·n = 2.
To relative order v 4 accuracy, the J/ψ LCDA is given in the NRQCD factorization formalism by − J/ψ| · Q(x)|0 = nc n (x) m dn−3 ε · J/ψ|O n |0 , (6.44) where d n is the dimension of the NRQCD operator O n , andc n (x) are perturbative shortdistance coefficients. To relative order v 4 accuracy, the sum over n involves the NRQCD operators with J P C = 1 −− of dimensions up to 7 that are listed in eqs. (4) and (6) of [137]. Additionally, NRQCD operators with J P C = 1 +− of dimensions up to 7 can be found in appendix C of [137]. Our goal is to compute thec n (x) at leading order in α s . The coefficientsc n (x) can be determined from the matching conditions On the left-hand sides of eqs. (6.45), we project onto the J P C = 1 −− state, because unlike the NRQCD operators on the right-hand sides of eqs. (6.45), the operator Q α (x) does not have a definite J P C . Here, ε is the polarization vector of the perturbative QQ or QQg state; a state |QQg is a state made of a heavy quark (Q), a heavy antiquark (Q) and a gluon (g). While the NRQCD matrix elements on the right-hand sides of eqs. (6.45) can be computed straightforwardly at orders g 0 and g 1 , respectively [112], the calculation of the QCD matrix elements on the left-hand sides of eqs. (6.45) is much more involved. FeynOnium can provide significant simplifications of the calculation: the complicated algebraic manipulations that are needed to compute the QCD matrix elements on the left-hand sides of eqs. (6.45) can be done by using FeynOnium's built-in functions.
The initial step of the calculation involves the computation of the tree-level diagrams that contribute to QQ|Q α (x)|0 and QQg|Q α (x)|0 , and expanding in powers of the small relative 3-momenta of the Q,Q and g. The temporal and spatial components of the 4momenta that appear in the matrix elements are made explicit in the expressions using the LorentzToCartesian command. In this way complicated 4-dimensional tensors can be rewritten in terms of Cartesian tensors. The gamma matrices are now expressed in terms of Pauli matrices, and the spins of the heavy quark and the heavy antiquark combine into spin singlets and spin triplets. Then, the expansion in powers of the small momenta can be done by standard Mathematica commands.
The calculation of the left-hand sides of eqs. (6.45) is completed by projecting onto J = 1, C = −1 and P = −1. The projection onto C = −1 is straightforward: since the charge conjugation of the operator The projection onto P = −1 is also straightforward. The parity transform involves reversing the signs of momenta and gluon polarization 3-vectors, as well as the spins of the heavy quark and the heavy antiquark in singlet and triplet combinations. This can be done simply by using Mathematica's Replace command. Alternatively, instead of projecting onto P = −1, we can also include NRQCD operators with P = +1 on the right-hand sides of eqs. (6.45).
The projection onto J = 1 can be complicated, as it requires the application of the Cartesian tensor reduction algorithm developed in [138] to reduce the tensors of rank up to 5 built from 3-vectors to tensors of rank 1. This reduction can be done automatically with FeynOnium's built-in command FMCartesianTensorDecomposition.
Once QQ(J P C = 1 −− )|Q α (x)|0 and QQg(J P C = 1 −− )|Q α (x)|0 have been expressed in terms of 3-vectors, the short-distance coefficientsc n (x) can be obtained from eqs. (6.45). These results for thec n (x) can then be plugged into eq. (6.44) to obtain the quarkonium LCDAs. In [137], the J/ψ and Υ(nS) LCDAs were computed to relative order v 4 accuracy using FeynOnium, and the LCDAs were then used to compute Higgs boson decays to a heavy quarkonium plus a photon.

One-loop running of the chromoelectric dipole interaction in pNRQCD
To illustrate the practical usefulness of the nonrelativistic integral manipulations presented in section 5.5, we can reproduce the 1-loop renormalization group equations (RGEs) for the running of the matching coefficients V A (r) and V B (r) in weakly coupled pNRQCD (cf. eq. (2.3)) [31,139].
The calculation proceeds by first entering 1-loop amplitudes for the pNRQCD processes S → Og (for V A ) and O → Og (for V B ). The corresponding two diagrams are shown in figures 7 and 8 of ref. [31], respectively. The FeynOnium shortcuts for the pNRQCD Feynman rules render this procedure less error-prone and more convenient than copying expressions from a pen and paper calculation. We choose to handle mixed integrals such as d D k /[k 2 k 2 (p − k) 2 ] by explicitly integrating over k 0 and closing the contour below, so that we enclose the pole at k 0 = |k| − iη.
To apply this operation in FeynCalc we employ the function FCLoopExtract, which gives us a list of all loop integrals present in the expression. Then we implement the residue integration in form of a custom function that automatically generates a suitable replacement rule for each of the loop integrals. Upon substituting these results back into the original amplitude we end up with purely Cartesian integrals.
Our next step consists of employing FCMultiLoopTID for the tensor reduction and using ApartFF to enforce some custom partial fractionings. To this end we multiply the amplitude with k 2 /k 2 and |k|/|k| as explained in section 5.5. After this, we end up with 14 resulting 1-loop integrals that need to be calculated by hand. Some of them are obviously scaleless and can be put to zero immediately. Notice that to obtain the running we care only about UV singularities, so that all UV-finite parts can be discarded. This significantly simplifies the evaluation of the master integrals.
We find that both for V A and V B the pole of the first diagram contributing to the running is canceled by the pole of the second diagram. 11 This implies that at O(α s ) it holds that which reproduces the results of [139].

Summary
In this work we have presented a software named FeynOnium that works on top of Wolfram Mathematica and FeynCalc. The main purpose of FeynOnium is to facilitate the application of EFTs to various particle physics phenomena at tree-and 1-loop level by providing a large set of useful routines for tasks that typically arise in such calculations. One of the highlights of the package is the ability to perform calculations in nonrelativistic QFTs, a feature not readily available in other public codes. To achieve this, it was necessary to perform extensive modifications of the FeynCalc package, which is now capable of directly manipulating Cartesian tensors and integrals.
FeynOnium is open-source, publicly available, flexible and easy to use. The large number of included examples illustrates how this software can be used to quickly reproduce many well known EFT results from the literature, in particular in the frameworks of NRQCD, pNRQCD, HQET and ChPT.
Conceptually, the usage of FeynOnium is very similar to a pen and paper calculation, in the sense that everything can be organized as a sequence of simple operations (contractions, expansions, algebraic simplifications etc.) and all intermediate expressions are easily accessible for plausibility checks or comparisons to existing results.
FeynOnium should not be confused with other packages that target EFT practitioners but attempt to hide the technical side of the calculation by presenting to the user only the final results such as matching coefficients or cross sections. Our design philosophy was not to create another "black box", but a "toolbox" that naturally assumes the familiarity of the user with the EFT methods and provides her or him with means to investigate the relevant questions in a very flexible and convenient way.
We readily admit that performance-wise FeynOnium can hardly compete with codes that were specifically tailored and optimized for a particular calculation, especially if they are based on FORM. However, while such high performance private codes are usually accessible only to a very small subset of EFT practitioners, our package is available for everyone. This promotes good scientific practice by sharing tools that are beneficial for the whole EFT community and in particular by encouraging researchers from other branches to embrace the EFT techniques.

A Useful formulas for Pauli algebra
In this appendix, we collect some formulas for Pauli algebra calculations in 4-and Ddimensions that we find particularly useful and that can be readily implemented in a computer algebra system.
Chains of Pauli matrices in 4-or D-dimensions can be simplified as follows A trace of an even number of Pauli matrices in 4-or D-dimensions can be evaluated via the following recursive relation A trace of an odd number of Pauli matrices is not well defined in D-dimensions, but in 4-dimensions one can use

B Tensors and matrices in FeynCalc
For the sake of completeness, we list here the implementation of various tensors with Lorentz and Cartesian indices in the internal and external representations of FeynCalc. These are useful for interpreting FeynCalc results and writing codes that rely on the package.

C Derivation of Feynman rules in NREFTs
The extraction of Feynman rules from a given Lagrangian is an important ingredient of every particle physics calculation. Since every EFT Lagrangian formally contains an infinite number of operators, it is clearly not possible to derive those rules once and for all, as it can be done e.g. for QED or QCD. Instead, every EFT calculation that aims for higher precision in the expansion parameter(s) must take into account new operators that show up at that accuracy. On the one hand, the procedure of deriving Feynman rules for new operators can be automatized using dedicated software packages such as LanHEP [140,141], SARAH [142][143][144][145] or FeynRules. On the other hand, those tools mainly focus on Lorentz covariant theories and are therefore less useful for nonrelativistic calculations. For example, operators containing Pauli spinors, Cartesian tensors, spatial and temporal derivatives or (chromo)electric and (chromo)magnetic fields are not supported out-of-the-box. For this reason Feynman rules for NREFTs are still often derived by hand. Although the technicalities behind this procedure are certainly very familiar to the practitioners, they are rarely discussed at length in the literature. In the following we would like to treat this subject in a more detailed and pedagogical way, including explicit examples and useful recipes for practical calculations.
Path integral formalism and canonical field quantization are the two most popular methods for deriving Feynman rules. We choose to employ the latter procedure, owing to its conceptual simplicity and the fact that it is straightforward to automatize in almost any symbolic manipulation framework. A very concise description of the method can be found in the FeynRules manual [51], which we will follow here.
However, when dealing with NREFTs we also need to account for fields that directly annihilate the vacuum i.e. satisfy φ |0 = 0. This happens if a field that contains both particle and antiparticle components is transformed in such a way, that both components decouple from each other and are then treated as separate fields that create/annihilate a single particle/antiparticle. Therefore, we will slightly modify the rules from [51], making them applicable to relativistic and nonrelativistic theories alike. The main recipe for deriving the Feynman rule associated to a given operator O can be then summarized in the following 3 steps: 1. For fields that do not annihilate the vacuum i.e. φ |0 = 0 and 0| φ = 0, multiply O by creation operators for the fields from the right. For fields that annihilate the vacuum, multiply O by creation operators for φ from the right and by annihilation operators for φ † from the left. This ensures that the matrix element under consideration does not vanish. The fermion creators and annihilators should be ordered in a reversed way as compared to the ordering of the corresponding fermion fields in the operator.
2. Put the resulting expression between the vacuum states 0| and |0 . Move the creation (annihilation) operators to the left (right) where they annihilate 0| (|0 ).
3. Replace 0|0 by unity, remove the overall exponential and the external states (e.g. spinors and polarization vectors) and multiply the rest by i. Finally, reverse the sign of each momentum that stems from an annihilation operator that was multiplying O from the left, so that all momenta are incoming.
Conceptually, this technique is very similar to the calculation of matrix elements in the socalled "old fashioned perturbation theory" that was the main way of doing QFT calculations even before the invention of Feynman diagrams. Regarding the treatment of operators with field derivatives, the FeynRules approach is to pull the derivatives outside of the matrix element and apply them to the exponentials after the second step. Alternatively, one can also work out the (anti)commutation relations for fields with derivatives in advance and use them already during the second step. From our experience, the latter is often more convenient in pen and paper calculations, while what is done in FeynRules is naturally more useful for automatic codes. The main ingredient of the provided recipe is the process of moving creation and annihilation operators past field operators, where we need to apply suitable (anti)commutation relations. Once we have quantized the free part of our EFT in the operator formalism, these relations are obtained straightforwardly.

C.1 Feynman rules for NRQCD
For definiteness, let us start with NRQCD. The free part of the NRQCD Lagrangian at O(1/m) reads where ψ (χ) is a Pauli field that annihilates (creates) a heavy quark (antiquark), while G a µν = ∂ µ A a ν − ∂ ν A a µ is the noninteracting part of the field strength tensor G a µν , with A a µ being the gluon field. Notice that ψ |0 = 0 and 0| χ = 0. The Fourier decompositions of the free NRQCD fields are where p 0 = p 2 /(2m) in eqs. (C.2a) and (C.2a), and p 0 = |p| in eq. (C.2c). Moreover, s, λ, c and d denote the spin, polarization and color quantum numbers respectively. A quark field carries one spinor index i and one fundamental color index n, whereas a gluon field has an adjoint color index a and a Lorentz index µ attached to it. The quantities ξ n i (s, c) and η n i (s, c) should be understood as a product of a 2-spinor and a 3-color vector, i.e.
with a possible choice for the spinors (ξ, η) and 3-color vectors v 3 being Likewise, for the gluon field we have where ε µ (p, λ) is an ordinary polarization vector, while v a 8 (d) describes an 8-color vector with e.g. v 8 (1) = (1, 0, 0, 0, 0, 0, 0, 0) T . In order not to clutter the notation further, here we chose to suppress the flavor indices of the quark fields. Those can be made explicit exactly in the same manner as the color indices, i.e. by introducing corresponding 6-dimensional vectors.
Furthermore, in order to avoid dealing with unphysical degrees of freedom of massless vector bosons, we let the gluon field possess only transverse polarizations (radiation gauge). This is not an issue here, since we employ the operator formalism only as a convenient shortcut to derive Feynman rules. Repeating the same exercise using the BRST construction [146,147] would only complicate the derivation but yield exactly the same Feynman rules, so that we do not consider it useful here. We also omit the treatment of the gluon-ghost interactions, since the corresponding Feynman rules are exactly the same as in ordinary QCD.
The creation and annihilation operators of the NRQCD fields satisfy following nonvanishing (anti)commutation relations As it is customary in NRQCD, we define the 1-particle Fock states to have nonrelativistic normalization, so that By differentiating the exponentials in eqs. (C.9) one can also obtain the corresponding relations with temporal or spatial derivatives applied to the field operators e.g. with ∇ i = ∂ i = ∂/∂x i . It is also easy to derive auxiliary formulas containing products of fields, which may be convenient for pen and paper calculations e.g.

C.2 Feynman rules for pNRQCD
Finally, we analyze pNRQCD in the same manner as we did it for NRQCD. The main reason for doing so is to highlight some interesting aspects of the theory that make pNRQCD conceptually similar to ordinary nonrelativistic quantum mechanics. In weakly coupled pNRQCD our degrees of freedom are bilocal color singlet and color octet fields as well as multipole expanded gluons. As it has already been explained in section 2, the bilocal fields depend both on r and R, while gluons are sensitive only to R. The free part of the pNRQCD Lagrangian at O(1/m) and O(r 0 ) in the multipole expansion is given by a µνĜ µνa , (C. 18) with S ij and O ij defined as in eq. (2.4). From these definitions it follows, in particular, that where i, j and a are fundamental and adjoint color indices respectively. By Fourier expanding free singlet and octet fields in terms of their creation and annihilation operators we find where the color structure V a (c, c ) is defined as with v 3 (c) being the color vectors that have already been introduced in appendix C.1. As far as the kinematics is concerned, we have P ≡ (P 0 , P ) T and R ≡ (t, R) T  Regarding the gluon fields, the corresponding formulas given in appendix C.1 still apply, which is why we do not repeat them here. The only difference is that the exponential e ip·x should be replaced with e iP ·R . The derivation of the pNRQCD Feynman rules is now straightforward. For example, we can work out the Feynman rule for the singlet-octet chromoelectric dipole interaction with one gluon emission at O(1/m 0 ) and at O(r). From the term gV A (r) Tr(O † r · E S) in the Lagrangian we obtain which is also the Feynman rule for gV A (r) Tr(S † r · E O). As far as the octet-octet sector is concerned, the treatment of (gV B (r)/2) Tr(O † {r · E, O}) is equally simple and boils down to