New Covariant Feynman Rules for Effective Field Theories

We provide a new and completely general formalism to compute the effective field theory matching contributions from integrating out massive fields in a manifestly gauge covariant way, at any desired loop order. The formalism is based on old ideas such as the background field method and the heat kernel, however we add some crucial new ingredients that greatly improve the simplicity and general applicability of the approach. We formulate our method in terms of Feynman rules, the resulting effective action is expressed in terms of local heat kernel coefficients. We also provide as supplementary material a mathematica code that facilitates the computation of these coefficients.


Introduction
Effective field theories (EFT's) constitute one of the most important tools in high energy physics, in particular in the perturbative regime. The fact that a theory with heavy particles can be represented in the infrared (IR) by a simpler one with only light degrees of freedom and a series of higher dimensional operators greatly simplifies the phenomenological analysis of the theory, providing an accurate way of obtaining low-energy observables directly related to experiment (see ref. [1] for a review). The problem is broken down into two simpler ones, the matching of the full UV theory to the EFT (the calculation of the coefficients of the higher dimensional operators), and the calculation of IR observables directly in the EFT, employing the renormalization group (RG) for improved precision.
In recent years, a lot of effort has gone into studying the Standard Model effective field theory (SMEFT). The full list of dimension-six operators has been obtained in ref. [2], and their complete RG equations have been computed in refs. [3][4][5]. Countless phenomenological studies exist relating the Wilson coefficients to experimental data, putting limits on their values which in turn sharpen our understanding of possible extensions of the Standard Model (SM). This very model independent analysis eventually has to be completed by the first step mentioned above, the so-called matching.
While matching is in general quite straightforward at the tree-level, in many perturbative UV theories not all operators are generated in this leading order, and loop calculations are necessary. Moreover, some of these operators are directly connected to very precise experimental data, so even loop-suppressed Wilson coefficients can be of great value for putting bounds on possible UV theories. Examples include dipole operators mediating µ → eγ transitions as well as contributing to electric and magnetic dipole moments of light particles.
Loop calculations are traditionally performed in a non-gauge covariant approach, breaking up covariant derivatives into partial derivatives and gauge fields. This allows for a simple evaluation of the diagrams in momentum space, but manifest gauge invariance becomes hidden and is ensured only by intricate relations between different diagrams (Ward identities). Diagrams containing a certain number of external gauge fields have to be combined with diagrams with fewer gauge boson lines but more powers of external momenta to assemble into gauge invariant operators. Usually a large number of such diagrams has to be computed in order to obtain all possible operators of a desired dimension.
In contrast to this canonical non-covariant perturbation theory, a covariant alternative has been developed a long time ago [6][7][8][9][10][11][12]. The first manifestly covariant formalism appeared in the seminal work of Schwinger [6] (see ref. [13] for a pedagogical introduction), and was later further pioneered by DeWitt [7,8], who also recognized its applicability to gravity. The approach is based on the background field method in which the effective action is obtained by integrating the quantum fluctuations over a classical background. Within this approach it is possible to maintain manifest background field covariance throughout all calculations. Schwinger and deWitt advocated the use of the so-called heat kernel (HK), allowing for an efficient calculation of the one-loop effective action of bosonic background fields (including scalars, gauge fields, and gravity). The Schwinger-DeWitt technique was considerably generalized for arbitrary one-loop graphs (including fermionic background fields as well as different spins and masses in the loop) in [10]. Most notably, the covariant HK method has been applied to chiral perturbation theory [14] and to the calculation of anomalies in various dimensions [15][16][17][18][19]. The method has also been used in the context of extra dimensions [20][21][22]. The background field method for higher loops was developed in refs. [23][24][25], and has been used in its covariant form to calculate β functions for gauge theories at various loop orders [26][27][28][29] by considering covariantly constant background field strengths. The covariant approach typically reduces by a lot the number of independent diagrams to be computed.
However, despite the decent amount of literature and recent progress on covariant perturbation theory, a completely universal treatment valid to any loop order is still missing. The purpose of the present work is to provide such a simple and universally applicable framework. Our main new idea is to factorize each Feynman diagram into an "n-point function" Γ(x 1 . . . x n ) that is gauge invariant at any of the vertices 1 x i of the diagram, and a momentum space loop integral which depends on n external momenta I(p 1 . . . p n ). The contribution to the effective Lagrangian is then simply given by I(i∂ 1 . . . i∂ n )Γ(x 1 . . . x n )| x i =x , which can be systematically calculated in an expansion of the background fields and their covariant derivatives, making essential use of HK techniques. This simple structure of loop graphs is quite powerful but to the best of our knowledge has never been exploited in the literature. Our method has no restrictions whatsoever for the types of interactions, loop order, types of background fields, loop masses etc. 2 We illustrate our formalism with examples at one and two loops.
This paper is organized as follows. In sec. 2 we review the concept of the HK representation of the propagators. In sec. 3 we present and proof our generalized Feynman rules. In sec. 4 we illustrate our formalism with examples at one and twoloop order. In appendix A we review the background field method and write explicitly the field dependent masses and vertices for gauge theories. In app. B we review the techniques for the calculations of the HK coefficients and tabulate the leading ones. In app. C, we summarize some basic integrals over Schwinger parameters. In the supplementary material to this article we provide a short mathematica notebook that computes HK coefficients with arbitrary number of covariant derivatives.
2 Heat kernel representation of propagators

Scalars
We work in Lorentzian signature. Consider the background-field dependent Feynman propagator for a scalar field in an irreducible representation of the gauge and (if present) global symmetry group: The mass m is field-independent and proportional to the identity, it may be zero. In contrast, the mass X is field dependent, but, since we consider irreducible representations, X does not mix different representations. Representation-mixing field dependent masses will be treated as mass insertions (bilinear vertices). The gauge fixing procedure may generate contribution to X, see eq. (A.9). In the following, we will suppress the i terms. The heat kernel (HK) representation of the propagator is defined as [6,49,50] The HK K s (t, m, X) satisfies the Schrödinger equation 3 Notice that for the case of a scalar, X is Hermitian. However, for nonzero spin (see below), there is a contribution to the effective mass matrix given by the term −S µν F µν , where S µν are the generators for Lorentz transformations which are non-Hermitian in non-Euclidean spacetime signature. Defining K s (t, m, X; x, y) ≡ x|K s (t, m, X)|y , one now makes the ansatz is the free HK (that is, the HK of a gauge singlet with vanishing X). The function B(t) has a regular expansion around t = 0: The b 2n are known as the heat kernel coefficients. In particular, the coefficient b 0 , which is independent of X, is equal to the Wilson line along the straight line segment connecting x to y (see app. B). The higher coefficients have no closed-form expressions, however all we will need are their various covariant derivatives evaluated at coinciding arguments (x = y), which we will refer to as the local heat kernel coefficients (LHKC's). These objects are local polynomials in the background fields and can be computed systematically [7-9, 12, 51]. If the field transforms in the r representation, B transforms as the r representation at x (that is, from the left) and as the conjugater representation at y (from the right), since G α β (x, y) = Φ α (x)Φ † β (y) . K 0 transforms as a singlet. We will be using the HK in a "mixed position-momentum representation" Sometimes an interaction contains a covariant derivative acting on a fluctuation field (this happens, for instance, in a gauge theory with scalars, see app. A). We will associate this derivative to the propagator directly, that is we also will need propagators and similar terms with a derivative acting on y as well as both arguments. Such terms are easy to deal with in our formalism as D µ G s = dt D µ K s (t), and etc. One may evaluate explicitly the integral over the Schwinger parameter t to obtain the propagator in terms of the HK coefficients: (2.10) Even though the Feynman rules can be given equally well in terms of this propagator instead of the HK, we find it more practical to perform the integration over the Schwinger parameters last (in particular, after the momentum integration), so eq. (2.10) is not relevant for us in practice. There is however one important observation that can be made from eq. (2.10). Suppose X is non-Hermitian, which, as we already mentioned, can occur for fields of nonzero spin. The representation of the propagator in terms of the HK in eq. (2.2) is not well defined for arbitrary complex eigenvalues. However, all we are interested in is the local part of the effective action, which is implied in the small-t expansion of eq. (2.6). If we define our propagators in terms of this expanded HK, we obtain the well defined object in eq. (2.10).
2.2 Spin-1 2 fermions We start by writing the fermion propagator in terms of the scalar one: 4 where S µν = i 4 [γ µ , γ ν ] is the Lorentz generator for spin 1 2 . We then define the corresponding HK such that as previously G f = dt K f (t). The HK in its mixed representation that we will be using in our Feynman rules is thus (2.14) We notice that the same function B appears on the right hand side, albeit with a different argument X.

Gauge fields
Finally, consider the gauge sector. It is easiest to use Feynman gauge ξ = 1 to avoid dealing with non-minimal operators. 5 Then from the terms quadratic in the fluctuations detailed in app. A one can read off the propagator In the presence of a field dependent mass Y + iγ 5 Z, we can instead write In most applications Y = Z = 0. 5 See ref. [10] for the HK with non-minimal operators. Like in the spin-1 2 case, the term in X v proportional to the field strength can be understood as −S µν F a µν t a , with the spin-1 Lorentz generator (S µν ) αβ = i(δ µ α δ ν β − δ µ β δ ν α ) and the adjoint generator (t a ) b c = −if ab c . The gauge HK can be decomposed as The gauge and Lorentz indices on G v , K v and B are suppressed here for clarity. Notice that in the limit of vanishing background fields we obtain B a µ b ν = δ a b δ µ ν , recovering the usual propagator in the Feynman gauge, including the correct negative sign.
The HK for the ghost is given by the scalar one, with the field dependent mass that is of the same structure as the one for the gauge fields.

Feynman rules 3.1 Summary of the rules
Let us first give a concise list of rules which can be directly applied to any given problem. We will provide a proof for these rules in sec. 3.2.
1. Draw all diagrams with with no external legs.
2. For each propagator write the expressions given in tab. 1. For propagators of derivatives of fields add additional factors of (−ik µ + D x µ ) and (ik µ + D y µ ) as explained in eq. (2.9). The expressions for the fermion propagator and the derivative scalar propagator assumes that the momentum k flows parallel to particlenumber.
3. Write the field dependent couplings, including the usual factor of i for each of them.
4. To each vertex, associate an ingoing momentum p i , and impose momentum conservation at every vertex. Perform the loop momentum integration in d dimensions.
5. Expand in the external momenta p i to the desired order, replace the momenta p i by the partial derivatives and take the limits x i → x everywhere.
6. Integrate over the Schwinger parameters, remove all poles via M S, and take the limit d → 4. The result, multiplied by −i, is the EFT matching contribution to the Lagrangian.
If there is a field dependent mass that does not mix representations (i.e., that commutes with the field-independent mass matrix), one may wish to include it into the propagators and heat kernels. Field dependent masses that do mix representations need to be included as mass insertions (bilinear vertices).
The partial derivatives appearing in step 5 always act on a total gauge singlet, essentially the product of all the B-factors and field-dependent couplings (see sec. 3.2). Once these derivatives are distributed over the individual (non gauge singlet) factors, covariant derivatives will reappear.
Observe the presence of the coincidence limits x i → x in step 5. This means that, as claimed in section 2, only the LHKCs, i.e., the coincidence limits of coefficients b 2n (and their covariant derivatives) appear. A technical simplification that we already would like to point out is that any string of partial derivatives resulting from the expansion in external momenta is automatically symmetric (since partial derivatives commute) where the parenthesis on the index denotes symmetrization with strength one. This symmetrization will thus appear also in the resulting covariant derivatives, which greatly simplifies the calculation. For instance, symmetrized covariant derivatives of the ubiquitous coefficient b 0 are zero, see app. B.
A last comment concerns one-loop diagrams with zero vertices, that is with one single massive particle in the loop. This case is somewhat special since it involves the log of the propagator, but is also the most studied one in the literature [51]. For bosons one needs the functional trace identity and for Dirac fermions These traces are straightforwardly evaluated in terms of traces over the HK coefficients b 2n (x, x).

Proof of the rules
We now proceed to prove the rules given in sec. 3.1. Rule 1 just states that we are calculating the effective action in the background field method. Lines correspond to propagators of fluctuations and not to background fields, which are hidden in the lines and vertices.
Rule 2 is evident from the HK representation of the propagators G = dtK(t) and the explicit expressions for the HK in eq. (2.7) and eq. (2.14). Notice that we are not writing the factors of e ikx , they will be combined into momentum conserving delta functions similar to the standard (non-covariant) Feynman rules (see below). We also delay the integrations over the momenta and Schwinger parameters.
Rule 3 is self-evident.
To show rules 4 and 5, let us examine the possible spacetime dependence for a given vertex at x resulting from the rules so far. For instance, consider a bilinear interaction where Φ and Ψ are two fluctuation fields (of possibly different spins), and C is a background-field dependent coupling. We suppress possible Lorentz indices and only display the gauge indices α and β, belonging in general to different irreducible representations. Such an interaction will give rise to factors of the kind This expression is gauge invariant at x, the gauge indices α, β are contracted in a gaugeinvariant way in the same way as in the Lagrangian. For this reason, the product of all propagator B factors and couplings C becomes a complete gauge singlet that we call where P denotes the number of propagators and V the number of vertices. Diagrams with only scalar propagators become (suppressing the t i integrations for clarity) where P is the total momentum exiting at the th vertex in the propagators. We have Fourier-transformed the "n-point function" Γ and performed the integration over the x i to generate momentum-conserving δ functions. 6 This Fourier transform of Γ introduces V new "external" ingoing momenta p , one for each vertex. The generalization to fermion propagators and derivative propagators is straightforward: in this case the objects Γ in eq. (3.7) may contain additional derivatives (from terms such as / D and D µ in eqns. (2.14) and (2.9)), and one gains additional factors of propagator momenta / k and k µ . These details however are not important for the following arguments. Using momentum conservation we get for our diagram where L is the number of loops, and the propagator momenta k i are now to be interpreted as the corresponding linear combinations of the loop momenta q r and external momenta p as dictated by momentum conservation. At this time, we can perform the integral over the loop momenta which we call such that (restoring the t i integrations): (3.11) Transforming back to position space shows that this is indeed local upon expansion of I: We note that the integral I is regular at vanishing external momenta, such that a Taylor expansion around zero external momenta is well defined. Eq. (3.12) is precisely what is stated in rule 5, including the coincidence limit x → x. Eq. (3.12) is known as the effective action of the full theory. To really obtain the matching contribution, we need to subtract the EFT contribution to the same operators. However, provided that the EFT does no longer contain any massive fields, all EFT integrals are scaleless, and hence vanish in dimensional regularization. For logarithmically divergent integrals, this is due to a cancellation of UV and IR poles. Following the pedagogical discussion of ref. [1], we make this cancellation visible by defining UV ≡ , and IR ≡ , and write schematically the result of the EFT contribution as (at one-loop for definiteness) where "c.t." stands for counterterm and A is some coefficient. Thus, the net effect of the EFT contribution are some poles in IR . The same poles must occur in the full theory calculation as the two theories by construction reproduce the same IR physics. 7 Thus, in dimensional regularization, the EFT contribution's only effect is to cancel the IR divergences of the full theory (of course, any UV divergences of the full theory are cancelled by its own counterterms). In practice, we can avoid dealing with counterterms by simply using M S for any divergence appearing in the full theory calculation (UV and IR), that is, we apply minimal substraction to all the poles in eq. (3.12). The above procedure works because we have expanded out any IR scales (fields and external momenta) in both the full and EFT theories. In the presence of explicit light masses in the EFT the preceding discussion is still correct as long as we expand in this mass parameter as well (in both the full and EFT contributions) before integrating over the Schwinger parameters. 8 This concludes the proof of the Feynman rules. At this point a few comments are in order • The object Γ is gauge invariant at each of the vertices, and hence the formula (3.12) is manifestly gauge invariant. Once the partial derivatives are distributed over the different factors of B i and C present in Γ, covariant derivatives will appear, since the product rule is valid for covariant derivatives.
• The new external momenta are not the momenta of single background fields, but rather of a product of many of such fields. The momentum conservation however is the usual one that can be read off of the diagram, that is, an external momentum p entering at a vertex has to be equal to the sum of the momenta exiting through the attached propagators.
• The momentum integrations are simple Gaussians, and there are no divergences (neither IR nor UV). However, divergences will appear once the t i are integrated over. UV divergences appear in the t i → 0 region, while IR divergences in the t i → ∞ region. However, if the momentum integrations are performed in d dimensions, the t i integrations will be automatically regulated. If the Schwinger parameters are integrated before the momenta, the UV and IR divergences appear as usual in the loop integrations.
• It is possible to evaluate the Gaussian momentum integral as well as the resulting integrals over the Schwinger parameters without any Wick rotations. However, it is convenient to perform the usual Wick rotation of the momenta together with the following rotation of the Schwinger parameters with real τ i > 0. This is allowed as it precisely leaves a real Gaussian momentum integration as well as an integral over Schwinger parameters that is exponentially suppressed for massive propagators ∼ e −m 2 i τ i .
• It is possible to perform the integration over the Schwinger parameters first as in eq. (2.10), in this case the expansion in the external momenta should be done before the loop integrations. However, we find it more practical to do the t i integrations at the end, for two reasons. Firstly, the Schwinger parameters serve as a tool to efficiently combine propagators (similar to Feynman parameters 9 ) and hence facilitate the treatment of external momenta. Secondly, delaying the expansion of the B factors in t and of the loop amplitude in the externa momenta to the end of the calculation means that the loop momentum integrations are Gaussians that can even be evaluated in closed form in full generality (see app. C.1).

A one-loop example
In this section, we integrate out (at one loop) a massive fermion that couples to a light scalar and a light fermion. Let us also focus specifically on the part of L eff containing two-fermion operators with additional covariant derivatives and field strengths but no scalar fields. In this case, the only relevant interaction at one-loop is where Ψ is the fluctuation of the heavy fermion, Φ the one of the light scalar, and C is a fermionic field-dependent coupling. As an example consider a heavy triplet fermion Ψ a coupled to the SM Higgs Φ i and the top quark doublet q j L , in which case C a i (x) = y(σ a ) ij q j L (x) where y is a real coupling. Applying the Feynman rules to the diagram in fig. 1 gives where the trace is over the gauge indices, and we have defined the shorthands B s ≡ B(X = 0), B f ≡ B(X f ). Performing the Gaussian momentum integration yields We have defined τ ≡ it and σ ≡ is, see discussion around eq. (3.14). Eq. (4.3) is the exact one-loop result to all orders in the operator dimension. To extract individual operators of a desired dimension, we expand the exponentials containing the partial derivative, as well as B f and B s according to eq. (2.6). In this particular example, the term proportional to m will yield odd-dimensional operators (zero for chiral C), while the other terms will yield even dimensional ones. For a nontrivial example, consider the dimension-six operators: where the semicolon denotes covariant differentiation, see eq. (B.6). Making use of the observation pointed out after eq. (3.2), many terms will be zero, the nonzero ones are: The calculation of these dimension-six operators in non-covariant perturbation theory would require the evaluation of 14 diagrams (not counting permutations of external legs), including some tedious Dirac algebra, plus the reconstruction of the covariant operators in terms of the noncovariant result.

A two-loop example
In this section we would like to apply our formalism to integrate out a massive scalar coupled to gauge theory at two loops. As should be clear by now, our formalism contains roughly three main calculational steps: the Gaussian momentum integration, the reduction to local operators, and the integral over the Schwinger parameters. We will see that the first two steps are virtually identical at all loops, in particular, the momentum integration can be given in closed form, see eq. (C.2). The integrals over the Schwinger parameters become more involved at higher loops, in the example provided below they can however still be computed analytically.
Assuming that the heavy scalar does not acquire a vacuum expectation value, we can set the scalar background field to zero, that is, all the bilinear vertices in tab. 2 vanish in any gauge. The two-loop diagrams with at least one heavy propagator are the three sunset graphs and the figure-eight graph shown in fig. 2.
The figure-eight diagram gives where the symmetry factor of 1 2 results from the tadpole loop of the gauge boson. Performing the momentum integration yields x, x) (4.8) Being a scaleless integral the τ integration gives zero in dimensional regularization.
Let us then evaluate the three sunset graphs. These diagrams also illustrate how derivative interactions are treated in our formalism. We mark a covariant derivative of a field by a star. The triple vertex connects to one scalar and one scalar with covariant derivative and hence each vertex carries one star, and the three diagrams represent the different ways of placing the stars in the diagram. The first two diagrams carry a symmetry factor of 1 2 because the two vertices are equivalent, whereas the third diagram has no symmetry factor (the two vertices are not even the same, but rather complex conjugates). 10 Their contribution to the effective Lagrangian is The three traces correspond to the three diagrams. Notice how extra factors of −ik µ + D µ etc. appear due to the propagators of derivatives of fields (denoted by stars in the diagrams), these derivatives are meant only to act on the B factor immediately to the right (in contrast the partial derivative in the exponent which acts on all three B factors). The Gaussian momentum integration is straightforward, see the master formula Eq. (C.2). Focusing for instance on the first diagram, one gets where ω ≡ τ σ+σ + τ . Again this is the all-order result in the operator dimension. The other two diagrams are completely analogous. The next step would be to expand this result up to a desired order in the operator dimension, that is, perform the expansion of the exponentials containing the derivatives as well as the expansion eq. (2.6). There is nothing new to say here, this expansions work exactly as in the one-loop example. The resulting intergrals over the Schwinger parameters can still be performed analytically, see app. C.3.

Conclusions
We have presented new covariant Feynman rules that allow for a simple and efficient calculation of the local effective Lagrangian from integrating out heavy particles. The formalism is universal with no restrictions to the type of interactions (including arbitrary derivative interactions), loop order, or type of particles in the loop (e.g. spin, massive, massless). Even though our main focus was on the EFT matching contributions, our formalism can equally well be applied elsewhere, for instance to the computation of renormalization group functions of general effective field theories.
A possible generalization is the inclusion of a gravitational background, many results on the gravitational LHKC's exist in the literature, see ref. [51] for an overview.
Our algoithm is amenable to full automatization and we have provided a first small step in that direction by including with this publication a mathematica notebook that computes the LHKC's.

A Background field method for gauge theories
In this appendix we will derive the field dependent couplings for gauge theories. We include scalar fields here due to their entanglement in the gauge fixing procedure. Let us thus consider the Lagrangian One decomposes the fields into backgrounds (φ, A) and fluctuations (Φ, A) 11 11 We work with the convention D µ = ∂ µ − iA a µ t a such that [D µ , D ν ] = −iF a µν t a . Furthermore the adjoint generators are (t a adj ) bc = −if abc To eq. (A.1) we add the following gauge fixing Lagrangian and ξ and η are real gauge fixing parameters. One finds the following terms quadratic, cubic and quartic in the fluctuations: For the ghost Lagrangian, one obtains The gauge choice ξ = 1 eliminates the non minimal terms in the gauge propagator. In addition, most often the choice η = 1 is adopted, in order to remove the bilinear term ∼ A µ D µ Φ. However, as is clear from the above Lagrangians, η = 0 comes at a price, as we are generating many new terms in the above Lagrangians. Since our formalism deals very naturally with derivatives on propagators, other choices, in particular ξ = 1, η = 0, seem at least equally motivated. The last line in eq. (A.5) are bilinear vertices, they cannot be eliminated completely for any choice fo ξ and η. From eq. (A.5) and eq. (A.8) one reads off the field dependent masses for scalars, vectors and ghosts:  Table 2. Bilinear vertices for gauge theory coupled to a scalar. Recall that G ia = g(t a ) i j φ j . An asterisk with index ν indicates that the attached propagator is a derivative one, i.e., Φ ;ν Φ † (for ingoing arrow) or ΦΦ † ;ν (for outgoing arrow). The right and left columns are complex conjugate vertices.
1. Propagators of gauge fields that start and end at a the same vertex produce a factor of 1 2 .
2. P equivalent propagators give a factor of 1 P ! 3. V equivalent vertices give a factor of 1 V ! . The first rule only exists for real fields, i.e. gauge bosons in the present context. It is the source of the factor of 1 2 in the last diagram of fig. 2. The other two rules are best illustrated by further examples. Consider the four diagrams shown in fig. 3. Diagram A carries a factor of 1 2 due to rule 2 (the two propagators are interchangeable). Diagram B has two equivalent vertices and two equivalent gauge boson propagators and hence has symmetry factor 1 4 (the two scalar propagators are non-equivalent since the particle number flows in different directions). The presence of the stars (marking derivatives of fields) typically removes some of the symmetry of the same diagram without stars (as appearing e.g. in non-covariant perturbation theory). Diagram C has two equivalent vertices and two equivalent propagators (the ones without the stars) and hence gets a factor of 1 4 . Finally, the last diagram D has 3 inequivalent propagators but two equivalent vertices, and hence receives a factor of 1 2 . As always, these factors can also be computed by counting Wick contractions, yielding the same results. Table 3. Trilinear vertices for gauge theory coupled to a scalar. Recall that G ia = g(t a ) i j φ j . An asterisk indicates that the attached propagator is a derivative one, i.e., Φ ;ν Φ † (for ingoing arrow), ΦΦ † ;ν (for outgoing arrow), A µ;σ A λ , or cc ;ν .

B Heat kernel coefficients
In this section we briefly review DeWitt's recursive procedure for the calculation of the LHKC's. Eqns. (2.3) and eq. (2.5) imply the following differential equation for B(t) Plugging in the expansion eq. (2.6) yields the recursion For n = 0 this amounts to a differential equation for b 0 with the boundary condition b 0 (y, y) = 1 which is a consequence of K(0; x, y) = K 0 (0; x, y)b 0 (x, y) = δ(x − y) and K 0 (0; x, y) = δ(x − y). The solution to this equation is the Wilson line where P denotes path ordering and the line integral is to be taken along the straight line segment connecting y to x. Let us introduce the notation for the coincidence limits, Furthermore let us denote covariant differentiation with a semicolon as follows 12 The coincidence limits for b 2n and its covariant derivatives can be computed from repeated covariant differentiation of eq. (B.2). Let us start with derivatives with respect to the first argument only. Differentiating N times and using the commutator [D µ , D ν ] = −iF µν , one can derive the generalized recursion relation for the LHKC's This relation allows one to compute the coincidence limit of the b 2n with an arbitrary number N of covariant derivatives in terms of the analogous quantities with smaller N and/or n. The cases with up to 2n + N = 4, are easily obtained by hand.
• Dimension-0 coefficients • Dimension-2 coefficients • Dimension-4 coefficients The coefficients with derivatives acting on the second argument can be obtained from Synge's rule, where in eq. (B.29) we also used the Bianchi identity eq. (B.25).

C Loop integrals C.1 Multi-loop momentum integrals
The most general L-loop momentum integral that we need to evaluate in our formalism is the following simple Gaussian integral: The matrices T , S, and U are linear in the Schwinger parameters t i , and in addition T and U are symmetric and positive definite. We have introduced sources ξ i in order to account for possible additional terms q µ i in front of the exponential which can be obtained by differentiating with respect to ξ µ i . As mentioned around eq. (3.14), in addition to the usual Wick rotation of the momenta we can make the rotation t i = −iτ i . In terms of these Euclidean parameters, we have where T ≡ iT , U = i(U − S T T −1 S), and R = T −1 S. The matrices T , U and R are now purely real in terms of the τ i , in addition, T is positive definite.
After the momentum integration and the expansion in the external momenta as well as the expansion eq. (2.6), the remaining τ i integrals to be evaluated are of the form where the a i > 0 and γ ≥ 0 are integers. Notice that det T is a homogeneous polynomial of degree L in the τ i . In the following we will evaluate the general one-loop case and some selected twoloop cases.

C.2 Schwinger integrals appearing at one loop order
At one loop order, there are only bilinear vertices. For P = P 0 + P m propagators (with P 0 massless and P m degenerate massive propagators), there are thus P such vertices. The most general one loop graph leads to the integral where σ ≡ P 0 i=1 σ i , τ ≡ Pm j=1 τ j , a ≡ P 0 i=1 a i and b ≡ Pm j=1 b j . Let us discuss a few special cases. where n is the order of the HK expansion that also determines the dimension of the operator, D = 2n.
• A graph with two bilinear vertices and only (degenerate) heavy fields can be written in terms of (C.6) • A graph with two bilinear vertices and one heavy and one light field is proportional to -25 -

C.3 Schwinger integrals needed at two loop order
Moving on to two loops, let us consider sunset graphs with exactly three propagators (that is, only two trilinear vertices and no bilinear ones). We will need integrals of the type The integrals for the case of three massive propagators are not expressable in a simple way in terms of Γ functions, and we will skip them here. The other possible two-loop topology is the figure-eight graph. The two loopmomentum integrations are independent and thus result in a product of two one-loop integrals f ({a i }, {b j }, c) already evaluated in app. C.2.