All two-loop scalar self-energies and tadpoles in general renormalisable field theories

We calculate the complete tadpoles and self-energies at the two-loop order for scalars in general renormalisable theories, a crucial component for calculating two-loop electroweak corrections to Higgs-boson masses or for any scalar beyond the Standard Model. We renormalise the amplitudes using mass-independent renormalisation schemes, based on both dimensional regularisation and dimensional reduction. The results are presented here in Feynman gauge, with expressions for all 121 self-energy and 25 tadpole diagrams given in terms of scalar and tensor integrals with the complete set of rules to reduce them to a minimal basis of scalar integrals for any physical kinematic configuration. In addition, we simplify the results to a set of only 16 tadpole and 58 self-energy topologies using relations in order to substitute the ghost and Goldstone-boson couplings that we derive. To facilitate their application, we also provide our results in electronic form as a new code TLDR. We test our results by applying them to the Standard Model and compare with analytic expressions in the literature.


Introduction
The Higgs boson mass has been measured to an accuracy of about O(100) MeV, making it an electroweak precision parameter. In the Standard Model (SM), this is used to extract the Higgs quartic coupling λ, which, through a significant amount of work, can now be done with high precision, where all of the relevant running parameters of the Lagrangian can now be extracted from calculations at full two-loop order [1,2] and partially at the three-and four-loop order [3][4][5][6][7][8][9][10][11][12][13][14][15][16]. A code for calculating λ in Landau gauge, SMH [17], and codes for extracting all relevant SM parameters (the gauge couplings, top and bottom Yukawa couplings, and Higgs quartic coupling), mr [18,19] and SMDR [20], exist. As a result of this effort, the uncertainty on the measurement of the top mass is now more important than the scalar self-energies in the SM.
However, in theories beyond the Standard Model (BSM), it is not possible to make full use of the Higgs-mass measurement because the theoretical uncertainty on the mass calculation can be much larger, owing to the new degrees of freedom. As a result, an enormous amount of effort has gone into refining the calculation of the Higgs-boson mass from a given set of physical or topdown inputs, in both generic and specific theories. This has typically been driven by the need for accurate predictions of the Higgs mass in supersymmetric models, where the Higgs quartic coupling is predicted from the gauge couplings (and other top-down parameters in extended models).
Some results now also exist beyond the gaugeless/effective-potential limit: a complete effective- The purpose of this work is to finally complete the set of scalar self-energy diagrams in the gauge coupling for general renormalisable theories, and provide the tadpole diagrams at the same time. This completes the set that was promised in Ref. [112]. In this paper we present the analytic expressions and the technical machinery that we have used, specialising to the Feynman gauge. However, since the final results (and therefore this paper) are rather long, we have created a new package TLDR where they are available in computer-readable form. Readers wishing to skip the details and apply the results are invited to download the code from: http://tldr.hepforge.org While some of the evaluation of spinor/Lorentz traces and tensor reduction for specific models could be accomplished using TwoCalc [120] and TARCER [121] (of which we have made use) we derived some reduction rules not available there with the help of the general relations of Refs. [122,123], so that all results can be reduced to a basis of just a few one-and two-loop scalar functions that can be numerically evaluated in TSIL [124]. Moreover, our results are already renormalised, and in particular we reduce the number of classes of diagrams by making extensive use of identities relating couplings of ghosts and Goldstone bosons to other couplings in the theory.
Our calculation can be used for: • Corrections to charged and/or coloured scalar masses. For example, in supersymmetric theories this would mean e. g. squark or sgluon masses.
• Electroweak corrections to the Higgs-boson mass ↔ extraction of Higgs/neutral scalar quartic couplings. These ought to be supplemented by a two-loop extraction of the electroweak expectation value and gauge couplings (which requires the two-loop corrections to muon decay and vector-boson self-energies, which we hope to return to in future work).
• EFT matching of the Higgs quartic coupling via the pole-mass matching technique [71]. As described at one-loop order in Ref. [94], although a priori this would seem to require a calculation of the Z-boson mass, in fact all of the necessary information is contained in the scalar self-energies/tadpoles, and thus the calculations here may be sufficient.

Treatment of tadpoles and application of our results
The main application of our results is expected to be the evaluation of the pole mass for scalar bosons in any given theory. For a general theory with scalars having indices i, j, . . . and masses m 2 i at the tree level, this corresponds to finding the (complex) solutions of the equations in p 2 , where Π ij (p 2 ) is the self-energy. We then write Π ij (p 2 ) = Π (1) ij + Π (2) ij + . . .
where the superscripts denote the order of perturbation theory. There are three main techniques used to solve this in practice. The first is to iteratively evaluate Eq. (1) by starting with p 2 = m 2 i ; this does not respect gauge invariance or perturbation order. The second is to perturbatively expand the momentum as p 2 = m 2 i + p 2 1 + p 2 2 + . . . and use matrix perturbation theory to extract the pole mass at each order; this gives a gauge invariant result which respects the order of perturbation theory, but can be tricky to implement in the cases of some masses being degenerate. The third method (see Ref. [112]) is to solve iteratively, then expand the one-loop self-energy, giving and insert into Eq. (1) to solve for p 2 2 , and so on. This method is gauge invariant but does not respect the perturbation order in p 2 , since p 2 1 will contain contributions from Π (1) ij (m 2 i ) 2 and higher powers, etc. In particular, if the aim is to use the pole-mass-matching procedure to extract threshold corrections, then the first or third methods will lead to uncancelled higher-order logarithms. We therefore recommend the use of the second approach, whereby the pole mass at the two-loop order is simply ii m 2 i + Π (2) ii m 2 i + Π (1) ii and where the one-loop self-energies have already been diagonalised on the subspace of states that are degenerate at the tree level; the sum over j includes all scalar and vector states with the same quantum numbers. In general, it is also necessary to include the contributions of tadpole diagrams. This is because, for each neutral scalar with a non-trivial expectation value, there is a non-trivial vacuum minimisation condition, which can be used to eliminate one parameter from the theory. Commonly, parameters of mass dimension 2 are substituted in this step, such as the µ 2 |H| 2 term in the Higgs potential of the SM. If we write the neutral component of the Higgs field as H 0 = 1 √ 2 (v + h + . . . ) and the quartic term as λ |H| 4 , then the necessary condition for the vacuum being a minimum is 0 = µ 2 v + λ v 3 + ∂∆V ∂h (6) where ∆V are the loop corrections to the effective potential; its derivatives correspond to tadpole diagrams. If we insist that v is the correct vacuum expectation value to all loop orders, then we can eliminate µ 2 wherever it appears in favour of λ v 2 and tadpoles. Formally, then we can expand where the superscripts denote loop orders. Since the tree-level Higgs mass is m 2 h = µ 2 + 3 λ v 2 , this means that the tadpole contributions modify the Higgs mass at higher orders.
Utilising this method for the computation of the mass of any particle in the theory that depends on µ 2 requires the calculation of self-energies and tadpole diagrams. In particular, the application to computing the mass of the Goldstone bosons leads to the GBC: in Landau gauge the "tree level" Goldstone mass squared becomes of one-loop order and is of indeterminate sign, which causes infra-red singularities in the two-loop tadpoles and self-energies [125,126]. The initially proposed solution was to resum Goldstone diagrams, and this was performed for the tadpoles of the MSSM in Ref. [127]. However, in general this is cumbersome to implement; and general solutions now exist for both, tadpoles and self-energies, where we can instead use an "on-shell" mass for the Goldstone bosons [114], or just perturbatively expand the generalisations of Eq. (7) to the loop order that we are working to in tadpoles and self-energies [115], known as taking "consistent tadpole equations"; Figure 1: Examples of one-particle-irreducible diagrams with "internal" propagators that we do not include in our list of topologies. All diagrams of these classes can be straightforwardly computed from the one-and two-loop self-energies and tadpoles that we include here.
for example, we would need to take Π (2) ij p 2 → Π (2) ij p 2 In this way, the infra-red singularities cancel between the two parts on the right-hand side, and this should continue order by order in perturbation theory. Another way of treating tadpoles is to work only in terms of running parameters, so our expectation values solve the tree-level vacuum-minimisation conditions only. That means that we must include tadpole diagrams as part of the self energies: these are one-particle irreducible but contain propagators carrying no momentum (referred to as "internal"). This was the approach used in the SM calculation of Ref. [18] and leads to a gauge-invariant result, without needing to perform expansions of the form of Eq. (7), at the expense of a proliferation of diagrams, such as those depicted in Fig. 1.
In this paper, we shall work in the Feynman gauge. Prima facie, one would think that our results do not suffer from infra-red singularities, and so we sidestep the GBC, and could avoid making an expansion of the form of Eq. (7): we could (as originally envisaged in Ref. [112]) just modify µ 2 at the tree level so that Eq. (6) is satisfied and dispense with any extra diagrams or shifts of the form of Eq. (8). A reader wishing to implement this program can use the results from the appendix (121 self-energy classes and 25 tadpoles) or the reduced set of 89 self-energy diagrams described in Sect. 4. However, such an approach does not respect the order of perturbation theory or gauge invariance. † Hence, we shall simplify the expressions where all couplings and masses are expressed in terms of tree-level parameters, so that the reader may use any of the other approaches. In this way, we are also able to make extensive use of relationships amongst the couplings in general theories, and so reduce the number of topologies to 16 for tadpoles and to 58 for self-energies for non-Goldstone boson scalars.
To apply the results in practice for all but the simplest of models is a task for a computer. Even for the Standard Model it would be far too tedious to implement by hand. Hence, we have provided our results in electronic form as part of TLDR. All reduction rules are included so that the renormalised expressions can be applied for any physically relevant configurations of masses and momenta, in the form of Mathematica modules and notebooks, with code to link from Mathematica to TSIL. A user manual is provided online. In the future we also intend to include c++ code to link to TSIL for use with packages such as SARAH, the (currently private) results for Landau and general R ξ gauges, and extensions to vector/mixed scalar-vector/fermion self-energies.

Guide to the paper
The paper is organised in the following way: • in Sect. 2 we introduce our nomenclature for the fields and couplings that appear in renormalisable theories, and we explain our method of computing the two-loop diagrams and counterterms. As our results are valid for real fields, we also show how to apply them to complex scalars and Dirac fermions.
• The purpose of Sect. 3 is to reduce the number of different diagrams by making use of relations among the couplings that are dictated by gauge invariance. In this way, ghosts and Goldstone bosons can be eliminated from the theory.
• The nomenclature that is used for the loop integrals is introduced in Appx. A. There, we also describe how each integral can be reduced for any kinematic configuration into a basis that can be quickly evaluated numerically.
• The full list of results for renormalised two-loop tadpoles and self-energies in terms of the previously defined couplings and loop integrals is given in Appx. B.
• The substitution of ghost and Goldstone couplings is applied to these results in Sect. 4. In the same section, we also compare to previously known expressions.
• Our conclusions are summarised in Sect. 5.

Notation and methods
In this section we shall give our definitions and methods, needed to understand the results presented in Sect. 4 and the appendix.

Coupling definitions
In Ref. [112], scalar self-energies were given using two-component spinors and a compact notation in terms of couplings in a general lagrangian. Such a lagrangian reads where The fields Φ i with indices {i, j, k, l} denote real scalars, ψ I with indices {I, J, K, L} Weyl fermions, A a µ with indices {a, b, c, d} gauge bosons, and ω a , ω a ghosts and antighosts (which carry gauge-boson indices). The Eqs. (10) slightly differ from the ones in Ref. [112] because we work with a different metric signature (+, −, −, −) and, since we work away from the gaugeless limit/Landau gauge, we have couplings between scalars and ghosts.

Scalar, vector and fermion couplings
In this work, due to the large numbers of diagrams and the complexity of the expressions, we perform the generic calculations using computer algebra, and we expect that the application of the results will be best accomplished by implementation on computers. Hence, we use coupling definitions that are more practical for that case; indeed, since we use FeynArts to generate the set of generic diagrams, we adopt an abbreviated form of the notation of the generic model file Lorentz.gen. FeynArts works in four-spinor notation and distinguishes particles and antiparticles; we remove this distinction in our results by transforming all fields to real scalars, gauge bosons and Majorana fermions. Our results are therefore given in terms of vertices which, in general, have more than one possible Lorentz structure; we denote this with an index as last argument of each coupling. Our vertices are named by their adjacent particles, S for scalars, F for fermions, V for vectors, and U for ghosts. The dictionary with Eqs. (10) is: For the terms with fermions, the Lorentz structure with index 1 or 2 refers to a left-or right-chiral projector, respectively; the FFV couplings contain a gamma matrix in addition.
The remaining couplings require more inspection. The SSV vertex is given by We note that we do not enforce the equality of the two terms for the counterterm vertex. Next, for the pure gauge-coupling terms, we have (η µν is the Minkowski metric) For the four-vector coupling, we have the generic vertex with a sum on e.

Ghost couplings and ghost flow
Finally, for the ghost terms, in the FeynArts generic model, no particular form of the couplings is enforced: the general vertex for the UUV coupling is equal to whereas from Eqs. (10) we see that in general R ξ gauge one of these terms is always vanishing because the vertex only contains a ghost and antighost (not ghost-ghost or antighost-antighost) and moreover only contains a factor of the momentum of the antighost. Hence we must preserve the distinction between ghost and antighost in our amplitude, and we should also include both directions of ghost flows. Then we have The signs here are confusing, because we expect that the ghosts are anticommuting, and so should be the vertices. This is an artefact of the algorithm used to construct the amplitudes, where the ordering of the indices in the couplings is not meant to be taken literally: it is assumed that for a given choice of particles in a vertex, either there is only one way of combining them into a vertex, and this is the ordering that is implied (e. g. for a scalar-ghost-antighost coupling there is only one correct choice), or the ordering does not matter. However, there are two good reasons that the reader does not need to worry about this issue: the first is that they can just take the above prescriptions and plug them into our results, given in the appendix. Importantly, for the ghost amplitudes they should sum over both signs of each ghost index, e. g. for diagram I (2) 025 we have the result that should be interpreted as However, the second reason to not worry about this is that in Sect. 3 we demonstrate how all of the ghost couplings can be removed from the amplitude, giving a much smaller number of classes of diagrams to evaluate.

Processing of diagrams
Here we describe our approach to generating and renormalising the diagrams.
Feynman-diagrammatic approach: The two-loop Feynman diagrams are generated with the help of FeynArts [128][129][130]. The different one-particle irreducible topologies of tadpoles and selfenergies are depicted in Fig. 2. Each of these topologies is populated with all possible combinations of fermions F (straight lines), scalars S (dashed lines), vector bosons V (wavy lines), and ghosts U (dotted lines) in renormalisable theories. The external legs are fixed to be scalars. Each particle in the diagram is assigned a unique index ∈ {i1, . . . , i7} for full generality.
BPHZ method: In general, the integrals of the two-loop diagrams in Fig. 2 are ultra-violet (UV) divergent. In order to regularise them, we follow the BPHZ prescription [131][132][133]. Each two-loop diagram contains UV-divergent sub-loops of one-loop order that can be regularised by appropriate counterterms. In general, an additional two-loop counterterm is necessary in order to regularise all UV poles. While the latter can be defined as a pure polynomial in the UV regulator 1/ , the former regularise all non-local divergences and in general give rise to UV-finite terms as well.
The realisation of the splitting of two-loop diagrams into the so-called forest of sub-loop diagrams is carried out in an automated way. For each sub-loop diagram, the corresponding one-loop diagram with counterterm insertion is generated. In addition, the appropriate counterterm topology that regularises the logarithmic divergence of this sub-loop is determined automatically. The algorithm is based on a graph-theoretical interpretation of Feynman diagrams: first, the closed cycles (loops) of each diagram are identified; second, for each cycle, the lines (internal propagators) that make up the loop are shrunken into a point (counterterm vertex); third, the adjacencies to the shrunken cycle in the original diagram (and the cycle itself) are used in order to determine the required counterterm insertion. An illustrating example of this procedure is given in Fig. 3.

Symmetry factors:
Analytical results in terms of amplitudes for the genuine, generic two-loop diagrams as well as the corresponding one-loop diagrams with counterterm vertex, and the oneloop counterterm insertions are generated with the help of FeynArts, FormCalc [134,135] and TwoCalc [120] (we also make use of OneCalc that is part of TwoCalc). Note that among the 121 different self-energy and 25 different tadpole diagrams that FeynArts generates, those diagrams with symmetries of internal propagators have already been removed, whereas diagrams with symmetries with respect to the external propagators are all kept. This processing fixes the symmetry factors of the genuine two-loop diagrams. The symmetry factors of the counterterm diagrams are modified accordingly in order to match the regularised two-loop diagram.     Fig. 2 are marked in red at the left column. After shrinking the lines of these sub-loops, the remaining one-loop diagrams with counterterm vertex are displayed at the middle column. The appropriate counterterm insertion for each counterterm vertex is given at the right column. The consecutive numbers at the lines label the propagators; the same number in the diagrams of one row refers to the same particle at that propagator. A signed number indicates the antiparticle to the unsigned one (note that external particles are defined as incoming by default).
Couplings: All occurring couplings are re-labelled by the sequence of acronyms for all adjacent fields and carry the corresponding propagator indices as argument. An additional numeric argument at the last position allows one to distinguish the couplings of the same fields with different Lorentz structure, as described in Sect. 2.1 where we give the dictionary to parameters in the Lagrangian. In FeynArts the propagator indices may be signed in order to refer to antiparticles. ‡ The different Lorentz structures of the couplings are summarised in appendix A of the FeynArts manual [130]; the argument refers to the -th entry of the vector of couplings. The only deviation from the default setup of the Lorentz structure applies to the coupling SSV[ia, ib, ic] that depends on the momenta k ia and k ib of the scalar fields: at the level of the counterterm vertices not all instances are proportional to (k ia − k ib ); instead, the dependence on k ia and k ib needs to be distinguished. For this purpose, the default generic model file (and the model file) of FeynArts is initialised with the modifications given in Fig. 4; for the couplings (rather than counterterms) we keep the original vertex SSV[ia, ib, ic, 1]. § Counterterm vertices: Since the same counterterm vertex can emerge from shrinking sub-loops of different two-loop diagrams, it is mandatory to store information about the two-loop topology along with the counterterm vertex in order to determine the appropriate insertion. Our choice is to stay close to the existing description in FeynArts, and to assign a canonically ordered list of edges to each Feynman diagram. The edges correspond to the propagators (including the field type) of the diagram. The propagator indices are stored as well since they are required for correctly labelling the counterterm insertions; for the purpose of sorting or identifying diagrams they are not considered. Due to the canonical ordering of the edges, the direction of a propagator might be reversed. In that case, the propagator index of that line receives a sign, indicating an antiparticle. An example of this correspondence is shown in Fig. 5. ‡ Note that at this stage the fields F, S, V, U are not yet specified and may themselves be physical antiparticles.
Therefore, a signed index of an antiparticle refers to the particle. § The generic model file Lorentz.gen is utilised together with the model file SM.mod. The latter is required in order to load all structures for FeynArts. It already contains all possible renormalisable generic couplings.    Counterterm insertions: All one-loop one-point, two-point, three-point and four-point processes were computed at the generic level using FormCalc, since they can in principle appear as counterterm insertions. Each of the Feynman diagrams that appears in these processes was evaluated separately and stored together with its list of edges. In this way, individual results can be looked up quickly and inserted into the correct counterterm vertex; they may in future be provided as part of TLDR. In fact, to evaluate our results we use a separate algorithm to calculate just the divergent parts of the counterterms (which are identical in MS or DR at the one-loop order) on the fly, and which works in any gauge.
Gauge fixing and regularisation: All two-loop self-energy and tadpole diagrams, the corresponding one-loop diagrams with counterterm vertex, and the UV divergences of the counterterm insertions have been determined in 't Hooft-Feynman gauge, Landau gauge, and the general R ξ gauge. However, diagrams with a large number of gauge bosons can initially be expressed in a much shorter form in 't Hooft-Feynman gauge. For this reason, we performed the complete integral reduction only for the Feynman gauge (see Appx. A). In addition, we computed our results using dimensional regularisation (for MS renormalisation), and dimensional reduction (for DR renormalisation). The extra terms in dimensional regularisation have a coefficient δ MS .

Results:
The outcome for the combination of each genuine two-loop integral with the corresponding counterterms for all sub-loops is given in Appx. B.1 for the tadpoles, and in Appx. B.2 for the self-energies; first, all previously known expressions of Refs. [112,113] are given in our nomenclature, and then all new results are stated. Note that the pure UV divergent two-loop counterterms are not contained in these results, as they simply add clutter to the expressions. After carrying out the integral reduction and extracting all UV divergences via the relations in Appx. A, polynomials in the regulator 1/ will remain, which simply correspond to the genuine two-loop counterterm of the diagram; but we may also find additional divergent and finite parts corresponding to any infra-red divergences: we do not introduce infra-red counterterms. As described in Sect. 1.1 any infra-red divergences should cancel in a full amplitude when combined with one-loop self-energy and tadpole diagrams.

Basis of loop integrals:
We have provided all of the integral reduction rules in Appx. A to extract the finite part of our generic renormalised amplitudes for any kinematic configuration in terms of the basis of one-and two-loop scalar functions that can be evaluated in TSIL. This is a basis of six two-loop scalar integrals and two one-loop ones. The TSIL basis are actually "renormalised" integrals, with specific subtractions; essentially they can be built up from Ref. [136]: where the definitions of the boldtype integrals are given in Appx. A.3.1, and the number of spacetime dimensions is d = 4 − 2 . In our results in the appendix we use the equivalent of the boldtype integrals, and so when taking the finite part of the diagram we must use where Fin[. . .] denotes the finite part as → 0, i. e. neglecting non-zero powers of , and In principle, since B(x, y) is known analytically to all orders in , it is not difficult to evaluate B | 0 (and it is even available in TSIL). However, its presence is a sign of an infra-red divergence: in the absence of infra-red divergences, all B | 0 integrals cancel in the renormalised amplitude. We have explicitly checked that this is the case for all the results in the appendix. So then the reader might be curious as to why we do not give the results in terms of the renormalised TSIL basis, as done in Ref. [112]; the reasons are twofold: • Diagrams with vector bosons contain many more tensor integrals. By listing them in unreduced form we are able to give our results for each diagram in only a few lines, whereas some diagrams could fill pages by themselves in expanded form, even just for the most general case where there are no vanishing or degenerate masses.
• The integral reduction depends crucially on whether there are coincident or vanishing masses. The simplest cases concern scalar integrals with repeated propagators, where for non-identical masses we can use partial fractions to reduce them (see Eq. (121)). In Ref. [112] there were generally one or two such cases per integral. In the implementation of the results of Refs. [114,115] in SARAH, one scalar integral contains 12 different reductions. However, in our results, we also have cases in the reduction where some vanishing gauge-boson or ghost masses naively lead to poles; furthermore, there are also special cases with the external momentum taken equal to the scalar mass that are relevant (in particular) for charged scalar propagators. Hence our results (for example for I (2) 086 ) have up to 47 different special kinematic configurations! In TLDR we provide the different special cases for each diagram, and our reduction rules transform them to the appropriate renormalised integral; clearly it would be impractical (and not especially useful) to list all of the final results here.

Charged scalars, Dirac fermions and γ 5
Our results are calculated using Majorana four-spinors, and in Sect. 2.1 we gave a prescription of how to translate them to Weyl-spinor notation-but under the condition that the spinors are also Majorana, i. e. have diagonal (and real) masses. In Ref. [112] the results were presented in Weylspinor notation without this condition, by allowing Dirac spinors to have off-diagonal masses M IJ where M IK M KJ = m 2 I δ IJ ; such mass terms simply link the left and right Weyl spinors that form a Dirac spinor. To translate our results to that notation is surprisingly easy, because we pull out linear factors in the mass for each diagram including fermions. For example, we can write for I For practical applications, however, it is likely most useful to retain the four-spinor notation; a translation to Dirac spinors as a sum of two Majorana spinors, and to complex scalars as a sum of two real scalars is straightforward. Naively, for each complex scalar or Dirac spinor the number of amplitudes would be doubled in this way, but actually there must be a symmetry preventing the two components of a complex scalar or a Dirac spinor from mixing or splitting their masses-hence for any three-point coupling there is a unique way of combining them into a coupling. For example, when considering three complex scalars φ 1 , φ 2 , φ 3 where the lagrangian contains or their complex conjugates are not allowed because they would violate the symmetries keeping the fields complex. Therefore, the notation with Majorana spinors and real scalars can be very easily applied to Dirac spinors and complex scalars: when inserting fields into our results we just choose for each set of couplings the one combination of fields that is allowed in the theory. Indeed, this is the algorithm used in SARAH (albeit currently for the diagrams in the gaugeless limit only).
Finally, another motivation for retaining the Majorana-spinor notation is the problem of γ 5 . For two-loop self-energies, we can use a naive anticommuting prescription for γ 5 , but for two-loop decays or three-loop self-energies this is known to be inconsistent and it would be necessary to choose another definition, for example involving the µνρκ tensor. In the two-component formalism, spinors are automatically split by chirality, corresponding to a naive γ 5 prescription, and it is not known how to make this consistent for such higher-order calculations.

Removing Goldstone and ghost couplings
In this section, we shall derive tree-level relations among the Goldstone and ghost couplings of general theories, which extend those in Ref. [137,138], and eliminate four-point couplings involving vectors. Indeed, the first such relation is already implicitly included in Eq. (10); the four-vector coupling is just related to the product of three-vector couplings, as can be seen from either the requirement of unitarity of the amplitude, or just read off from the standard kinetic term. We shall derive all the necessary couplings using the second approach; starting with fields in the gauge eigenstate basis of scalars R i and gauge bosons V a Once we break the gauge symmetries, we must diagonalise the vector masses. Naively this just involves orthogonal rotations and therefore the sum over intermediate states is not affected; however, in the presence of kinetic mixing of U (1) gauge bosons we must first make a non-orthogonal transformation. We must first unravel the kinetic mixing via but then, since only U (1) gauge bosons may mix, we must have f abc Z cd = f abd . We subsequently make an orthogonal transformation to diagonalise the vector masses; and the relation between quartic and cubic gauge couplings is clear; also that g abc is antisymmetric. For four-point scalar-vector interactions, we must look at the covariant derivative of the scalars in the gauge-eigenstate basis: where θ a ij are real antisymmetric matrices. They obey the group algebra [θ a , θ b ] = −f abc θ c . This then yields The scalars are rotated by an orthogonal transformation (there is no kinetic mixing at tree level) so we have lj , g abij = g aki g bkj + g akj g bki .
The couplings g aij are antisymmetric on the exchange of the two scalars. It should be noted that the assumption that the scalar rotation is orthogonal will be violated if the running parameters do not sit at the minimum of the tree-level potential, e.g. if we work with running parameters that sit at the minimum of the full loop-corrected potential. Such an approach, however, leads to many complications in the calculations and we do not recommend it. Alternatively a choice of finite counterterms (such as using an on-shell scheme) may cause the identity above and in the following to be violated.

Ghosts and gauge fixing
To derive the Goldstone and ghost couplings, we first need the gauge-fixing terms. Once we give expectation values to the scalars, so R i = v i +R i , and defining the scalar kinetic terms contain We thus have the R ξ gauge-fixing terms defined so as to remove tree-level kinetic mixing between scalars and vectors.
Rotating the gauge transformations in the original gauge basis α a so that α a ≡ Z ab α b , we have This gives From this we can read off the ghost mass matrix and the ghost couplings. For ghost-vector couplings, we have Also as expected, is the mass matrix for the gauge bosons too, and so we can diagonalise both ghosts and vectors with the same orthogonal rotation O ab , and for each massive vector of mass m a there will be a ghost with mass √ ξ m a . However, it is important that, since the ghosts and antighosts are not identical, we can treat them as complex fields, and we actually have the liberty of defining them with an additional phase.
From Eq. (30), after diagonalising the scalars we read off that, as noted in Ref. [139], but sinceĝ abi has an antisymmetric piece we cannot simply invert this relation. In order to writeĝ abi in terms of the other couplings of the theory we will have to consider the Goldstone bosons.

Goldstone bosons
The gauge-fixing terms are also expected to give mass to the Goldstone bosons (except in Landau gauge); they contain scalar mass terms To see that these concern just the Goldstone bosons, consider the standard perturbative proof of Goldstone's theorem: we expand the potential without gauge-fixing terms and differentiate the relation once: If we work at the minimum of the potential, this becomes This is true for any α a ; the F a i are null eigenvectors of the mass matrix until we add the gauge-fixing terms. Adding the gauge fixing terms we have Now let us use the singular value decomposition of F a i and define, suggestively: ji are orthogonal and F D is a diagonal-but not in general square-matrix. Since ji is arbitrary when acting on Goldstone-boson indices; it can be chosen to simultaneously diagonalise both matrices in Eq. (40), splitting the scalars into would-be Goldstone bosons and a remaining set. Then, in the diagonal basis Φ i , F D becomes a projector onto Goldstone bosons: where N G is the number of Goldstone bosons/massive vectors. Armed with this, we can writê i. e. we can exchange the scalar-ghost coupling for scalar-vector couplings. However, because the coupling has a different form depending on whether the scalar is a Goldstone boson or not, this introduces some complications in calculating amplitudes: we should either sum separately over Goldstone indices and the remaining scalars (which is necessary for finding a gauge-invariant result), or just use these pieces to remove the ghost couplings. In this work we take the latter approach.
Another point is in order: there is actually some ambiguity in the above definitions because we have the freedom to introduce signs/phases of the Goldstone bosons and ghosts. In all calculations such signs should drop out. This is in particular notable because implementations of gauge-fixing in the literature are defined via the standard procedure but, in order to verify the above relation (and, indeed, those that we shall introduce below), it is necessary to introduce such signs/phases-we checked our identities against the FeynArts model file of the SM, for example, where we need to introduce a sign for the Z-boson Goldstone and a factor of ı for the W -boson Goldstone.

Eliminating all Goldstone-boson couplings
Now that we have given explicit forms for the ghost couplings in terms of scalar and vector couplings, and found that the form depends on whether the scalar is a Goldstone boson or not, we can also consider relating the couplings of Goldstone bosons to other couplings in the theory. Partial results for these can be found in Ref. [137,138]. The general strategy will be to use our projector F D and invert the relation in Eq. (41), and use the identity Throughout we shall distinguish would-be Goldstone bosons from ordinary scalars S by using the letter G to represent them, and we use G a , G b , G c , . . . for their indices (instead of i, j, k, . . .): the subscript a, b, c, . . . of course indicates that they correspond to the gauge boson of that index.

Scalar-vector couplings
GGV Inserting our projector into g aij we have: Next we use and therefore the coupling of two Goldstone bosons to a gauge boson is: This is the expression found in Ref. [138] by requiring that high-energy scattering amplitudes of theories with massive gauge bosons have the correct behaviour.
SGV For a general scalar, coupled to a Goldstone boson and a vector, we can derive GVV Consistent with the previous two expressions, we can derive

Scalar-Goldstone couplings
The pure-scalar interactions involving Goldstone bosons can also be related to couplings involving vectors, thus allowing them to be eliminated. To derive the required relations, we continue to apply derivatives to Eq. (38): SSS From the first line of Eq. (51) we get the GSS coupling where we have denoted by m 2 0,i the eigenvalues of ( These are equal to the scalar masses in the full potential, except that they do not include contributions from gauge fixing: wouldbe Goldstone bosons have m 2 0,Ga = 0. It is also useful to have the expression for an GGS coupling: where we note (as shown in Ref. [114]) that a triple-Goldstone coupling vanishes. To avoid using m 2 0,i we can write and therefore SSSS From the second line of Eq. (51) we retrieve the GSSS coupling Here there is a sum over all scalars i including Goldstone bosons, and indeed j, k, l can be Goldstone fields. To eliminate couplings with more Goldstone bosons we then just need to insert the formulae that we have given above into this equation; for example (as we shall later require) a four-point coupling of two Goldstone bosons and two scalars that are not Goldstone bosonsk,l is Also particularly interesting (but not needed for this work) is the four-Goldstone coupling, which can be written as (summing over all non-Goldstone scalars i)

Fermion couplings
To complete the removal of Goldstone boson couplings, we would require the couplings of fermions to Goldstone bosons. We will not actually use these in this work, and they were given in Ref. [137]; we give them again here in our notation for future reference:

Results
In this section we shall describe our results. The renormalised expressions for all of the basic classes of diagrams are given in Appx. B.1 for the tadpoles and Appx. B.2 for the self-energies, since they are rather long; initially, there are 25 tadpole and 121 self-energy classes which is a much larger set than in the gaugeless limit. Therefore we shall describe how we can reduce this set to 89 or 92 for generic self-energies (depending on whether we choose to exchange the ghost-ghost-vector coupling for a triple-gauge coupling) or 58 for non-Goldstone scalars; and just 16 for the tadpoles.

Tadpole diagrams 4.1.1. Unreduced diagrams
To present our results in a readable way, and to make the connection with the diagrams in Ref. [113], we shall denote tadpole topologies 1, 3 and 2 with all scalar propagators as T SS , T SSS , T SSSS ; the subscripts are modified for different fields accordingly. The total tadpole can be written as where the superscript (2, n) indicates O(2n) in the gauge couplings. Then the unreduced set of tadpole topologies is: The integrals I N can be found in Appx. B.1. Of these, the expressions in T (2,0) are equivalent to equations (2.32), (2.33), (2.34), (2.36) and (2.37) in Ref. [113], since they are independent of the gauge fixing (and the results in the gaugeless limit were given there).

Combined diagrams
The diagrams with fermions are irreducible, but we can exchange all four-point couplings including vectors, and all ghost couplings, for three-point couplings involving vectors and scalars using the identities in Sect. 3. This means that we can reduce the number of topologies by combining the loop functions together. To do this, we need some notation: we write each integral in the appendix as a sum over the different combinations of Lorentz structures multiplied by a loop function. Suppose we have an integral with n propagators and m couplings with p indices in total; then let us write the couplings generically as where {L 1 , . . . , L m } denote the Lorentz structure of the couplings. The diagrams can be written as In the cases where there is only one function, we omit the superscript with the Lorentz indices. For example, we can write but For fermions, there are several Lorentz structures and the loop functions differ, while for amplitudes without fermions there are at most three, and the loop functions are typically equal or differ very little (in this example, t ). Armed with this notation, we can then combine the amplitudes. For all the combinations, we will also be able to reduce each class to just a single Lorentz structure. It is therefore straightforward to convert our expressions in FeynArts-based coupling notation to the notation of Eqs. (10) using the identities in Sect. 2.1. The result for the tadpoles is: We see that the original set of 25 topologies is reduced to just 16. A more detailed description of the subtleties in the reduction involving scalar-ghost couplings is given in the next section. Here we simply present the results for the reduced expressions in turn.
The simpler combinations are while the more complicated combination is 4.2. Self-energies

Unreduced diagrams
To make the connection with the diagrams up to second order in the gauge coupling in Ref. [112] and Ref. [113], we write the total self-energy as where the superscript (2, n) indicates O(2n) in the gauge couplings. The minus sign reflects the difference in our conventions compared to FeynArts. The symmetrisation on the external legs is to account for the fact that some diagrams are not explicitly symmetric (for example, we eliminate I since it is identical to I 077 with i1 ↔ i2, i4 ↔ i7) and it is far faster to symmetrise the total final result rather than evaluate twice as many diagrams.
The expressions for Π (2,0) and Π (2,1) are known, and we can write them as All of these simple terms except Π (R)SV,g 2 are irreducible.
The diagrams without vectors are while the diagrams with only one vector propagator are In Ref. [112], the expressions are given in unreduced form; those in Π (R)SV,g 2 can be simplified by removing the four-point scalar-vector interaction, as we shall do in the next section.
Next we can write The letters I, R represent irreducible and reducible classes respectively; the final two contain special classes. The expressions are For the third-order terms, we have Π (2,3) i1,i2 = Π (I),g 6 + Π (R)SV,g 6 + Π ‡,g 6 .
The expressions read Finally, the fourth-order terms are given by Together, these describe 92 classes. However, we can immediately reduce this to 89 by removing the ghost-ghost-vector couplings inside Π ‡,g 4 and Π ‡,g 6 , as we describe in the next sections.

Combination of classes
As in the case of the tadpole diagrams, we shall exchange four-point couplings involving vectors, and all ghost couplings, for three-point couplings involving scalars and vectors. We therefore use the same notation as in the tadpole case to represent our loop functions and their couplings; the self-energy diagrams can be written as In the cases where there is only one function, we omit the superscript with the Lorentz indices. For example, we can write but The first trivial application of this is to remove the diagrams with ghosts that do not couple to scalars. In the above notation we can write where we note that f , and Π ‡,g 4 consists of two classes corresponding to diagrams I (2) 099 and I (2) 119 : This reduces the number of classes of diagrams to evaluate to 89, which will speed up evaluation (since the loop functions are all of the same class, they require no substantial extra time to evaluate, whereas performing loops over the couplings and evaluating the loop functions each time is slow). In fact, by using Eq. (28) we can combine these two classes of diagrams into one graph of topology V SV V V V ; but we shall apply this systematically in the next section.

Ghostbusting
In the previous section we eliminated diagrams containing only ghost-ghost-vector couplings and no ghost-ghost-scalar couplings. In this section we shall remove all ghost couplings from the amplitude, and also apply Eq. (28) and Eq. (14d) in order to obtain only 58 classes to evaluate in total (compared to the 28 for up to O(g 2 ) terms) at the expense of requiring that the external scalars are not would-be Goldstone bosons.
The key equation that we shall apply is Eq. (44) which shows how we can relate SUU couplings to scalar and vector couplings. However, the form of those couplings has an extra contribution for would-be Goldstone bosons: We therefore have to make some distinction between would-be Goldstone bosons and the other scalars in the summation. One approach to dealing with this would be to introduce the Goldstone bosons as a new class of fields (separate from other scalars) from the start, and remove them via the identities in Sect. 3 afterwards. This would be necessary to obtain an explicitly gauge-invariant result, and should lead to a faster evaluation of the final result, but comes at the expense of complicating the expressions. We leave this to future work. Instead, we deal with the above problem by noting that the first term on the right-hand side of Eq. (82) is universal for all scalars whether they are Goldstone bosons or not, so we can split any diagram containing an SUU vertex into two (or more) and explicitly sum over the second part, since it effectively becomes a vector propagator. As an example, consider I (2) 027 : Two combinations here have dropped out: Schematically we then have and we see that the diagrams that drop out would not have fit with the above picture (note that almost exactly the same pattern reproduces for diagram I 061 ). In the first diagram on the left, the sum over scalar propagators is indeed a sum over all scalars in the theory. We restrict to the case that the external states are not Goldstone bosons here, because otherwise the couplings for the diagrams on the right hand side of the above relations would include gauge bosons as external legs, and we would therefore need to combine those amplitudes with mixed scalar-vector and vector-vector amplitudes.
There are also complications when the "Goldstone" leg attaches to a triple or quartic scalar vertex. The reason for this can be traced back to Eq. (52): masses appear in the GSS coupling relation that are different for Goldstone bosons than for other scalars (in any gauge except Landau gauge), and this then feeds into the relation for the GSSS coupling. There would also potentially be a similar problem for the SGV vertex, but in all our examples this is avoided because we assert that we have no external Goldstone bosons. There is only one diagram with a quartic scalar coupling and ghosts-I (2) 112 -that we will discuss in more detail below, but diagrams I (2) 025 , I (2) 060 and I (2) 071 contain triple scalar couplings, that we treat by using Eq. (55) and again splitting the diagram into sums over all scalars and sums over just Goldstone indices, explicitly trading the Goldstone bosons for vectors.

Final simplified results for self-energies
The final result for our set of 58 self-energy topologies can be expressed as: The pieces without bars are unchanged from above, the ones with bars are explained in the following.

Reduced diagrams of O(g 2 )
The combined diagrams start at O(g 2 ): The only new diagrams with fermions enter at order O(g 4 ) and are given by: The other combined diagrams at order O(g 4 ) are: These new topologies are: The 10 combined diagrams involving scalar and vector couplings of O(g 6 ) are given by The topologies of type "M " are: The last term in the above appears fearsome, but actually its finite part simplifies dramatically: This expression is untroubled if any masses are coincident or vanishing. The topologies of type "V " are: For the last expression, the factors of m 2 i4 and m 2 i7 in the denominators may at first appear to lead to divergences if the corresponding gauge bosons are massless. However, it is straightforward to show that, in that case and for i1, i2 not would-be Goldstone bosons, the SVV[i1, i3, i4, 1], SVV[i2, i3, i7, 1] vertices also vanish. This expression simplifies to a relatively short form similarly to M V V V V V , except that the reduction is different for m 2 i4 = m 2 i7 as compared to the non-degenerate case. The nominally highest-order diagrams in the gauge coupling consist of just two topologies, reduced from the original four: The expressions for these are: 4.2.9. Special topologies: Π §,g 4 Recall that Π §,g 4 = I 107 + I 112 . These two diagrams contain reducible couplings; the first has a four-point scalar-vector coupling, and the second contains ghosts. However, when we remove the four-point vector/scalar interactions and the ghosts we find topologies that are either not 1PI or contain internal propagators. The first of these is I 107 : The first diagram contributes to topology M V SSV S or I 041 , while the second is a non-1PI diagram; schematically: The second interesting topology is I 112 : Here we have used Eq. (54) and substituted Eq. (55) into Eq. (56). The topologies on the right-hand side correspond to I 116 (W SSV V ) and an extra diagram, schematically: There is therefore some ambiguity in how to treat these two special classes. In the first case, it is probably simpler to retain a definition for the four-point scalar-vector interactions and work with I 107 rather than reducing it. For the second case, the simpler approach depends on the treatment of tadpole diagrams: if we use a standard approach and do not include "internal" propagators, then it may be easier to retain ghost couplings/propagators for just this one class and work with I 112 . On the other hand, if we treat tadpoles by including internal propagators, then all of the topologies on the right hand side of the above diagram already exist as sets of couplings that we should evaluate. In this case it must be simpler/faster to work with the right-hand side of Eq. (97).

Comparison with Landau gauge
To make the connection with the results of Refs. [113,114] completely explicit, we shall present here the expressions for the tadpoles in the Landau gauge and make the connection with our diagrams in the appendix.
The set of diagrams that survive in the Landau gauge are where the superscript (2, n) indicates O(2n) in the gauge couplings. Then we have where The topologies (and therefore labelling of the couplings) are the same in Landau gauge as for our Feynman-gauge results in the appendix, but we add the superscript ξ = 0 here. We do not give the explicit expressions for these, since several are rather long-they contain many more integrals than the Feynman-gauge case-but they are provided as part of TLDR.
Ignoring for the moment cases where some masses vanish, and using the loop functions defined in Ref. [110] and the notation we find the following expressions for the tadpoles: Tadpole topology Label Finite Part Our function The fermionic diagrams are: In Refs. [113,114] the topology T ξ=0 SV V was missing. Our new approach, however, pays immediate dividends when we consider some masses to be vanishing. We find, for example, that where R SS (x, y) is defined in Refs. [114,127]; using the expressions in the appendix we can find the finite part of all of the integrals even in the case of vanishing masses. The second term on the right-hand side of Eq. (103b) clearly arises from an infra-red divergence; such divergences are the counterparts of the expansions found in Ref. [114], where the equivalent for Eq. (103b) is with the renormalisation scale Q. Instead of an expansion in the mass of the Goldstone bosons, we use dimensional regularisation in order to regularise the infra-red divergences. Then we should find that the infra-red divergent parts in the two-loop calculation cancel exactly against equivalent parts coming either from putting the Goldstone-boson masses on-shell in the one-loop parts, or using "consistent tadpole solutions" (by using tree-level masses in the one-loop calculation). We intend to elaborate further on this connection elsewhere.

Standard Model calculation
Our self-energies up to second order in the gauge coupling have been compared with the expressions in Ref. [112], and, as we have shown above, the tadpoles in the Landau gauge exactly match those found in Ref. [113]. In order to explicitly cross-check our Feynman-gauge result, and the combinations of classes, it would be desirable to have models to compare to. However, as described in the introduction, the only example where scalar electroweak corrections have been computed is for the Higgs mass in the Standard Model. Since the results are very long, these have been made public in the form of computer codes SMH [17], mr [18,19] and SMDR [20]. The result used in mr has been computed including tadpoles via internal propagators in a general gauge, and is excellent for a numerical calculation. However, the expressions in SMH are actually written in a sufficiently convenient form to be translated into code readable by Mathematica, and so we were able to use it for an analytical cross-check of our results. On the other hand, the result there is in Landau gauge, and combines the tadpoles with self-energy contributions via Eq. (7). Therefore we have cross-checked our results against the Standard Model by the following procedure: • We calculated the tadpoles and Higgs self-energies in the Landau gauge using the FeynArts model file for the SM (with some modifications that will be described elsewhere), and using TwoCalc, TARCER and our own reduction to reduce the basis of integrals to those that can be evaluated by TSIL. We compared the analytic expressions with those in SMH and found perfect agreement in all terms. For this calculation we computed self-energies and counterterms separately, rather than using the BPHZ method (in the SM all necessary counterterms can be reduced to just evaluating two-point functions since there are sufficiently many tree-level relationships among the parameters).
• Using the same method, we performed the same calculation in Feynman gauge. The results differ from SMH only by the treatment of tadpole diagrams once we take p 2 equal to the treelevel Higgs mass. This leads to 899 self-energy diagrams based on the initial 121 generic classes in FeynArts, and 75 one-loop self-energy diagrams with one-loop counterterm insertions (plus 161 tadpoles and 25 one-loop tadpoles with one-loop counterterm insertions).
• We performed the same calculation using only the topologies of our 58 combined classes with the combined loop functions as given in this section. This yielded only 425 self-energy diagrams (and, of course, includes the counterterm contributions). Upon reducing to the TSIL basis we find exact agreement of all terms with the above "brute force" approach.

Conclusions
We have given expressions for scalar self-energies and tadpoles at the complete two-loop order in Feynman gauge. With the aim of being flexible, we have given the results for all renormalised component diagrams, and also a much shorter version where the diagrams are combined into just 16 tadpole and 58 self-energy topologies, provided that the external scalars are not (would-be) Goldstone bosons. The results are provided in the most compact analytical form that we could give, but are also available electronically as a new package TLDR at http://tldr.hepforge.org (where more tools and calculations will be added over time) so that they can be easily applied. We also include online all of the replacement rules and degeneracies for each diagram with routines for extracting the finite part. This work closes a gap in the literature that has existed for at least sixteen years. However, thanks to the technology and techniques that we have developed, it is just the first step in a new program of completing the electroweak corrections in generic theories. The logical next step is an implementation in a general code such as SARAH so that they can be applied to any model automatically. The following steps are: 1. Gauge-boson and mixed gauge-scalar self-energies. In particular, the latter are required for decay processes, while the former are essential for the extraction of electroweak parameters.
2. Threshold corrections to match general theories onto the Standard Model via a pole-matching technique. Naively this requires matching gauge-boson masses, but as shown in Ref. [94] we only need the derivatives of the scalar self-energies, and gauge threshold corrections at zero Higgs expectation value.
3. Muon decay in general theories. This is required for extraction of the Higgs electroweak expectation value in fixed-order calculations. It is a four-point fermion amplitude, so is combinatorially much more complicated than two-point functions, but since it is at zero external momentum we do not have any new problems related to γ 5 , and the loop functions all reduce to the same ones used in the tadpoles.

Fermion self-energies, in particular to calculate the top-quark mass.
We have applied our results to the Standard Model and compared to the existing expressions in the literature, finding perfect agreement. However, since the only analytically available calculation was performed in Landau gauge, it was not possible to compare our Feynman-gauge calculation with an independent check. However, by combining the diagrams with all tadpole diagrams it should be possible to give an explicitly gauge-independent result that could also be compared with mr. Moreover, using the results from Sect. 3 it should also be possible to implicitly sum over all Goldstone-boson propagators in our general result, transforming their couplings and masses into those of gauge bosons. We leave these developments to future work.

A. Integral relations
In this appendix, we introduce the notation for the integrals that are used in the list of results for tadpoles and self-energies in Appx. B. After mentioning some general relations among these integrals, the reduction rules that are necessary in this article are displayed. Finally, the UV-divergent parts of the basic set of scalar integrals are given.

A.1. Notation and symmetries
We mostly follow the notation of Ref. [120], with the indices i 1 , . . . , i n , j 1 , . . . , j o ∈ {1, 2, 3, 4, 5}, the dimension d = 4 − 2 with the ultra-violet regulator , the dimensional regulator µ, and the kinematic variables that depend on the loop momenta q 1 , q 2 and the external momentum k 6 ≡ p. In addition, we use sub-indices a and b in order to distinguish denominators of the same kinematic type, but with possibly different masses. The set of different two-point integrals can be reduced by using shift and mirror symmetries of the loop momenta (see Ref. [120]). They can be applied by exchanging the indices of the kinematic variables in the integrals in the following ways, 1 ↔ 2 and 4 ↔ 5 , 1 ↔ 4 and 2 ↔ 5 , 1 ↔ 5 and 2 ↔ 4 , We also introduce the following notation for one-loop integrals appearing in counterterms: with i 1 , . . . , i n , j 1 , . . . , j o ∈ {1, 2}. For scalar one-loop integrals with at most two propagators, we adapt the well-known notation of Passarino and Veltman [140]: The reduction of one-loop integrals with more than two propagators will be addressed together with the two-loop integral reduction in Appx. A.3. For convenience, we also introduce the following symbols in order to denote a sum of one-or two-loop tensor integrals: with integer coefficients c 1 , . . . , c r of the products of kinematic variables in the numerators with indices j 11 , . . . , j ro .

A.2. Reduction of tensor integrals
Tensor integrals that carry identical indices in the super-and subscript can be reduced easily via After the repeated application of Eq. (111) to a given integral, the following cases appear: two loop integrals whose integrand does not depend on both loop momenta are equal to zero; the reduction of scalar integrals with a numerator equal to 1 is described in Appx. A.3; the reduction of the remaining tensor integrals is described e. g. in Refs. [120,141]. Reduction rules for the integrals appearing in this article that are valid for all different kinematical configurations (with the exception of vanishing external momentum p) are given in the following: When applying Eq. (111) to Eq. (112i) also the following integrals appear after taking into account the symmetry relations in Eqs. (107):

A.3. Reduction of scalar integrals A.3.1. General relations
Our aim is to reduce all scalar two-loop integrals into the basis of integrals that can be evaluated numerically by the code TSIL [124]. We begin by reducing to the unrenormalised basis where we include the one-loop integrals A and B that are known analytically. While TSIL is capable of delivering these bare integrals, it finds them by solving differential equations for "renormalised" versions; and indeed it is these renormalised integrals that we actually need for practical applications. This finite basis is [136] S(x, y, z) = lim I(x, y, z) ≡ S(x, y, z) The integral T 134 is also known analytically, and is therefore available in other private codes (such as SARAH). Closed analytic expressions for some integrals with special kinematic configurations can be found e. g. in Refs. [142,143].
while the one-loop integrals are ∦ and we define where Q is the renormalisation scale.
In this appendix we shall give the expressions for diagrams in terms of the bare integrals. As described in the text, this is because the reduction is different for different configurations of masses (coincident masses, vanishing masses, special kinematic configurations, certain on-shell conditions).
After applying the symmetry relations of Eqs. (107), the following different scalar integrals remain: The integrals in Eq. (119a) differ from the others since no index appears repeatedly. From those, the integrals that are not part of Eqs. (114) can be decomposed into a product of two one-loop integrals, The integral T 1a1 b 2 is equal to zero in each kinematical configuration; the remaining integrals with the repeated indices 1 a and 1 b can be decomposed via partial fractioning if the corresponding masses m 1a and m 1 b are not equal to each other (the same rule applies to the one-loop integrals of Eq. (108)), Instead, if the masses are equal, a dedicated reduction of the integral is required. The same rule applies to the massless propagators with index 1 . In addition to these obvious qdegeneracies, special kinematical configurations for the other masses and/or the external momentum need to be taken into account in order to avoid manifest poles in the reduction formulas. We do not take into accout accidental thresholds when one mass is equal to the sum of two others (unless one of these masses is equal to zero).

A.3.2. Reduction rules
We have derived the reduction formulas for all integrals that appear in 't Hooft-Feynman gauge with the help of TwoCalc, TARCER [121] and the recurrence relations in Ref. [123]. Instead of writing out the full expression for each integral with multiple index 1 , we present the respective results in terms of recurrence relations (see Ref. [122]) in the following. In some cases closed forms are given. ¶ We make use of the shorthands The reduction of all scalar one-loop integrals is given by ∩ (124e) ¶ Reduction rules for integrals with higher powers of massive propagators can be computed by evaluating the derivative with respect to the mass at this propagator while keeping all other masses as independent variables. ∩ The extension of these rules to two-loop integrals that factorise into a product of one-loop integrals is straightforward.
The required reduction rules for the scalar two-loop integrals follow from T 1···1 n 3 and T 1 ···1 n 3 : T 1123 and T 1 ···1 n 23 : T 1···1 n 34 and T 1 ···1 n 34 : T 11234 : : Additional integrals are introduced by the previous reductions. They can be further reduced by the following relations: : : The integrals with multiple index 4 can be reduced after using the symmetry relation 3 ↔ 4. All remaining integrals belong to the set of Eq. (114) after applying the symmetry relations of Eqs. (107).

A.4. Vanishing external momentum
In the limit of vanishing external momentum, p → 0, many of the reduction rules cannot be applied directly due to manifest poles in p 2 . However, if this condition is imposed onto the integrals in the first place, their reduction becomes much easier. Symbolically, the implementation of this limit amounts to performing the index substitutions 2 → 1 and 5 → 4 (the masses are relabeled accordingly, but in order to distinguish the possibly different masses of repeated indices, new subscripts are introduced). For the tensor integrals, only the first reduction in Eq. (112f) requires some special care, since it contains an explicit pole in p 2 . The momentum-free equation reads which corresponds to the first term on the right-hand side of Eq. (112f) with p → 0. For all other tensor reductions in Eqs. (112) the limit p → 0 can be taken explicitly.
In addition to the types of integrals in Eqs. (119), the following scalar integrals appear after taking into account the symmetry relations of Eq. (107) and partial fractioning via Eq. (121): Further integrals with multiple massless propagators with index 1 can be solved with the recurrence relations given above. In general, the solutions for the integrals of Eq. : T 1···1 : T |m 2

A.5. UV divergences
The remaining scalar integrals are UV divergent in general. In dimensional regularisation, these UV divergences can be parametrised by the regulator . For the two-loop two-point integrals appearing in this article, the terms that diverge in the limit → 0 are known analytically; the finite terms are evaluated numerically with the help of TSIL. In the following we list these UV divergences explicitly with the considered order in indicated in the superscript of the integrals: T 234 : T 1234 : T 2234 and T 2 2 34 : T 11234 :

B.2.2. Known results with scalars and fermions
Known results with scalars and one vector

B.2.4. Known results with fermions and one vector
B.2.5. Known results with scalars, fermions and one vector B.2.6. New results with vectors