Functional methods for Heavy Quark Effective Theory

We use functional methods to compute one-loop effects in Heavy Quark Effective Theory. The covariant derivative expansion technique facilitates the efficient extraction of matching coefficients and renormalization group evolution equations. This paper pro- vides the first demonstration that such calculations can be performed through the algebraic evaluation of the path integral for the class of effective field theories that are (i) constructed using a non-trivial one-to-many mode decomposition of the UV theory, and (ii) valid for non-relativistic kinematics. We discuss the interplay between operators that appear at intermediate steps and the constraints imposed by the residual Lorentz symmetry that is encoded as reparameterization invariance within the effective description. The tools presented here provide a systematic approach for computing corrections to higher order in the heavy mass expansion; precision applications include predictions for experimental data and connections to theoretical tests via lattice QCD. A set of pedagogical appendices comprehensively reviews modern approaches to performing functional calculations algebraically, and derives contributions from a term with open covariant derivatives for the first time.


Introduction
The modern perspective on Effective Field Theories (EFTs) takes them to be well-defined quantum field theories in their own right, with all the attendant structure that implies. In particular, one can define an EFT at the quantum level using a path integral. However, the standard approach to extracting observables is to sidestep this intimidating object by organizing perturbation theory through the use of Feynman diagrams. The recent reintroduction [1] of the covariant derivative expansion (CDE) [2][3][4] in a form tailored to tackling a broader set of EFTs has sparked a renaissance in the extraction of physics directly from the path integral using functional techniques. A noteworthy example of this progress has been a general matching formula at one-loop [5,6] that captures the far offshell fluctuations of heavy fields. This result can be used for the extraction of the effective operators and their Wilson coefficients that characterize the effects of the heavy physics, encoded in the so-called Universal One-Loop Effective Action (UOLEA) [1,[7][8][9][10]; one of its features is that it elegantly packages together what would be multiple independent Feynman diagram calculations. This has seen ready application to beyond the Standard Model (SM) scenarios by facilitating one-loop matching onto the so-called SMEFT. However, there exist a large number of EFTs whose relationship to UV physics cannot be captured by simply integrating out an entire heavy field. Here, we present the first application of the modern functional approach to such EFTs. Specifically, we work in a particular low-energy limit of the SM [11][12][13][14] expanded around a non-trivial background, the Heavy Quark Effective Theory (HQET) [13,[15][16][17], see also the reviews [18][19][20][21][22].

JHEP06(2020)164
There are many novel features of HQET when compared to the SMEFT. Perhaps the most striking is that -as a model of the long distance fluctuations of a heavy particle -HQET is a non-relativistic EFT. Obviously, the path integral is a valid description of a non-relativistic theory (after all, it was invented as an alternative description of quantum mechanics [23,24]), but the concrete demonstration that the functional approach is useful for precision non-relativistic field theory computations has been lacking until now. 1 Another intriguing aspect of HQET is that it invokes the concept of a mode expansion. A single heavy quark field is decomposed into a pair of fields, which model short and long distance fluctuations. Then, assuming a particular kinematic configuration, only the long distance modes can be accessed as external states. The short distance fluctuations can therefore be integrated out, which generates a tower of local EFT operators. It is this description that is called HQET. This type of one-to-many correspondence plays a role in many modern formulations of EFTs, and for the first time here we put such models on even firmer theoretical footing by computing observables directly from the path integral. Furthermore, a theoretically appealing aspect of this work is that the covariant derivative expansion of the functional integral manifests the symmetries of the theory in a transparent way. This is obviously true for gauge invariance, but we will also be able to approach the residual Lorentz symmetry known as Reparameterization Invariance (RPI) from a novel vantage point. In particular, we will identify an intermediate stage in our calculations where RPI becomes manifest. This provides a nice contrast to the Feynman diagram approach, where this invariance only holds once one sums the full set of relevant diagrams.
Through our application of the functional approach to HQET, we will expose some important features of the general formalism. Matching a UV theory onto an EFT, first done for HQET in ref. [29], can be performed diagrammatically by equating matrix elements computed with the two descriptions of the theory at a common kinematic point. This requires knowing the EFT operator expansion and identifying the operators that are relevant to the EFT matrix element calculation. By contract, using functional methods the generation of operators and matching of Wilson coefficients occurs in a single step. There is no need to specify the structure of the EFT before performing a matching calculation. Working through the example of HQET will show that the problem for how to marry the mode decomposition with functional methods has a nearly universal solution. Specifically, given an implementation of the mode decomposition using operator-valued projectors, one can derive a functional equation of motion for the short distance modes. This can be used to integrate out these modes, yielding a non-local EFT description that encodes the complete dynamics of the full theory in the relevant kinematic limit. This justifies the construction of the EFT as a full path integral over a well defined field.
Deriving concrete predictions using conventional techniques typically involves several steps and many subtleties, but a functional approach makes the procedure more algorithmic. The complicated multi-mode matching calculation has a simple form whose structure is elucidated by analyzing the resulting integrals using the method of regions [30,31]. As 1 To our knowledge, the only time the path integral has been discussed in the context of HQET was in [25] where the tree-level formulation of the theory was first derived, and in [26][27][28] where RPI was briefly discussed in the context of the path integral.

JHEP06(2020)164
is well known, matching in HQET only receives support from diagrams that have loops with both a short distance mode and a mode that propagates in the EFT. Understanding how to access this exact class of diagrams using functional methods has been the subject of some confusion. Showing that we can derive matching and running in HQET with these methods is a conclusive demonstration that functional methods provide a complete framework at one-loop. These results will be summarized below as a simple master formula that encodes the matching of QCD onto HQET at one-loop and to any order in the heavy mass expansion, see eq. (5.4).
There are important practical implications of this work. One of the primary purposes of HQET is to provide efficient methods for precision calculations within the SM. There are a number of EFTs widely used to facilitation precise predictions including Soft Collinear Effective Theory [32][33][34], theories of non-relativistic bound states such as nrQCD [35][36][37], and others. Improving on the precision of a calculation often requires the application of novel theoretical approaches, with a goal to provide computational benefits over a naive perturbative expansion evaluated using Feynman diagrams. Specifically, functional technology has seen little use in this context, despite the fact that some of the simplifications it provides are arguably most relevant to the questions these kinematic EFTs are designed to answer. By showing we can reproduce non-trivial matching and running results for HQET here, we open the door to understanding how to apply functional techniques in these other contexts. We additionally lay the foundation for performing new calculations within HQET itself. In particular, one can now perform matching calculations to higher order in the heavy mass expansion, which would be relevant to high precision measurements made at experiments such as LHCb and Belle II, and for connecting lattice gauge theory calculations performed in the heavy quark limit to their continuum limits [38,39].
The rest of the paper is organized as follows. In the rest of section 1, we provide a summary of the known results that we reproduce in a novel way throughout this paper. Next, section 2 provides a condensed introduction to HQET. Section 3 provides a summary of the traditional method of calculating the simplest piece of the one-loop HQET matching, the residue matching. (Readers familiar with HQET could skip these two sections.) In section 4, we review how to extract matching and running from a functional determinant. (Readers familiar with functional methods could skip this section.) In section 5, we introduce the use of functional integration to construct kinematic EFTs, and clarify how the method of regions simplifies the derivation of the resulting master matching formula. Section 6 is then dedicated to an explanation of how such matching is performed to one-loop order. This is followed by section 7, which strengthens the case for functional methods by providing additional matching calculations, along with an example that shows how operator running can be derived in the formalism as well. Section 8 then concludes.
An extensive set of pedagogical appendices provide an introduction to many of the relevant technical details. The "covariant derivative expansion" technique used here was originally invented in 1980s [2][3][4], we refer to this as "original CDE" in appendix B. This was reintroduced in the context of modern EFT calculations in ref. [1], and has been applied by refs. [7][8][9][10]40] to develop one-loop universal effective actions. A closely related variant of the original CDE, which we call "simplified CDE" in appendix B, was proposed in ref. [5].

JHEP06(2020)164
This more rudimentary version of the CDE turns out to be significantly more convenient for extracting operators that do not involve a gauge field strength. Most of the results in the main text of this paper were derived using this simplified CDE. In appendix B, we clarify the relations between these two versions of the CDE. For completeness, in appendix B we also provide some simple universal results (tabulated in appendix B.4.2) for functional traces derived using the CDE. These include the famous elliptic operator, see eq. (B.84a), which is the central object of study in the development of the UOLEA. Additionally, we provide the first computation of the contributions to the UOLEA from a term with an open covariant derivative (truncated at dimension four). This result is provided in eq. (B.84b), and was previously unknown as emphasized by refs. [10,41,42]. A variety of RGE example calculations are given in appendix C, and a generalization of the Heavy-Heavy matching calculation is provided in appendix D.

Summary of results
In what follows, we will present our methods by way of a few canonical matching and running calculations. First, we derive the high scale HQET Wilson coefficients by matching QCD onto the EFT for the heavy-light currentsq γ µ (γ 5 ) Q and the heavy-heavy currents Q 1 γ µ (γ 5 ) Q 2 using purely functional methods. We also present a functional derivation of running effects, using the first subleading operators in the HQET Lagrangian as a concrete example.
In order to fix our notation and make comparison to standard results straightforward, we provide a brief compendium of the results we will reproduce in this paper. The conventional derivation is presented in much more detail in standard references, e.g., refs. [19,22]. 2 This list provides a useful context for our goals here.
Propagator residue. When performing a matching calculation in any off-shell scheme, such as MS, it is critical to track the difference in propagator residues (necessary for obtaining the desired S-matrix elements using the LSZ reduction procedure) when moving from the full theory to the EFT. Taking the functional point of view, the resulting effect shows up as a rescaling of the kinetic terms for the EFT fields. This can then be moved to its canonical position in the Wilson coefficients by a field redefinition. Thus, the first step performed in what follows is to derive the one-loop corrections to the kinetic terms from QCD: where h v models the long distance fluctuations of the heavy quark field, v µ is the reference vector defined in eq. (2.2) below, D µ is the covariant derivative, α s is the strong fine 2 Notation. We have mostly chosen to follow the notation of ref. [22], but have made a number of minor changes, which is partially why we provide this summary here. For dimensional regularization (dim. reg.), our convention is to work in d = 4 − 2 dimensions. We have also chosen to use a different standard notation of the EFT heavy quark fields. Additionally, we have reserved the superscript numeral in parenthesis notation to denote loop order, R (0) is tree-level, R (1) is the one-loop correction, and so on. This is different from the standard notation in the HQET community, e.g. in ref. [22] the one-loop correction to the residue is denoted as R1. Additionally, we note that we will only denote the loop order of the renormalized terms, i.e., counterterms are implicit.

JHEP06(2020)164
structure constant, and ∆R (1) is the matching correction for the propagator residue. In a traditional matching calculation, the correction comes from the difference of the residues between the full and EFT descriptions: While this result is sensitive to the choice of matching scale µ due to presence of a UV divergence; the cancellation of terms that depend on the IR regulator serves as a check that the IR behavior of the two theories is identical.
Heavy-light current. At leading order in the heavy-mass expansion, two HQET operators have the same quantum numbers as the heavy-light vector current of QCD, implying that they can appear in the matching: where C V,i is the to-be-calculated matching coefficient for the vector current, which is a function of the heavy quark mass m Q and the strong coupling; there is an analogous expression for the axial current derived by replacing The answer up to one-loop (see eq. (3.48) in ref. [22]) can be written as HL,2 α s + · · · , (1.4b) where ∆R (1) is defined in eq. (1.2), and the other terms (see eqs. (3.66) and (3.73) in ref. [22], respectively) are where 1 − γ E + ln 4π is the standard factor that is subtracted when using the MS scheme, with γ E denoting the Euler-Mascheroni constant. Then the axial matching coefficients are given by C A,1 = C V,1 and C A,2 = −C V,2 . The one-loop matching coefficients are thus (compare with eq. (3.74) in ref. [22])

JHEP06(2020)164
Heavy-heavy current. In case of the heavy-heavy currents, for simplicity we will take the special kinematic choice of zero recoil, corresponding to v 1 = v 2 in the EFT (the generalization to v 1 = v 2 , first done in ref. [43], is discussed in appendix D). In this limit, all possible HQET operators at leading order in the mass expansion are equal by the equations of motion, and the matching between QCD and HQET is simply 3 with an analogous expression for the axial current given by the replacement η V → η A and γ µ → γ µ γ 5 . The results can be parameterized as (compare with eqs. (3.98) and (3.101) in ref. [22]) where m 1,2 corresponds to the mass of Q 1,2 . The matching coefficients then take the form in agreement with eqs. (3.97), (3.99), and (3.101) of ref. [22].
β-functions. Finally, we reproduce the expressions for the running of the HQET matching coefficients at one-loop and at O(1/m Q ) (compare with eq. (4.8) in ref. [22]): The Renormalization Group Equations (RGEs) are [44][45][46] µ d dµ where C A denotes the Casimir factor for the adjoint representation In what follows, we will show how to reproduce all of these results using functional methods equipped with the CDE technique.

JHEP06(2020)164 2 Heavy Quark Effective Theory
One of our main goals here is to initiate the study of functional methods for higher-order calculations in EFTs that are derived by performing multi-modal decomposition of full theory fields. Such EFTs are quite common; they occur when the kinematics of the process being studied selects a preferred reference frame due to, e.g. a conservation law preventing the decay of a heavy particle, or a measurement function that forces the external states into a particular region of phase space. These restrictions imply that there are full theory modes which cannot be put on-shell within the EFT regime of validity, and so it is sensible to integrate them out. This procedure breaks the full theory space-time symmetries to some subgroup, while potentially also introducing new internal ones. A theory of this type is the Heavy Quark Effective Theory (HQET), which describes the fluctuations of a heavy quark (m Q Λ QCD ) in the presence of light QCD charged degrees of freedom. The simplifications and universal behavior of QCD in the heavy-mass limit were first appreciated by refs. [11,12] and especially refs. [13,14]. A non-covariant EFT making this behavior manifest was later developed and shown to be well-behaved in perturbation theory [15,17], and finally given a covariant formulation [16]. We provide a review of HQET here, with a particular emphasis on the equation of motion due to the critical role it plays in what follows. The reader familiar with HQET can skip ahead to section 4, while more details can be found in, e.g. section 4.1 of ref. [22].
The full theory (QCD) Lagrangian for a heavy quark includes where Q is our heavy quark, and the covariant derivative only includes the QCD interactions, D µ = ∂ µ − ig s G a µ T a . In the following, additional interactions of the heavy quark will be modeled through the introduction of current operators as needed.
Naively, it does not seem that the Lagrangian of eq. (2.1) has a good expansion about the m Q → ∞ limit. The key insight that allows one to circumventing this issue lies in an appropriately chosen phase redefinition of the field, whose purpose is to cancel the mass term for certain components. The full quark field can then be separated into so-called short distance and long distance fields, where the latter become approximately mass-independent. Physically, this is motivated by the realization that the heavy quark cannot be pushed very far off-shell by degrees of freedom for which |q| Λ QCD . To make this manifest, we decompose a heavy quark's momentum as where v µ is a unit time-like vector, and k µ is the residual heavy quark momentum which models small fluctuations about its mass shell. For kinematic configurations such that all Lorentz invariants depending on k µ are small compared to m Q , a truncation at finite order in the |k|/m Q expansion is justified. 4 Since the theory is expanded around the m Q → ∞ JHEP06(2020)164 limit, the structure of HQET does not know about the mass of the heavy quark except through various non-dynamical quantities. In particular, no sensitivity to m Q appears in any of the calculations beyond what is encoded in the structure of the matching coefficients. As a brief aside, we note that the decomposition in eq. (2.2) is not unique. In particular, a simultaneous transformation of k µ and v µ by a fixed vector: leaves p µ unchanged. Enforcing that v µ remains a unit vector implies that δk µ must satisfy v · δk = δk 2 /(2m Q ). Reparameterization invariance (RPI) is then the statement that physical observables cannot depend on δk µ , thereby enforcing the residual Lorentz invariance of the underlying full theory. We will explore the interplay between RPI and the functional approach in section 6.1 below. In order to make use of eq. (2.2), we decompose the heavy quark field into two fields by first extracting the rapidly-varying phase that remains unchanged by low-energy interactions, and then using v µ -dependent projectors to split the remaining field into short distance and long distance components: or equivalently Plugging eq. (2.5) into eq. (2.1) yields The Lagrangian in eq. (2.6) makes the interpretation of H v as the heavy mode manifestthis field has an effective mass of 2m Q , permitting a description at lower energies in terms of h v alone. To formally integrate out H v at tree-level, we solve for its equation of motion and plug it into eq. (2.6), which yields

JHEP06(2020)164
Provided that the momentum of the field h v satisfies |k| m Q , eq. (2.9) can be expressed into a convergent series of local terms (2.10) Since the leading term is then insensitive to the Dirac structure of the spinor components of the field, the theory has gained an approximate SU(2) heavy-spin symmetry, which can be expanded to SU(2n h ) in the presence of n h heavy flavors. This is an example of an aforementioned emergent symmetry of the EFT. To understand how this procedure for deriving the EFT Lagrangian can be interpreted at the path integral level, we note that the projection operators in eq. (2.4) imply that h v and H v can be treated as two orthogonal projections of Q. As such, the path-integral measure factorizes: Since the resulting Lagrangian is quadratic in H v , the Gaussian integral over the shortdistance field immediately yields eq. (2.9). Therefore, the procedure described here allows us to literally integrate out the short distance modes of our original quark. This elegant derivation of the tree-level HQET Lagrangian by way of the path-integral measure was first presented in ref. [25]. Furthermore, it also justifies the use of the functional approach to matching and running developed in ref. [5]. Note that there would be no significant obstructions even if our projectors were quantum operators instead of simply being functions of the kinematics as above, since all objects in the path integral are treated as operators acting on fields. As long as the mode decomposition that is used to define an EFT can be written in terms of operator-valued projectors, we expected the methods developed here to be universally applicable.

Decoupling the heavy quark
The decomposition of the full QCD Lagrangian given in eq. (2.6) makes it clear that H v can be identified as a short distance mode which one can integrate out in the limit that |k| m Q for all fields in a given process. This procedure yields the non-local Lagrangian given in eq. (2.9), which makes predictions that are equivalent to QCD for processes involving only h v modes in the external states. Expanding the non-local Lagrangian using eq. (2.10) and truncating the series at some order in 1/m Q yields an EFT which is valid for momenta |k| m Q . In this subsection, we will briefly review the connection between this procedure and the fact that the heavy quark should not contribute to the running of the QCD gauge coupling at scales below m Q . This both has a conceptual benefit, and will also be of practical importance since similar arguments will be used below when we derive our master formula for one-loop matching given in eq. (5.4).
For this argument, we will work directly with the h v and H v fields, and we will use the terms that are diagonal in these fields to derive propagators, while the mixed terms will be treated as interactions. Two features are of critical importance, the fact that the kinetic terms are linear in the momentum of the state, and the relative minus sign between the

JHEP06(2020)164
h v and the H v kinetic terms. The later fact implies that the i factors in the propagators will take the opposite sign such that the same Wick rotation can be used to Euclideanize diagrams involving both h v and H v propagators.
The QCD Lagrangian written as eq. (2.6) has three types of couplings between the heavy quark modes and the gluon: The integral that results when computing the contribution to the vacuum polarization for the gluon from the two diagonal loops schematically take the form where m is an IR regulator for the h v loop and is equal to 2m Q for the H v loop. Due to the linear nature of the kinetic terms and the sign on the factor of i , when integrating over q 0 , the poles reside on only one side of the real axis, and one can deform the contour away from all them yielding zero contribution. However, the situation is different for the mixed loop, where the integral takes the form Now we see that the opposite sign on the i terms implies that there is a pole in both the positive and negative Im q 0 half-planes. The contour will enclose a pole for any possible deformation, yielding a non-zero contribution to the β-function.
The argument above makes it clear why in HQET the heavy quark does not contribute to the RGEs of the gauge coupling. Once we construct HQET by integrating out H v and expanding, the heavy field H v is non-propagating. Therefore, there are no diagrams that yield contributions of the type in eq. (2.13). The mode h v only yields potentially relevant integrals of the type in eq. (2.12), which vanish as we have argued. Similar reasoning will be critical to the derivation of the master formula for matching in eq. (5.4) below.

Residues using Feynman diagrams
In this section, we will review the diagrammatic approach to calculating matching coefficients, using the simple example of the residue of the on-shell propagator for concreteness. Along the way, we will encounter an order of limits issue, which is a manifestation of the IR divergence structure of QCD. We will then revisit the calculation by relying on the so-called method of regions [30,31] that will avoid the need to deal with this subtlety directly. This has the additional benefit of providing a familiar setting to review the method of regions, which will be a critical tool in the derivation of our master formula for HQET matching coefficients below.
When performing matching calculations, one typically equates matrix elements as calculated in the full theory and the EFT. Particular care must be taken to include the appropriate residue factors for the external states to ensure that the LSZ reduction is correctly implemented. One way to extract the residue for the heavy quark Q is to take the JHEP06(2020)164 derivative of the 1PI corrections to the propagator −iΣ( / p), and evaluate it on the mass shell. We will compute this factor for a quark in QCD R Q = 1 + R (1) Q α s + · · · , and for a quark in HQET R h = 1 + R (1) h α s + · · · , from which we get the quantity that appears in matching calculations ∆R (1) h , see eq. (1.2).

Residue in QCD
The one-loop QCD residue R (1) Q can be obtained by computing the two-point 1PI function Here as well as throughout this paper, "c.t." denotes the counter term contributions, and we take the Feynman gauge ξ = 1 for gauge boson propagators. Performing standard manipulations, we derive 5 where the MS counter terms A ct and B ct are derived by taking the expansion of A p 2 and B p 2 . This yields These results are finite but non-analytic at p 2 = m 2 Q . In particular, their derivatives with respect to p 2 are divergent when evaluated at p 2 = m 2 Q . These are a manifestation of IR divergences that appear when taking on-shell kinematics.

JHEP06(2020)164
One way to side step this issue, thereby allowing us to extract the residue, is to keep = 0 until after taking the derivative. It is then straightforward to derive where is specifically regulating the IR divergence. Next, we will derive the residue in HQET, where we will see that the same IR divergent terms appear. Therefore, the object of interest ∆R (1) is IR finite. This is to be expected, since one of the standard tests that one has correctly implemented the matching procedure, namely that one is working with the correct low energy description, is to check that the IR of the full theory and EFT have the same divergence structure. We will provide a procedure for directly extracting ∆R (1) that avoids this IR subtlety utilizing the method of regions, see section 3.3 below.

Residue in HQET
Next, we perform the two-point function calculation in HQET. Diagrammatically, the structure is identical to the QCD calculation where the relativistic quark propagator is replaced by the HQET propagator, yielding which evaluates to 6 where the counter term is again determined in the MS scheme. 7 Noting that to zeroth order in 1/m Q the on-shell condition for h v is v · k = 0, we evaluate 6 See eqs. (3.67) and (3.69) in ref. [22]. 7 Note that when computing the counter term in dim. reg., one must be careful to isolate the UV divergence. This is done here by keeping v · k = 0 as an IR regulator at intermediate steps. This is why we must take a derivative of eq. (3.6b) before sending v · k → 0 to derive eq. (3.7b). If instead we took v · k → 0 first, we would effectively be using dim. reg. to regulate the IR divergence as well.

JHEP06(2020)164
which yields Similar to above, the lim →0 Σ HQET is not analytic at v · k = 0 due to an IR divergence, and so we had to defer taking the → 0 limit until after taking the derivative with respect to v · k.

Residue difference from the method of regions
The IR divergences in eqs. (3.4) and (3.8) are the same, as they must be if the EFT correctly captures the dynamics of the full theory below a certain scale. This implies that the residue difference is IR finite: Operationally, we were forced to maintain = 0 until we evaluated this difference. Therefore, it would be convenient to have an approach that would allow us to compute ∆R (1) directly. This can be accomplished by exploiting a technique known as the method of regions [30,31].
To begin, we rewrite some expressions in order to make the comparison between R h more obvious. Within QCD, we set p µ = m Q v µ + k µ and define (3.11) Note that we have changed the evaluation condition from v · k = 0 to / p = m Q in the first line, which is valid up to O(1/m Q ). Similarly, we define the HQET quantity This allows us to simply express the difference as (3.14) In order to evaluate this difference, we take the integral expressions for Σ QCD and Σ HQET given in eqs. (3.1) and (3.5) Note that the integrands of Ξ QCD and Ξ HQET are equal in the heavy quark limit: Naively, one might be tempted to conclude that Ξ QCD = Ξ HQET at leading order in 1/m Q , which would imply that ∆R (1) vanishes. This would be in conflict with the result derived in eq. (3.9). This is due to an order of limits issue: since the integral is taken over all momenta, one cannot expanding the integrand for large m Q before integrating.
The method of regions [30,31] is a technique for consistently expanding the integrands such that non-analytic dependence on small parameters is correctly reproduced after integration. The key is to isolate regions of the integration domain that are dominated by single scale contributions. For example, the QCD integral contains two physical scales that are separated by a large hierarchy |k| m Q . To expand the integrand, we introduce an intermediate cutoff scale Λ such that |k| Λ m Q , which can be used to split the integral over q µ into two pieces. The first receives support from the soft region |q| < Λ, while the second is non-zero due to the hard region |q| > Λ: The domain of integration is now bounded, so that one can expand the integrand according to the assumed scaling in each region, while maintaining |k| m Q fixed. The soft region can be isolated by assuming |q| < Λ m Q holds, while the hard region where |k| Λ < |q| yields a different expansion. Once the integrand has been expanded, the domain of integration can be restored to infinity when using dimensional regularization so that the explicit cutoff Λ no longer appears -the contribution from the extended integration limits is scaleless and therefore vanishes. This is reflected by our notation since we drop the Λ dependence for the expressions where the domain of the integral is taken to infinity.
Applying this procedure to the two-point function integrals, we recognize that the soft expansion is identical to the naive approach in eq. (3.16). Therefore, we conclude that This tells us that the residue difference is fully determined by the hard region from QCD:

JHEP06(2020)164
Note that eq. (3.18) holds to all order in 1/m Q . In eq. (3. 16), it was demonstrated only at the leading order, i.e. m Q → ∞. To explicitly verify this at higher orders in 1/m Q , one needs to include additional diagrams that contribute to eq. (3.15b) due to higher order operators in HQET. Then the result in eq. (3.19) still holds, with an appropriate modification of the evaluation condition v · k = 0. The practical implication of eq. (3.19) is that we no longer have to track the individual IR divergences. Instead, the residue difference ∆R (1) is determined by a single loop integral Ξ QCD,hard , and the only divergences that appear are from the UV, which can be subtracted using counter terms.
For completeness, we evaluate this hard region integral. Starting with eq. (3.15a), we expand the integrand assuming |k| |q| and |k| m Q : where we have truncated the expansion up to the linear order in k, since this is all that is required to isolate the residue. Evaluating the integral yields which reproduces the result in eq. (3.9).

Matching and running using functional methods
This section reviews how to utilize the functional approach for calculating matching coefficients and RGEs. To set the stage, we first briefly review the one-particle irreducible (1PI) effective action, which is a key object in functional methods that one could compute directly by evaluating the path integral. The general matching condition can be compactly expressed by equating the 1PI effective actions of the UV theory and the EFT, whose solution gives us a direct link between the Lagrangians of the two theories, e.g. see ref. [5]. Finally, we also provide a brief review on how to extract RGEs from the 1PI effective action; more details for RGE calculations are provided in appendix C.

1PI effective action from a functional determinant
In modern functional methods, the central object of study is the so-called 1PI effective action Γ[φ], where φ collectively denotes all the fields in the theory. It is a generating JHEP06(2020)164 functional for the 1PI correlation functions: One can in principle extract any perturbative quantum field theoretic prediction from Γ[φ]. Concretely, the 1PI effective action up to one-loop order can be organized as a loop expansion: where each piece can be extracted from the Lagrangian as , and "Sdet" is the so-called super-determinant, which tracks the minus sign difference between fermionic and bosonic loops. The normalization factor i 2 assumes that we are tracing over "real" degrees of freedom: a complex scalar field should be separated into its real and imaginary parts, and Dirac fermions should be decomposed as discussed in section 4 of ref. [5].

Matching from a functional determinant
In general, one performs a matching calculation by equating the EFT L EFT (φ) with a UV theory L UV (φ, Φ) at a matching scale, where φ (Φ) denotes the light (heavy) particles. The general matching condition is that all the one-light-particle irreducible (1LPI) diagrams agree at the matching scale [47,48]. The natural definition of this statement when using functional methods is to enforce that the so-called 1LPI effective action Γ L [φ] -the generating functional for all the 1LPI correlation functions -of each theory coincides at the matching scale: One can relate the 1LPI to the 1PI effective action for the UV theory and the EFT: The relation for the EFT is trivial. By contrast, to derive Γ L, UV , one must integrate out the heavy field by plugging in the solution Φ = Φ c [φ] to its equation of motion.
Solving eq. (4.4) in general is nontrivial. However, ref. [5] was able to make significant progress towards this goal by deriving the following general expression for the EFT

JHEP06(2020)164
Lagrangian up to one-loop order: 8 , (4.7a) (4.7b) A few clarifications are in order: • In this approach, we do need to know the EFT operators O i (φ) in advance. We simply obtain them (with the appropriate one-loop coefficient) by evaluating the right-hand sides of eqs. (4.6) and (4.7).
• The one-loop contributions to the Wilson coefficients derive from two classes of loops in the UV theory. The "heavy" Wilson coefficients C i,heavy collect contributions from loops where only Φ appear, while the "mixed" Wilson coefficients C (1) i,mixed are generated by loops with both Φ and φ.
• The non-local Lagrangian appearing in eq. (4.7b) is given by where the heavy field is integrated out using the tree-level solution to its equations of motion Φ c [φ]. No expansion in the heavy masses is performed at this stage, and hence the resulting Lagrangian contains non-local terms, e.g. see eq. (2.9). This is in contrast to the tree-level EFT action S , which is derived by expanding the non-local Lagrangian in the heavy mass limit; therefore, it only contains local terms, e.g. see eq. (2.10). As we have discussed extensively in section 3.3, expanding the action in the heavy particle limit yields a critical difference between the two descriptions, where the non-trivial effects result from being careful about the order of limits -performing the heavy mass expansion does not commute with taking the functional determinant (which is essentially equivalent to performing the loop integration). Although, the two terms in the first line of eq. (4.7b) are not equal, JHEP06(2020)164 they are intimately related in such a way that allows for the following simplification: using method of regions, 9 we can identify the second term as equivalent to the soft region of the first term. Their difference leaves behind the hard region alone, as shown in the second line of eq. (4.7b).

RGEs from the 1PI effective action
In this section, we briefly review the general procedure for using the 1PI effective action to compute the RGEs [5]. Given a Lagrangian we are interested in deriving the RGEs for a coupling λ, where O λ is the corresponding operator, O K denotes the kinetic terms, and the "· · · " allow for the possibility of additional interactions. The first step is to compute the 1PI effective action, which takes the form where a K and a λ encode the one-loop corrections. Next, we renormalize the kinetic terms to their canonical value by rescaling the fields φ: where b λ is derived by Taylor expanding the version of Γ[φ] with canonical kinetic terms as a series in λ. Finally, the RGE for λ is obtained: For reference, a number of standard RGEs are derived using this method in appendix C; additionally, functional methods were used to compute the bosonic dimension-6 SMEFT RGEs in ref. [49]. In section 7, we will apply this formalism in HQET to show that the kinetic term does not run, and will derive the RGEs for the Wilson coefficient of the HQET magnetic dipole moment operator.

HQET matching from a functional determinant
In this section, we apply the general matching result eqs. (4.6) and (4.7) to HQET. As we will see, salient features of HQET will result in further simplifications for matching calculations, which will be summarized by our master formula eq. (5.4) for HQET one-loop matching coefficients. In this case, we identify Φ = H v , H v as the short distance modes of the heavy quark, and φ collectively as all the propagating long distance degrees of freedom JHEP06(2020)164 modeled by HQET. There are two special features of matching QCD onto HQET that simplify the resulting master formula. First, taking second functional derivatives of eq. (2.6) with respect to H v andH v , we see that the matching coefficient C (1) i,heavy vanishes at one-loop: This functional determinant is zero due to the same contour arguments presented in section 2.1 (see also [25]), where we discussed decoupling the heavy quark. 10 This indicates that all the one-loop matching contributions for HQET are mixed loop contributions of the type given in eq. (4.7b). Therefore, the non-local HQET Lagrangian which is given explicitly in eq. (2.9), is equivalent to QCD with respect to the dynamics of the light fields φ.
Another useful feature of HQET is that loop corrections are scaleless [50], and therefore they vanish when using dim. reg. . This means that the second term in the first line of eq. (4.7b) is also zero: We emphasize that this expression is valid when (i) we are computing on-shell S-matrix elements, and (ii) when dim. reg. is used to regularize both UV and IR divergences (see section 3.2 for a calculation where we regulated the IR by keeping the long distance fluctuations off shell). Taking both simplifications into account, we arrive at our master formula for computing the Wilson coefficients of HQET 11 In particular, note that we have restored the vanishing region, implying that when using this result we do not need to perform any method of regions style expansion of the integrals that result from taking the functional trace, and can simple evaluate the integrals that result from this procedure directly. In sections 6 and 7, we will demonstrate how to apply this formula by working out a number of explicit examples. We will additionally see that this formula makes symmetry properties such as gauge invariance and RPI more manifest. 10 There is an even simpler (albeit gauge dependent) argument also given in ref. [25]: if one takes the v · A = 0 gauge, the determinant no longer depends on any field, and its evaluation simply yields a constant that is absorbed into the path integral measure. 11 There is one caveat to keep track of when applying this formula, which is that it is valid when one sets all the light masses that appear in a loop integral to zero. If one is interested in computing power corrections proportional to a light mass, then eq. (4.7b) should be used: one should keep the light masses non-zero, and the hard region must be isolated before integrating.

JHEP06(2020)164 6 Residue difference from a functional determinant
This section provides a first non-trivial example of the procedure for matching between QCD and HQET at one-loop using functional methods: the calculation of the difference between the propagator residue for QCD and HQET. Along the way, we will highlight many of the simplifying benefits of performing these calculations directly from the path integral. Then before moving on to a number of additional (more complicated) examples in the next section, we briefly discuss how RPI manifests in the functional language. Before getting into the details, we provide a simple road map highlighting the critical steps of the calculation: (i) Beginning with the UV Lagrangian, we will integrate out the short distance mode using the tree-level equations of motion. This will result in a non-local Lagrangian (using background fields as necessary to implement gauge fixing) that encodes a complete description of the dynamics of the light modes. See eq. (6.1).
(ii) Given this non-local Lagrangian, we will then derive the matrix of the second-order functional derivatives (see eq. (6.5)), whose determinant is the central object to evaluate according to our master formula eq. (5.4).
(iii) Next, we will rewrite the determinant into an efficient form for explicit evaluation via row reduction. See eq. (6.8).
(iv) Finally, we will evaluate this trace using the methodology of the covariant derivative expansion as reviewed in appendix B. From this expansion, we will isolate the operators of interest. This yields the second line of eq. (6.12), which is an expression for one-loop Wilson coefficients multiplied by the appropriate operators. The last step will be to explicitly evaluate the loop integral. See the final line of eq. (6.12).
Note that along the way, we will drop any terms which do not contribute to the operator of interest for simplicity. The rest of this section is devoted to the explicit calculation of the residue difference. Our goal is to apply eq. (5.4) in order to derive the one-loop correction to the two-point function of heavy quarks. This implies that we only need to take variations of the Lagrangian with respect to the gluon and the heavy quark field. The first step is to derive the non-local HQET action. Following section 2, where we reviewed the derivation of the HQET Lagrangian, we start with QCD (including the gluon kinetic term, gauge fixing, and ghost contributions) and integrate out the short distance quark modes: where the explicit gauge-fixing term L gf and ghost term L gh are specified in eqs. (B.74) and (B.75). Next, we take functional variations of S non-local HQET . We will treat gauge bosons using the background field prescription described in appendix B.3: where G B,µ is a background field and A µ encodes the fluctuations. Covariant derivatives are expanded similarly: where D B,µ is the covariant derivative that includes the background gauge field. After taking functional variations with respect to A µ , we replace We only need to take the second variation of the Lagrangian in eq. (6.1) with respect to A a µ , h v , andh v to compute the residue difference: where terms not relevant for the residue operator in U µν,ab , Γ µ,a , andΓ µ,a are discarded: The reality of the Lagrangian implies that B † = γ 0 B γ 0 andΓ µ,a = (Γ µ,a ) † γ 0 must hold, and it is straightforward to verify these with the explicit expressions in eq. (6.6).
Next, we evaluate the functional determinant by row reducing eq. (6.5): Following the prescription in section 4.2, we extract the residue operator part of the one-loop HQET Lagrangian using where the subscripts denote the field that is being traced over, i.e., the state that is running in the loop, and starting with the third line, we have truncated the series and dropped terms that do not contribute to the residue operator. Note that in the last line, the trace over the gluon space is taken (i.e. µν and ab indices are contracted). Next, we evaluate this expression using the explicit objects given in eq. (6.6), and simplify the two traces as follows: where we used T a h T a h = C F = 4/3, and In performing these manipulations, we set [D µ , D ν ] = 0 since this commutator returns the field strength G a µν , which does not contribute to the residue difference. Putting these two results together, we obtain (6.11)

JHEP06(2020)164
Finally, we evaluate this functional trace using the simplified CDE technique described in appendix B.2.2 to derive 12 The first line follows from eq. (6.11) by taking the steps shown in eq. (B.46). In the second line, we have expanded in the covariant derivative D µ -this is the explicit step that uses the Covariant Derivative Expansion. Note that in taking this expansion we need to invoke D m Q to justify only keeping the residue operator. Its Wilson coefficient is given by a loop integral with a nonzero hard region but a vanishing soft region. This verifies our general statement made in section 5, in particular eq. (5.3). In the third line, we have evaluated the loop integral, with the MS scheme counter terms. The result agrees with eq. (1.2), providing our first demonstration for using functional methods to compute HQET loop effects. Before moving on to additional examples in section 7, we will briefly discuss RPI in the context of the residue difference calculation.

When does RPI become manifest?
In this section, we briefly explore the interplay between RPI and the calculation detailed in the previous section. 13 Our goal here is to understand at what point RPI holds when computing with the functional approach. For simplicity, we will only explore this question to leading order in 1/m Q in the RPI transformations. The full RPI transformations are significantly more complicated, since the eigenstates h v and H v rotate into each other at subleading order, see [26-28, 52, 53] for a discussion.

JHEP06(2020)164
The RPI transformation shifts 14 We will now show that this is a good symmetry of HQET by noting that the non-local HQET Lagrangian is built using the operator B defined in eq. (6.6). In particular, this object transforms covariantly as and hence . Alternatively, one can check the interplay between the two terms in the first line of eq. (6.14). We begin by analyzing the explicit transformation of the local term. Keeping all terms to O(1/m Q ), we find This change is compensated by a corresponding shift in the non-local piecē where we used the identities Then clearly eq. (6.15) holds.
14 To reduce notational clutter, we have opted to use k as opposed to δk here to track the change induced by RPI when comparing with eq. (2.3). Given the form of the relevant expressions, there is no ambiguity.

JHEP06(2020)164
The loop-level Lagrangian is given by the super-determinant of the second variation of eq. (6.1). Therefore, it is also invariant under RPI. Let us explicitly check this to the order O(1/m Q ) for the residue difference calculation presented in the previous section. The various components defined in eq. (6.6) shift under RPI as From here, we can check the transformations of each term that appears in the argument of the functional trace in eq. (6.8): where in the second derivation we used This demonstrates an elegant feature of functional methods when applied to HQET. Specifically, in the last line of eq. (6.8) the terms within the square brackets are manifestly invariant under RPI before evaluating the trace. This is in contrast with the Feynman diagram approach, where one must sum the full set of Feynman diagrams before the RPI symmetry becomes apparent.

JHEP06(2020)164
7 More matching and running calculations We have provided an explicit demonstration of a functional calculation in HQET for the simplest non-trivial example above. However, we have not yet made contact with an actual observable quantity. The purpose of this section is to do so by exploring more examples of matching calculations and a derivation of the RGEs for some Wilson coefficients. This will provide overwhelming evidence that these techniques capture all of relevant physics. Using the formalism developed in refs. [29,54], once the one-loop matching of operators is known, one-loop relations between exclusive quantities, such as decay constants and form factors, are straightforward to extract. Furthermore, explaining these examples will provide us with the opportunity to highlight some additional subtle aspects of applying functional methods.

Heavy-light current matching
This section is devoted to our next example, matching the heavy-light current. The heavylight operator in QCD is defined as where q is a light quark and Q is a heavy quark that will be treated as an HQET field. Since one should use the reference vector v µ when constructing HQET operators, two operators can be written down at leading order in O(1/m Q ) that manifest the same symmetry properties as the heavy-light QCD operator: The matching condition requires that matrix elements derived using QCD and HQET at a convenient matching scale µ are equal, 15 where the Wilson coefficients for the HQET operators C 1,2 = C 1,2 (m Q /µ, α s (µ)) are functions of µ and α s . At one-loop order, the two matrix elements can be expressed schematically: HL,1 α s γ µ + V HL,2 α s v µ u p, s , h α s = R h − 1, and R q denotes the residue for the light quark propagator, which is the same in QCD and HQET. Equating the two expressions in eq. (7.4) lead to the simplified form given in eq. (1.4). In particular, R q drops out and the matching coefficient only depends on the residue difference ∆R (1) h , computed in the previous section.
In order to extract the heavy-light matching coefficient from the master matching formula given in eq. (5.4), we follow the same steps outlined at the beginning of section 6. The set of fluctuating fields that contribute are A a µ , h v ,h v , q, andq. We need the equation of motion for the H v field, which can be derived from the UV Lagrangian (where sources J ± µ for the heavy-light current are now included) Note that the sources J ± µ are not dynamical fields, but must be included to ensure that the desired operators appear when matching the full theory and EFT actions. We split the heavy quark field as in eq. (2.5) above, and derive the equation of motion for the short distance mode: Since J + µ is merely a source, we are free to redefine it to absorb the phase, J + µ → J + µ e −i m Q v·x . Plugging the equation of motion eq. (7.6) into the UV Lagrangian eq. (7.5) yields the treelevel non-local HQET Lagrangian: We then take the second variation of the non-local action with respect to the relevant fluctuating fields, again following the prescription described in appendix B.3 for the gluons: 9g) Again, only relevant terms are kept here. Since the Lagrangian is real, we expect which is explicitly satisfied by the expressions in eq. (7.9).
In order to evaluate this functional determinant, we row reduce the matrix and obtain the one-loop HQET action S (1) Next, we use eq. (7.9) to derive This functional trace can be converted into integral expressions using the method described in appendix B.2.2 (since we do not need any operators involving the field strength). This

JHEP06(2020)164
procedure yields where we have implicitly defined the integral I α HL in the last line. Evaluating it yields and hence S

Heavy-heavy current matching
The last matching example we will present is the derivation of the coefficient for the heavyheavy current. The steps are essentially identical to the previous calculation, but the details are a bit more involved here since there are two flavors of heavy quarks to track. Note that while we will allow the two heavy quark masses m 1 and m 2 to differ, we will make the simplifying kinematic assumption that v 1 = v 2 = v. The velocity labels on the fields are thus to be taken as flavor indices to keep track of masses only, and not to be taken as arbitrary. Otherwise, various simplification relations, e.g.
In appendix D, we work out the generalization allowing v 1 = v 2 . We start with the UV Lagrangian: where J ± α and J ± 5α are sources for the vector and axial heavy-heavy currents respectively. To reduce the clutter and absorb the phase, we introduce the shorthand

JHEP06(2020)164
where ∆m = m 1 − m 2 . Note that due to the gamma matrices, their conjugation relation is A new feature of the heavy-heavy current matching calculation is that now the two heavy quarks mix: For our purposes here, we only need to solve these equations to linear order in J: Plugging this solution into the UV Lagrangian, and taking the relevant second variations yields with again irrelevant terms for the heavy-heavy operator are discarded in Γ µ,a i andΓ µ,a i C µν,ab = η µν D 2 ab − 2 U µν,ab 21e)

JHEP06(2020)164
Note that this matrix takes the same form as eq. (7.8), so the row reduction is identical. The functional determinant is then Our notation (1 ↔ 2, J + ↔ J − ) here (as well as later in the paper) represents only one term, which results from exchanging both 1 ↔ 2 and J + ↔ J − in the previous term.
Applying the CDE prescription in appendix B.2.2 to evaluate these functional traces, we get This can be simplified down to where the loop integrals are which agree with eq. (1.8) and eq. (1.9). We note that the case of heavy-heavy matching at finite recoil (v 1 = v 2 ) is provided in appendix D. The next section turns to calculation of the RGEs, focusing on the particular examples of the Wilson coefficients for two operators that appear at O(1/m Q ).

HQET RGEs
In this section, we will derive the RGEs for the two operators that appear at subleading order in the 1/m Q HQET expansion: The Wilson coefficient c kin is related to the leading order kinetic term through RPI, and as such it will not be renormalized; we will derive this explicitly at one-loop in this section. The Wilson coefficient for the chromomagnetic moment operator c mag does run, and has phenomenological consequences such as predicting the mass splitting between the ground-state vector and pseudoscalar mesons that contain a heavy quark. A precision determination of this mass splitting requires evolving c mag between the bottom mass and the charm mass. It is straightforward to match onto these Wilson coefficients at tree-level by expanding eq. (2.9). Note that becausē we can replace D µ ⊥ with D µ when it is multiplied by σ µν : These results serve as the boundary conditions for integrating the RGEs whose derivation are the subject of what follows.

JHEP06(2020)164
We follow the procedure outlined in section 4.3. The first step is to compute the 1PI effective action. We need the second variation of the tree-level action with respect to A µ , h v andh v . This yields an equation of the same form as eq. (6.5). As opposed to the residue difference calculation, which uses the non-local form of the Lagrangian, here we are computing the 1PI effective action for HQET directly. Therefore, we must Taylor expand in 1/m Q and truncate. In particular, we only keep terms up to order 1/m Q , and we drop everything that does not include quark fields: Note that we have used a non-standard notation in the above -a bracket with a subscript x is to indicate that the covariant derivative is closed. This is in contrast with the other terms, where the covariant derivatives are open. A detailed elaboration of this notation is given in appendix B.1.5, in particular eq. (B.19). Reusing the row reduction result eq. (6.7), we follow the same intermediate steps given in eq. (6.8) to obtain the one-loop effective action Γ (1) where again we emphasize that we are not including higher order 1/m Q terms than above. Next, we use the expressions in eq. (7.33) to derive Γ HQET ⊃ i Tr

JHEP06(2020)164
where we have used with C F and C A denote the Casimir factor for the fundamental and adjoint representations respectively. Note that we defined the six terms in eq. (7.35) as Γ i=1,··· ,6 in the order that they appear. Since we want to derive the RGEs, we need to isolate the UV divergence. To this end, we regulate the IR using a non-zero gluon mass m 2 , which is sufficient for any of the loop integrals we encounter below. For the first two terms we use eqs. (B.81b) and (B.83a) to obtain For the third term, none of the results in appendix B.4.2 directly apply, so we provide some explicit steps: In the above, the loop integrals V 1,2 and V µν 2,2 can be found in appendix A.4. For Γ 4 and Γ 5 we can use eq. (B.81c) and eq. (B.83b) to find

JHEP06(2020)164
For the sixth term, we again provide the details: Summing all the six Γ i yields the one-loop level effective action: Note that in the last line, we have dropped the operatorh v (iv · D) 2 h v since it is redundant due to the equations of motion. Combining this result with the tree-level effective action Γ (0) where in the second line, we have performed the field redefinition

JHEP06(2020)164
to canonically normalize the tree-level kinetic termh v (iv · D) h v . Finally, we can read off the RGE equations: These results agree with the literature as summarized in eq. (1.12) above.

Conclusions
This paper provides the first application of functional methods augmented by the covariant derivative expansion to a kinematic EFT. In particular, we studied HQET -a description that emerges in the limit that a Dirac fermion is very heavy compared to the scale being probed by the propagating fluctuations. We focused on the particular example of heavy quarks coupled to QCD, and reproduced a variety of matching and running results that had been previously computed using Feynman diagrammatic methods. A critical component was the existence of an operator valued set of projectors that are used to derive an equation of motion for the long distance modes that is valid to all-orders in the heavy mass expansion. Furthermore, the point at which RPI becomes manifest was emphasized. The primary result of this work was the matching master formula given in eq. (5.4). There are two clear directions for future progress. First, one can now use these efficient methods to take one-loop matching calculations to higher order in the heavy mass limit. This would be relevant for high precision applications to either experiments such as LHCb and Belle II, or to theory explorations comparing with lattice QCD calculations in the heavy mass limit [38,39]. Another interesting direction would be to understand how to use these methods for other kinematic EFTs, such as SCET, nrQCD, and others. Furthermore, these methods provide a compelling impetus to revisit the issue of higher-order functional integration directly in the path integral. If the relevant technology could be developed to the stage where two-loop (or higher) calculations could be implemented in a straightforward manner, 16 then these techniques could contribute to extend precision in α s as well. This paper makes the benefits of the functional approach to HQET clear, and furthermore opens the door to using these methods across a broader class of EFTs than had been previously explored.

JHEP06(2020)164
for Astro-and Particle Physics (MIAPP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2094 -390783311. MF thanks the Aspen Center for Physics, supported by National Science Foundation grant PHY-1607611, for hospitality while parts of this work were carried out.

A Loop integrals
In this appendix, we provide a number of frequently encountered loop integrals that are relevant for the calculations performed in the main text. All loop integrals are regulated with dimensional regularization d = 4 − 2 , and parameters are renormalized using the MS scheme.

A.1 Combining denominators
A general loop integral contains multiple propagators in the integrand. In order to facilitate the use of the general results presented in what follows, we use the Feynman parameter technique to combine propagators into a single term: The first integral in the above expression is convenient for combining propagators that are quadratic in the loop momentum. For situations where propagators are linear in the loop momentum, which frequently occurs for HQET calculations, it is more convenient to use the second integral form of this trick. Note that it is straightforward to obtain this last line from the more commonly used form given in the first line by making the variable change x k≥2 = x 1 λ k≥2 and then integrating over x 1 using the delta function.

A.2 Single scale relativistic integrals
After combining the propagators using eq. (A.1), if the numerator of the integrand has factors of the loop momenta with open Lorentz indices, one can convert them by relying on the symmetry of the integral:

JHEP06(2020)164
This reduces all loop integrals (with at least one relativistic propagator) to the following essential form with a single scale ∆: This integral is straightforward to evaluate after performing a Wick rotation to Euclidean space. Then to derive MS expressions, one expands these integrals for small , identifies the 1/ poles and cancels all the occurrences of 1/ − γ E + ln 4π using counterterms, where γ E is the Euler-Mascheroni constant. For completeness, we provide explicit results that are frequently used here. For divergent integrals, we use an arrow to indicate that contributions from counterterms defined in the MS scheme have been added. The integrals for n − m = 1 are quadratically UV divergent: For n − m = 2, the integrals are logarithmically UV divergent: For n − m = 3 (and greater), the integrals are UV finite: and so on.

A.3 Two scale relativistic integrals
Here, we explicitly discuss a commonly appearing integral in relativistic theories that depends on two scales:

JHEP06(2020)164
Note that when this integral contains a nontrivial numerator p µ 1 · · · p µs , one can use eq. (A.2) before combining the propagators, to reduce it to the above form. We evaluate this integral as follows: In the first line above, we have combined the propagators using eq. (A.1). In the second line, we have evaluated the single propagator loop integral using the result in eq. (A.3). In the last line, we have performed the Feynman parameter integral which yields a Gaussian Hypergeometric function 2 F 1 (a, b; c; z). If the arguments a, b, and c are integers (or half integers), the Gaussian Hypergeometric function can be expressed in terms of elementary functions. A non-exhaustive set of explicit evaluations are listed below for some frequently encountered integrals. Note that I r n,k (m 2 , M 2 ) = I r k,n (M 2 , m 2 ), so we will only provide results for integrals with n ≤ k. Taking n + k − r = 2 yields the log divergent integral where for convenience we have defined For brevity, we will we simply state the results for the other evaluations of use here. Taking n + k − r = 3 yields UV finite integrals proportional to 1/M 2 (A.12a) (A.12c)

A.4 Integrals with relativistic and linear propagators
A class of integrals that commonly appear in HQET calculations are built from a relativistic massive propagator (the mass is often included to regulate the IR) and an HQET propagator: (A.14) Note that for this integral, we cannot apply the numerator reduction formulas in eq. (A.2) before combining the propagators, because more tensor structures are available due to the presence of v µ . However, the same logic obviously applies, and we simply need to contract the integrand with both η µν or v µ to determine the possible forms. For example, Results for numerators with more factors of p µ can be reduced in the same way. Applying these relations, we can eventually reduce the integral in eq. (A.14) to the trivial numerator JHEP06(2020)164 case with r = 0, which we can evaluate as follows: 16) where in the second line we combined the propagators using the λ version of the Feynman parameters in eq. (A.1), and then replaced λ → 2mλ for convenience. In the third line, we made the replacements µ ≡ p µ + λmv µ and ∆ = λ 2 m 2 . In the fourth line, we integrated over ∆, and to derive the final result, we used eq. (A.3) to obtain the last line. As a cross check, we see that taking s = 0 in this expression agrees with eq. (A.3): For reference, we provide a non-exhaustive set of explicit formulas with non-zero values for s: where as before, the arrow in eq. (A.18a) denotes that we have included the counterterm contributions using the MS scheme.

B Covariant derivative expansion
In this appendix, we give a comprehensive review of the algebraic approach to evaluating functional traces (or determinants) known as the Covariant Derivative Expansion (CDE). A typical object of interest is the one-loop correction to the 1PI effective action which derives from a Lagrangian of the form where φ is a field, D µ is a covariant derivative, V (φ) is the potential encoding additional interactions, and U = d 2 V /dφ 2 for this explicit example. Here we will present the formalism capable of computing generalizations of eq. (B.1) such as 17 where U i (x) can be any functions of interest. This technique was originally invented in the 1980s [2][3][4], and recently reintroduced in the context of modern EFT calculations in ref. [1], followed by the generalization given in ref. [5]. This technology has been utilized to compute universal one-loop matching results for relativistic theories [7][8][9][10]40]. Our review of the CDE methodology is organized as follows. In appendix B.1, we provide a number of useful tools and tricks for manipulating covariant derivatives and functional traces. In appendix B.2, we review three different algebraic approaches to evaluating functional traces, contrasting them against each other, and highlighting the features of two formulations of the CDE. Then in appendix B.3, we show how the background field method can be used to take functional variations with respect to gauge bosons in a way that manifestly preserves gauge covariance. Finally, we provide a compilation of useful results obtained by applying CDE in appendix B.4.

B.1 Tools and tricks
In this subsection, we spell out many of the manipulations that are required to perform the calculations utilizing the CDE. This has the additional benefit that it allows us the opportunity to define notation and conventions. Particular emphasis is placed on subtleties that can arise.

B.1.1 Baker-Campbell-Hausdorff formula
In deriving the CDE, we will make frequent use of the Baker-Campbell-Hausdorff formula (see e.g. eqs. (B.57) and (B.68)): Here we have introduced a successive commutator notation A n [B] defined as

B.1.2 Gauge coupling
As is well known, one can redefine the gauge field to move the gauge coupling from appearing within the covariant derivative to the coefficient of the gauge kinetic term. In this paper, we take the normalization such that the covariant derivative is where the generator matrix T a depends on the representation of the field that it acts on: We will mostly suppress the representation subscript φ, as it should be clear from the context. Note that if there are successive covariant derivatives, their representation subscripts must all be the same due to the "covariant" nature of the covariant derivative, i.e., it preserves the representation: For convenience, e.g. see eq. (C.2) below, we define a convenient form of the field strength from the commutator of two covariant derivatives where as usual the field strength is Note that if there are multiple gauge sectors, G D µν is the sum of the relevant field strengths.

B.1.3 Functional traces
We use the notation "tr" to denote a trace over the internal symmetries, and "Tr" to denote a trace over both the internal symmetries and the functional space of the operators:

JHEP06(2020)164
Here the subscript R denotes the internal symmetry representation that is being traced over, which we will often suppress when it is clear from the context. Throughout this appendix, we work in d = 4 dimensional spacetime, except that when we confront divergent integrals we will use dimensional regularization with d = 4 − 2 .
Another remarkable property is that the functional trace Tr is invariant under cyclic permutations of operators/functionals, but this is generically not the case for the internal trace tr, because it only traces over a subspace of the functional. However, when all the objects in an internal trace tr are functions (i.e. they are all diagonal in the basis |x , see more elaborations in appendix B.1.5), the cyclic permutation rule holds, see e.g. eq. (B.26a) below.

B.1.4 Product rule for covariant derivative
We highlight an important feature of the covariant derivative D µ -the product rule takes the same form as for the partial derivative when it is acting on fields: This can be justified by noting that the direct product of the two fields A i ⊗ B j forms a tensor representation, for which we have the generator relation Suppressing the representation subscript, this reads and leads to eq. (B.12).
Although it might be clear, we emphasize that eq. (B.12) can be contracted with any non-field constants (that do not transform under the symmetry group in consideration). This leads to some familiar results such as For the last equation above, one can also check the non-trivial consistency between the last two expressions, where in the final form D µ is written out explicitly noting that φ † T a φ φ forms an adjoint representation. A similar but maybe slightly less intuitive consequence is that the covariant derivative can be pulled out of a trace over the internal symmetries: (B. 16) JHEP06(2020)164

B.1.5 Territory of covariant derivative
In what follows, we will be extensively manipulating functionals -typically referred to as operators in quantum mechanics. Note that a special class of functionals are just functions. For example, in quantum mechanics the HamiltonianĤ is built out of two functionals, the kinetic termT = −∇ 2 /(2m) and the potentialV (x). This latter functional is of this special kind, in that it is just a function with the propertyV (x)|x = V (x)|x . Because of this, the territory of a covariant derivative D µ could be ambiguous. When an object like D µ U appears, where U = U (x) is a function, it is not always clear whether this is a function with D µ acting only on U (x), or a non-trivial functional with D µ acting on everything that appears to the right. Usually, one would use a parenthesis (or a bracket of any kind) to remove this ambiguity. However, the calculations to be presented are cumbersome enough that brackets are often used for grouping functionals, as opposed to specifying the territory of D µ . As an explicit example, let us consider the following two toy Lagrangians where we have introduced a vector field U a µ transforming in the adjoint representation. The two interactions are different. In the first Lagrangian, D µ acts only on U a µ and the whole thing D µ U a µ is a function, which can be treated as a field (an adjoint scalar). But in the second Lagrangian, D µ is acting on both U a µ and φ. In the above, we successfully used the parenthesis to distinguish the two scenarios. However, in calculating the 1PI effective action using eq. (B.1), we will get Note that in the second line, a bracket has been introduced to group the three terms, but we do not mean to "close" the territory of D µ (i.e. to restrict the territory of D µ to U a µ only). In fact, this is precisely an example of a term with an "open" covariant derivative, a universal evaluation of which we will provide in eq. (B.84b) up to mass dimension four. This is in contrast with the case of the first line, where the parenthesis is intended to indicate "closing" the territory of D µ . Needless to say, carefully tracking this distinction is critical to our ability to calculate with functional methods. Practically, when we evaluate functional traces with the CDE, an open covariant derivative will get shifted due to e ip·x iD µ e −ip·x = iD µ + p µ (see eq. (B.46)), while a closed covariant derivative does not. Therefore, to remove this ambiguity in the meaning of brackets, we put an additional subscript "x" on the brackets when they are used for specifying the territory of the covariant derivatives. This was explicitly done in the last expression of the first line above. 18 18 When there is no ambiguity in the meaning of brackets, we will drop the explicit subscript x and go back to our usual way of addressing this issue. In particular, for Lagrangian expressions such as eq. (2.6) and final results such as eq. (6.12), the subscript x is often dropped. However, note that we have carefully kept all the subscript x explicit in the results presented in appendix B.4.2.

JHEP06(2020)164
We emphasize that a covariant derivative D µ = ∂ µ − igG a µ T a has two parts -the partial derivative and the gauge fields. Therefore, the meaning of closing its territory by (D µ U ) x is twofold. First, the partial derivative is only acting on U ; second, the generators T a are in the representation of U : where every functional is projected onto the identity function 1(x). Since this notation is non-standard, and the second aspect above can be easily overlooked, we provide a few explicit examples. First, when a D µ is followed by a field strength G D ρσ (acting on a field φ(x)), we have Note the difference in the representations of T a . Second, since the field strength is related to the covariant derivative through G D µν = [D µ , D ν ], it also has a notion of territory. So similar to the case above, for G D µν we also have To further demystify this notation, we can use the product rule in eq. (B.12) to derive the following operator/functional equation for a generic functional A: This shows that our notation defined in eq. (B.19) can be written as a commutator: In addition, one can repeatedly use this relation to derive expressions like This will be extensively used in our CDE derivations (see e.g. eqs. (B.57) and (B.68)). Again, because G D µν = [D µ , D ν ], there is a similar expression for the field strength: A potentially confusing consequence is

JHEP06(2020)164
Note that in the first line above, we are allowed to use the cyclic permutation property for the internal trace, because both objects G D µν (x) and U (x) are functions. On the other hand, the trace in the second line is generically non-zero. For example, taking U (x) to be another field strength, we have where d R is defined as tr R T a T b = d R δ ab . This term will be used in appendix B.2 as a benchmark result for contrasting the three algebraic approaches to evaluating functional traces (see in particular eq. (B.32)). Some additional non-intuitive manipulations are

B.1.6 Integration by parts for covariant derivative
Another useful property of the covariant derivative is integration by parts: To derive this relation, we first pulled the covariant derivative out of the trace using eq. (B.16), and then used the fact that tr(U ) must be a group singlet if this term appears in the Lagrangian. Eq. (B.29) implies

B.2 Algebraic evaluation of functional traces
In this subsection, we present three algebraic approaches to evaluating functional traces: the Partial Derivative Expansion (PDE), the Simplified CDE, and the Original CDE. Presenting the PDE approach will highlight the benefits of CDE. Since both versions of CDE are used in the text (matching uses simplified CDE and running uses original CDE), we will present them both in detail here. For concreteness, we will contrast the three approaches by evaluating the simplest functional trace T 0 in each framework up to mass dimension four: For notational simplicity, we will suppress the representation subscript R, with the understanding that the derivation holds for arbitrary representations. A catalog of results for more involved functional traces are given in appendix B.4.2. JHEP06(2020)164

B.2.1 Partial Derivative Expansion
To begin, we will derive the Partial Derivative Expansion (PDE). This is the brute force approach to evaluating functional traces. The idea is to ignore the covariant derivative structure, i.e., by treating partial derivative and the gauge fields independently, and to perform an expansion in terms of the gauge fields. This yields a large number of terms that are not gauge invariant individually, but that all combine into gauge invariant quantities at the end. In fact, PDE is nothing but the Feynman diagram approach represented algebraically. In particular, since each term generated in the PDE step corresponds to the contribution of a Feynman diagram, these approaches share the disadvantage that they lack manifest gauge invariance, which requires dealing with a large number of terms at intermediate steps. We review this approach to establish that any functional traces can be evaluated purely algebraically, without invoking any of the additional tricks required to derive the CDE. This already demonstrates the benefit of functional methods: calculations are organized into mindless algebraic expansions, which automatically takes care of all the Feynman rules, relative signs, symmetric factors, etc. . In later subsections, we will derive the CDE approaches, which are dramatically streamlined by comparison. As a simple example, consider a gauge theory with an associated coupling g. This implies there is a covariant derivative: In order to compute T 0 up to mass dimension four, we should expand the argument of the trace up to four powers in G µ : (B.36b) 19 Here we have used JHEP06(2020)164 We see that even for this simplest functional trace T 0 , the PDE generates many terms: Operator Terms It is straightforward to evaluate each term through the repeated insertion of the functional identity element

JHEP06(2020)164
For concreteness, we will provide some details for calculating part of T 2G 0 . It contains five terms 20 In order to demonstrate every detail, we provide the evaluation of T 2G1 0 : To evaluate this expression we used the definition of functional trace in eq. (B.11) to derive the second line, followed by insertions of the unit operator to derive the third line. In the fourth line, we extracted the eigenvalues, and used the fact that x|p = e −ip·x . In the last line, we make a convenient change of variables p 2 = p 1 + k. We are then left with a "loop JHEP06(2020)164 which gives Note that T 2G 0 is still not gauge invariant. We need terms with higher powers of the gauge fields T 3G 0 and T 4G 0 , which can be evaluated using the same procedure. A quite tedious evaluation yields Putting everything together, we obtain the final gauge invariant result in eq. (B.32) where we used eq. (B.9) in the last line. This completes our example for using the PDE to evaluate a functional trace. It should now be clear that this approach is very tedious, due to the many terms that are generated at intermediate steps. This can be traced back to the fact that we split the covariant derivative into a partial derivative and the gauge fields to perform the calculation, so gauge invariance appeared broken until all the pieces were finally assembled. Seeing that all the terms ultimately combined together to yield a very simple gauge invariant answer motivates finding a method for evaluating the functional traces without splitting the covariant derivative -this is the covariant derivative expansion.

B.2.2 Simplified CDE
Having motivated the desire for a covariant approach to evaluating functional traces, we now turn to the most straightforward implementation of a covariant derivative expansion, which we call "simplified CDE." This is the CDE method introduced in [5], which is used JHEP06(2020)164 to derive the matching results in the main text of this paper. This will also set the stage for a modified approach -"original CDE" to be presented in the next subsection.
Again, we evaluate T 0 as a concrete example. The CDE relies on manipulations like the following: For this derivation, the key step is going from the first line to the second line, which relies on the fact that x |f (x µ ,p ν )| p = f (x µ , i∂ ν ) x|p . After this step, the integrand in the second line is a function of x and p, and the partial derivative in D µ = ∂ µ − igG a µ (x) T a is interpreted as an operator that acts on any function of x to its right. In this example, this function of x is simply e −ip·x . We then use the fact e ip·x iD µ e −ip·x = iD µ + p µ to obtain the third line. Then in the last line, we have flipped the sign of the loop momentum for future convenience. Again we emphasize that the integrands in lines three and four are no more mysterious than simply being functions of x and p; the territory of the partial derivative ∂ µ (contained in D µ ) is closed by the square bracket, to the end of which there is an implicit identity function 1(x). However, this is not the same as the action implied by our notation defined in eq. (B.19), because we are not closing the territory of the generator T a contained in the covariant derivative D µ -its representation is determined by the representation R (which is suppressed in the expressions above) of the trace tr being evaluated.
Next, one can express the integrand as an expansion in the covariant derivative D µ : Truncating eq. (B.46) up to the fourth power of D µ , we obtain Note that the string of covariant derivatives, e.g. D α 1 D α 2 D α 3 D α 4 , are acting on the identity, which projects out the gauge fields thereby breaking the manifest gauge covariance. If one were to carry this out for every term in the first line of eq. (B.47), the resulting manipulations would simply revert to the PDE. Instead, the CDE approach treats these terms abstractly as D µ , and manipulate the covariant derivatives directly. This yields significant organizational improvements. Not only does it group many terms together, it also allows us to drop quite a few terms before evaluating the loop integral (another way of seeing that they would all conspire to cancel). Take for example the second line of eq. (B.47), where we have dropped terms that do not involve any covariant derivatives, as well as terms that are odd in the loop momentum p since these vanish once the loop integral is carried out. Additionally, simplifying the numerator using eq. (A.2), we can further reduce the remaining seven terms down to only four: Rewriting the last term using we get

JHEP06(2020)164
Now we are ready to evaluate the loop integrals. Note that the D 2 , D 4 , and D µ D 2 D µ terms do not yield gauge invariant objects after acting on the identity function 1(x), so their coefficients must vanish. We can see this magic happen explicitly by evaluating Here we have used the loop integral notation defined in eq. (A.3) and integrated results given in appendix A.2. So the only term in eq. (B.51) that survives is the last line, which gives This reproduces eq. (B.32). This shows how to perform calculation using the most naive implementation of a CDE. The benefit of this approach is that it is very straightforward. In particular, if we are not interested in operators that involve a field strength, we can discard many terms at intermediate steps, which dramatically simplifies the calculations. The downside is that it does generate terms that are not manifestly gauge invariant, but these do not ultimately contribute. As we saw, checking that the coefficients are zero requires tracking cancellations among various loop integrals, e.g. eq. (B.52). It would be desirable if this could be avoided, which is the purpose of the original CDE method presented next.

B.2.3 Original CDE
Having worked through the PDE and simplified CDE approaches, we will now explain the "original CDE" proposal for performing functional traces that was developed in [2][3][4] and recently reviewed in [1]. The simplified CDE generates terms that are not manifestly gauge invariant, which can be shown to not contribute. One can understand the origin of this attribute by studying the origin of the loop integral cancellations in eq. (B.52). We notice that integration by parts lets us organize the integrand into a total derivative of the loop momentum. For example, which is a surface term that integrates to zero. This is why the non-gauge-invariant terms in the naive CDE approach do not contribute. This CDE approach discussed next is JHEP06(2020)164 formulated to ensure that these contributions do not appear at intermediate steps, thereby further streamlining functional calculations. Practically, it is more efficient for extracting the matching coefficients at higher mass dimensions, where multiple insertions of the field strength could appear in many of the terms. We start by performing the same steps that led to eq. (B.46), followed by the pair of operator insertions that are marked in red: (B.55) These expressions are equal. The factor inserted on the right e −iD· ∂ ∂p acts on the constant unit function 1 p thereby evaluating to the identity. The factor on the left e iD· ∂ ∂p also evaluates to the identity since Taylor expanding it only generates total derivative terms in p that vanish upon integration. As we will see, sandwiching our integrand between these two factors will eliminate all the non-gauge-invariant terms that appeared when using the simplified CDE.
To understand how this works in detail, we note that the operator being sandwiched is a function of the combination f (iD µ − p µ ). As the two inserted factors are inverse of each other, they can be brought inside this function such that they act directly on the argument: e Next, we manipulate this into a more useful form using the Baker-Campbell-Hausdorff formula given in eq. (B.4):

JHEP06(2020)164
In the above, we have introduced ∂ α as a shorthand for ∂ ∂pα . We will adopt this notation from now on unless otherwise specified, since the position partial derivative ∂ ∂xα will always be contained in D α and as such will never explicitly appear. Note that to obtain the last lines in the above two equations, we have used D α , G D µν = D α G D µν x repeatedly, see eqs. (B.20b) and (B.24). Combining these two equations, we get Here we have introduced the "CDE field strength" function G CDE µν (x). The benefit of the additional insertions in eq. (B.55) is that they systematically converted all covariant derivatives into commutators, G CDE µν (x). Furthermore, since this object is an explicit function of x, there are no gauge covariance breaking terms appearing due to D µ acting on the unit function. The price is that we now have to keep track of the momentum derivatives that appear in eq. (B.59). In other words, we have systematically traded the position derivatives for momentum derivatives, which have the feature that they do not have any impact on the gauge transformation properties.
Putting this all together, eq. (B.55) becomes The next step is to perform the CDE expansion, which requires performing two simultaneous expansions. First, we expand the "CDE propagator" which is a Taylor expansion in small G CDE µν (x). This expansion is performed up to the mass dimension of interest, noting that each G CDE µν (x) contributes mass dimension two or more. Second, G CDE µν (x) itself should be expressed as a series using eq. (B.59), which should also be truncated up to the desired mass dimension.
For our example, we keep all terms up to mass dimension four:

JHEP06(2020)164
Having worked out this simple example in detail, we will briefly comment on evaluating more involved functional traces. For example, we frequently encounter functional traces of the following form First performing the simplified CDE step gives Then converting this into an original CDE expression using the additional insertions as in eq. (B.55), we have Note that the function U k (x) should also be sandwiched by the insertions, which defines the series The evaluation of several frequently used generic functional traces is given in appendix B.4.2.

B.3 Functional variations for gauge bosons
In the previous section, we showed how to use CDE to evaluate functional traces while maintaining manifest gauge covariance. In doing so, we assumed that the functional trace is already expressed in terms of the covariant derivative D µ . This is the case if the variation is being taken with respect to a field that is not a gauge boson. In HQET, we are typically interested in loops involving gluons. In order to take the functional variation with respect to the gauge fields, we would naively be forced to split up the covariant derivative, which implies a loss of the manifest gauge covariance. By using the background field method (see e.g. [55,57,58]) we can avoid this issue. The first step is to decompose the gauge field into a background field component G a Everything is manifestly gauge covariant with respect to the background gauge field G B,µ , i.e., every term is composed of either D B,µ or G a B,µν . Factors of A µ appear at intermediate steps. However, they will be integrated over perturbatively, which requires gauge fixing. The key to the background field method is that we can implement a gauge fixing condition for A µ that is manifestly gauge covariant with respect to the background gauge field G B,µ : Note that both the fluctuating part A a µ and the ghosts c a andc a transform under the adjoint representation of the background gauge symmetry. Once we have preformed a functional variation, we will drop the subscript B from D B,µ and G a B,µν .

B.4.1 Expansions up to mass dimension six
The definitions and values of the loop integrals in this expression can be found in appendix A.3.
Relativistic and linear propagators. A few functional traces containing both relativistic propagators and HQET propagators are important for the calculations in the main text.

JHEP06(2020)164
been studied extensively [7-10, 40, 42]. Note that the second determinant contains an open covariant derivative, and we have provided the evaluation truncated at dimension four. This result was previously unknown as emphasized in [10,41,42].

C More on RGEs with CDE
In section 4.3, we discussed how one can obtain RGEs using functional methods, first introduced in ref. [5]. The purpose of this appendix is to provide a pedagogical introduction to the methodology as applied to RGEs. First, we discuss a minor complication when computing the RGEs for gauge couplings in appendix C.1, and emphasize how to maintain the manifest gauge invariance in this case. Then in appendix C.2, we provide detailed calculations that reproduce the RGEs for some classic examples: the real scalar φ 4 theory in appendix C.2.1, a theory with two real scalars φ and Φ in appendix C.2.2, and gauge theories in appendix C.2.3.

C.1 RGEs for gauge couplings
The procedure described in section 4.3 applies to any coupling λ. However, there is a minor complication when considering a gauge coupling g. If the 1PI effective action is written in terms of gauge invariant operators, 21 such as what we will obtain from a CDE procedure, there will be no explicit operator O g -the gauge couplings are all hidden in the covariant derivatives. A brute force approach would be to abandon the manifest gauge invariance by splitting a gauge invariant term and isolating a part of it as a candidate O g . In doing so, one typically gets multiple terms that can serve as O g . For example, the QCD Lagrangian is There are no explicit factors of the gauge coupling in the first line, while expanding the gluon field strength and covariant derivatives exposes the g s dependence. Clearly, there are three terms that could be identified as O g . This is not a problem -one can choose to work with any of these three operators since gauge invariance will ensure that the results are the same. However, there is another approach, which avoids breaking up the gauge invariant terms. This is additionally convenient here since the CDE approach yields a form of the 1PI effective action that is already packaged in this way. 21 Note that this is not generically guaranteed. For example, in the traditional computation of the 1PI effective action for a non-Abelian gauge theory, the result will not be composed of gauge invariant operators due to the gauge fixing procedure, e.g. see chapter 16.5 of Peskin and Schroeder [57]. However, when using functional methods, we can benefit from background field gauge fixing eq. (B.74) described in appendix B.3. With this choice of gauge fixing, the 1PI effective action will be expressed in terms of gauge invariant operators, as emphasized in e.g. chapter 16.6 in Peskin and Schroeder [57].

C.2 Example RGE calculations
In this section, we provide detailed calculations of RGEs for a number of well known examples.

C.2.1 Scalar φ 4 theory
We start with the simplest example, φ 4 theory for a real scalar field φ. The Lagrangian is We will demonstrate how to compute the RGEs for the couplings m 2 and λ. Starting with the Lagrangian in eq. (C.8), we take the second variation with respect to φ: Making use of the universal result eq. (B.84a) with U = U φ = 1 2 λφ 2 , we see that the one-loop 1PI effective action contains In the second line above, we have only kept terms up to mass dimension four. Combining it with the tree-level 1PI effective action Γ (0) The kinetic term is already canonically normalized, so no additional field redefinition is required. We then simply use eq. (4.12) to find the RGEs:

C.2.2 Heavy-light scalar theory
As a second example, we consider a theory with two real scalar fields φ and Φ, with masses m and M respectively. This toy model was analyzed extensively using Feynman diagrams in [59]. The Lagrangian is

JHEP06(2020)164
The trace over gluons gives the adjoint Casimir tr G T a T b = N c δ ab . (C.20) We denote the quark piece as q and assume that all N f flavors transform under the same representation Q of SU(N c ). Accordingly, we have tr q T a T b = N f d Q δ ab , (C. 21) with the factor d Q depending on the representation Q. For example, d Q = 1 2 if Q is in the fundamental representation, and d Q = N c if Q is in the adjoint representation.
We will take variation with respect to the gauge fields using the background field method described in appendix B.3. The Lagrangian is Taking the second variation, we obtain (C.23) where we have defined C µν,ab = η µν D 2 ab − 2gf abc G µν,c = η µν D 2 + 2G D,µν ab , (C.24a) B = i / D , Γ µ,a = gγ µ T a q ,Γ µ,a = gqγ µ T a . (C.24b) To derive the RGE for the gauge coupling, we only need to keep the kinetic term of the gauge boson. Dropping all term with quark fields, we find (C.25) The one-loop effective action is then holds for the results obtained in ref. [43], as explicitly stated in eq. (A3) therein. So this verifies our derivation of C 2 (z, w) as well. We note that our coefficientC 1 (z, w) defined above is slightly different from the C 1 defined in ref. [43], because in our approach, the residue difference contribution needs to be added on top of the result obtained in eq. (D.7), cf. eq. (1.8). We leave the evaluation ofC 1 (z, w) and I α HH,5 as well as the comparison of them with ref. [43] for future work, as we believe the evidence that functional methods can be applied to HQET is sufficient.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.