Two-loop Higgs mass calculation from a diagrammatic approach

We calculate the corrections to the Higgs mass in general theories restricted to the case of massless gauge bosons (the gaugeless limit). We present analytic expressions for the two-loop tadpole diagrams, and corresponding expressions for the zero-momentum limit of the Higgs self energies, equivalent to the second derivative of the two-loop effective potential. We describe the implementation in SARAH, which allows an efficient, accurate and rapid evaluation for generic theories. In the appendix, we provide the expressions for tadpole diagrams in the case of massive gauge bosons.


I. INTRODUCTION
With the discovery of the Higgs boson in the range of 125 -126 GeV the standard model (SM) has been completed [1,2]. The uncertainty in the Higgs mass measurement has continuously decreased and is well below 0.5 GeV today [3]. This small uncertainty is currently much better than the theoretical prediction in any scenario beyond the Standard Model (BSM). Therefore, more precise calculations are necessary to better confront BSM models with the Higgs mass measurement. This has two motivations of particular weight: (1) In the minimal supersymmetric standard model (MSSM) and many extensions thereof, radiative corrections are required to be at least as significant as the tree-level contribution, so higher-order corrections are especially important. (2) In the Standard Model and non-supersymmetric extensions thereof a precise calculation is required to extract the parameters of the model, which when run to high energies gives information about the stability or lifetime of the potential -which may point the way to new physics if, as appears currently, the potential is metastable. Beyond these motivations, for a generic model of new physics with boundary conditions fixed from the top down (such as supersymmetric models) it is important to know what regions of parameter space are allowed, compatible with the measured Higgs mass. For example, a one-loop calculation may naively lead to excluding certain constrained scenarios, whereas with a two-loop calculation the Higgs mass may be large enough; this is related to the difficulty in estimating the error in the Higgs mass calculation, since at two-loop order there are new contributions from particles that have no direct couplings to the Higgs, and a simple variation of the renormalisation scale as an estimate of the error is not sufficient.
In general, there are three approaches to tackle the problem of finding the Higgs mass precisely: (i) effective potential calculations, (ii) diagrammatic calculations, (iii) renormalisation group equation methods. We shall concentrate in the following on the first two options. Calculations from the effective potential suffer from a larger uncertainty compared to diagrammatic calculations because of the missing momentum contributions. However, these are only really pronounced at one-loop level, and it is already possible to calculate the full one-loop Higgs mass inlcuding momentum dependence for generic models using SARAH [4][5][6][7][8] to produce SPheno [9,10] output or SOFTSUSY [11] output via FlexibleSUSY [12]; explicit results have been known for some time for specific models such as the MSSM with real parameters [13][14][15][16] and complex [17][18][19]; and for the NMSSM with real parameters [20][21][22] and complex [23]. At two loops the momentum effects are expected to be small: according to recent calculations for specific models they are comparable to the experimental uncertainty [24][25][26] and since the momentum-dependent corrections due to new physics scale at best as m 2 H /M 2 New Physics relative to the effective potential contribution, we expect this to be a general result. Hence effective potential calculations, with their concomitant great simplification over the diagrammatic approach, should be useful at two loop order and beyond (even if the inclusion of the momentum dependence will ultimately be necessary to reach the experimental accuracy).
The state of the art in these calculations is however somewhat suprising given that the complete generic expressions for the two-loop effective potential, valid for a general renormalisable quantum field theory, have been available for more than ten years by the work of Martin [55]. These were applied to a complete two-loop calculation of the light Higgs mass in the MSSM in the effective potential approach in Ref. [39]. Furthermore, generic results for the diagrammatic calculation including the momentum dependence up to leading order in gauge couplings have been available in the literature for almost as long [56][57][58]. Unfortunately the results of Ref. [39] suffered the so-called "Goldstone boson catastrophe" (recently re-explored in [59,60]) due to the presence of tachyons in the tree-level spectrum so were numerically unstable. Perhaps due to this no public code was made available to exploit these prior to Ref. [61] where an implementation in SARAH/SPheno was presented. Currently, the only generic two-loop results relevant for the Higgs mass calculation still not present in the literature are the all-electroweak loops and the corrections to the Z-boson mass relevant for determining the electroweak expectation value. These will be the subject of future work. Here we shall instead continue the process started in Ref. [61] of making the pioneering generic results of Martin available in a public code -which entails performing some new calculations.
As we stated above, calculating the two-loop corrections to the Higgs mass in the effective potential approach is expected to be a good approximation. However, there is more than one way to actually perform even this calculation: either we can calculate the potential and numerically take the derivatives, as done in Refs. [38] and [61], or we can perform the calculation diagrammatically and set the external momentum to zero. In this work we shall exploit this equivalence: we shall analytically take the derivatives of the effective potential, producing expressions equivalent to the diagrammatic calculation and having the same structure, but with much simpler loop functions. The advantages of this over the first method are that the results are numerically stable 1 ; it is in principle a faster computation for more complicated models where the numerical method must make several passes to ensure stability; it can later be extended to a full diagrammatic calculation by simply changing the loop functions -but at zero momentum the loop functions are much simpler and therefore significantly faster to evaluate. We shall therefore compute the analytic expressions for the first and second derivatives of the two-loop effective potential and implement them in SARAH. As in Refs. [61,62] we shall ignore broken gauge groups, and adopt the same ansätze regarding the contribution of the electroweak gauge couplings to the tree-level Higgs mass matrix, to which references we refer the reader; the reasons for restricting to the so-called "gaugeless limit" are (a) partial circumvention of the Goldstone boson catastrophe (complete evasion in the case of the MSSM or any theory where the electroweak gauge couplings entirely determine the Higgs quartic potential); (b) significant simplification in the expressions and therefore speed in calculation; (c) the electroweak contributions are expected to be small, of the same magnitude as the momentum-dependent contributions. On the other hand, in the appendix we provide just the tadpole contributions in the case of broken gauge groups, and will return to the full expressions in future work.
The layout of the paper is as follows. In Sec. II we explain our procedure to take the derivatives of the effective potential to extract the two-loop tadpole functions for a general theory with massless gauge bosons in a form convenient for automation. In subsection II C we summarize our results for the tadpole diagrams; we present the full set of second derivatives in appendix B. The implementation of these results in SARAH, including some technical details of the translation of the generic results into an an algorithm, is explained in Sec. III before we conclude in Sec. IV. Impatient readers interested in using our implementation of the results might want to jump directly to subsection III B.

II. DERIVATIVES OF THE EFFECTIVE POTENTIAL WITH MASSLESS VECTORS
In this section we derive the expressions for the two-loop tadpoles in a general quantum field theory with massless gauge bosons in Landau gauge. To do this, we analytically take the dervatives of the expressions in [55]. Writing the couplings in the notation of that paper, the theory is defined by real scalars φ i , Weyl fermions ψ I and massless gauge bosons A a µ where the gauge covariant derivative for the fermions and scalars are The structure constants θ a ij are imaginary antisymmetric matrices that obey the gauge algebra but, since we are writing in terms of real bosonic fields, for complex representations they will have twice the dimension of the equivalent generators T a (so e.g. a U (1) generator is two-dimensional). We define as usual tr(θ a θ a ) = i d(i)C(i) where d(i) is the dimension of the representation of field i and C(i) the quadratic casimir, and similarly for T a .
The lagrangian is then composed of the normal kinetic terms of the scalars and fermions using the above covariant derivatives supplemented by purely scalar and scalar-fermion interactions y is in general a dimensionless complex tensor with y IJk = y JIk , while λ ijk , λ ijkl are real, symmetric tensors.
A. Effective potential We can simplify the expression for the effective potential in the Landau gauge given in Ref. [55] for our case; with all gauge groups unbroken some diagrams do not contribute. The non-vanishing diagrams are shown in Fig. II A. The contribution of each diagram to the effective potential is given by: Here, y and λ are the trilinear and quartic couplings of eq. (II.2), g is a gauge coupling, and M are fermion masses. The loop functions, given in terms of standard basis functions given in appendix A, are the same for M S and DR in the case of f SSV =(x + y) 2 + 3(x + y)I(x, y, 0) + 3J(x, y) − 2xJ(x) − 2yJ(y) (II. 10) but differ for those with vectors and fermions: 6I(x, y, 0) (II.11) All of these functions are symmetric on the substitution of their first two indices, but may not be so with the third. In fact, we can then combine the vector-fermion diagrams to give where δ MS is zero for DR and one for MS.

B. Derivatives of the potential
Here we shall analytically take the derivatives of the potential. In a generic model, we may want the derivatives of the Higgs potential in terms of some unrotated fields, i.e. in a basis where the mass matrix is not diagonal; let us say that we start with such a case. We write the tree-level potential in terms of some expectation valuesv i and the associated real fluctuations S 0 i as Then there is a tree-level rotation to diagonalise the mass matrix; we then obtain We can write φ i = v i + S i and work with the couplings of eq. (II.2). We then need the quantities which enter in the effective potential calculation, which are masses and couplings depending on {S i }.
In general, we have with the shorthand notation S ≡ {S i }. Hence we can write Similarly for the fermions we have However, for the purposes of the effective potential, we then require a further diagonalisation for m 2 ij (S): we put We name the couplings in this basis asm i (S),λ ijk S (S),λ ijkl (S). We then express the effective potential by inserting the couplings and masses in the basis {S i } into the formulae of eqs. (II.3)-(II.9). However, to take the derivatives we rewrite the expressions in terms of the basis {S i }, and use the trick (with (II.21) For similar expressions we write by abuse of notation (using C ≡ 16π 2 µ 2 (2π) −d [57]) (II.22) For fermion propagators we can write (II.23) Let us demonstrate our method on a brief example: Note that the R 0 ra are not functions of S i and so do not present complications if we want to take futher derivatives.
We can similarly take the derivatives of all the remaining loop functions; these are derived from the basis: The derivative of a typical term in the effective potential will have the form where, generalising the above, it is straightforward to show that for a generic function appearing in the effective potential f α composed of polynomials (even containing monomials with negative exponents) multiplying the loop functions above, we can write and similarly for permutation of the indices. On the other hand, this explicit expression is often inconvenient in practice due to the need to carefully take the smooth limit when x = u; it is instead more practical to rewrite the right hand side in terms of our basis of loop functions multiplied by suitable polynomials. In the following we present explicit expressions for the first derivatives (and, in the appendix, the second derivatives) which have been appropriately simplified to remove the apparent singularities.

C. First derivatives
Here we gather the set of two-loop tadpole diagrams. There are only three topologies, two of which only apply to the all-scalar case. All generic possible diagrams are given in Fig. II C. Here we present our results for the first Tadpole diagrams at the two-loop level which don't vanish in the gaugeless limit. derivatives in the basis {S i }, so without the rotation matrices R; the full tadpole in the original basis is then given by with the tadpoles on the right-hand side to be defined below.

All scalars
We start with the purely scalar diagrams which are in the first row of Fig. II C. The entire contribution is given by The new loop functions are defined as to W SSSS of Ref. [56] in the limit of zero external momentum.

Scalars and fermions
We have, first, the diagrams with two scalar propagators: Then we turn to one scalar and three fermion propagators: Here we have defined

Diagrams with vectors
Finally, the two generic diagrams involving vectors are given by

D. Second derivatives
To find the second derivatives of the potential we can apply the same technique. However, in principle, we can simply use the results of [56], which computed (diagrammatically) the two-loop scalar self-energies at leading order in gauge couplings. Since we want the self-energies for neutral scalars, this comprises all of the contributions, and if we want the results in the effective potential approach we can simply set the external momenta to zero (and, for this work, the masses of the gauge bosons to zero). In fact, for the majority of the diagrams, these yield the same result. However, in a few cases we find that by taking the derivatives of the potential we find simpler results (which are of course entirely equivalent). The full result is given by where the superscripts correspond to the numbers of scalars, fermions and vectors with types of topology listed for diagrams (M ) and (V

III. IMPLEMENTATION IN SARAH
We have implemented the new routines in SARAH. By including the first and second derivatives of the effective potential using the analytic expressions here, rather than numerically taking the derivatives of the couplings and masses as performed in the previous version [61], we can guarantee greater numerical stability, accuracy and speed improvements. Moreover, this approach allows a straightforward upgrade to the pole mass calculation by simply changing the loop functions called to those defined in Ref. [56] based on loop functions implemented in TSIL [63], which will be made possible in a future version.

A. Method
For any given supersymmetric 2 model $MODEL, once the user has specified the particle content and their symmetries, SARAH calculates all of the vertices and masses. It then writes a Fortran code (placed in the suggestively named 2LPole $MODEL.f90) which implements our expressions, linking to a static Fortran code (named 2LPoleFunctions.f90) of the basis functions for the generic first and second derivatives of the effective potential defined in section II and the appendix. These two pieces of Fortran code are called by SPheno during the calculation of the loop corrections to the Higgs mass. Here we shall give a few details of how SARAH writes 2LPole $MODEL.f90. The overall algorithm is to 1. Generate masses and couplings for all relevant particles in the gaugeless limit.
2. Populate and classify all tadpole topologies according to particle content.
3. For each tadpole topology, pass the set of diagrams along with information specifying the symmetries to a generic writer function.
4. Rotate the total tadpole vector by the Higgs rotation matrix to the non-diagonal basis, c.f. eq. (II.30).

5.
Populate and classify all second derivative topologies according to particle content.
6. For each mass topology, pass the set of diagrams along with information specifying the symmetries to a generic writer function.
7. Rotate the mass matrix to the non-diagonal basis (as with the tadpoles).
The writer function is actually identical for tadpoles and mass diagrams with a switch to adjust the number of Higgses. It cycles through the list of diagrams and applies the following process for each: (a) Determine symmetry factor of diagram due to permutations.
(b) Determine the colour factor; for diagrams with a gluon propagator this is simply d(I)C(I) whereas otherwise we must trace over colour indices of the vertices. In principle, for four-point vertices there can be two colour structures for the vertex which superficially leads to more than one colour factor for such diagrams. However, as we can simply see by inspecting the expressions in the appendix, or by considering that the colour factors have to be inherited from a corresponding vacuum diagram (since differentiating with respect to neutral Higgs fields cannot introduce any additional colour factor), for diagrams with a four-point vertex consisting of four coloured fields the colour factor is given by a trace over the indices in pairs. Hence such four-point vertices are saved with the colour factor of the pairs of indices traced over. To be more explicit, such vertices can only contribute if they come from differentiating V SS which contains the coupling λ iijj . We then must simply take care that the indices of the vertices correspond correctly to the indices that are traced over.
(c) Write a nested set of loops to sum the diagram over the generations of all particles and, for the inner loop, the external Higgs legs, since the most computationally demanding aspect is evaluating the loop functions and this can be evaluated before calling the Higgs loop -indeed, we also check that the coupling multiplying it is non-zero first too.
There are subtleties in translating our results into a form usable by SARAH, both stemming from the fact that in the SPheno code the couplings are stored in terms of either real or complex scalars, and four-component spinor fermions, while, since our results are based on those of Refs. [55] and [56] and for economy we use real scalars and two-component fermions. The translation between the two bases as required by SARAH and SPheno is largely as described in Ref. [61], however here we have the additional complication for fermions of translating chains of couplings and masses such as In SARAH, the interactions of Weyl spinors ψ are derived from the corresponding Dirac spinor interactions Ψ as and so c L (I, J, m) = c * R (J, I, m). For a given topology, SARAH populates the diagrams using Dirac propagators (i.e. links fermions with its conjugate) and so for the coupling C 1 above we will find sets of particles Suppose each of the fermions is a Dirac spinor with Weyl spinors ψ I L,R etc, then to construct coupling C above we must sum over the left and right-handed Weyl fermions (which have opposite representations of all gauge groups) and thus (noting that SPheno always internally stores the fermion masses as real positive definite) where M = 1 for Majorana fermions and zero otherwise. A final point regarding the translation into a basis of real and complex scalars is that the new routines assume that there is a unique way of constructing a gauge-and global symmetry-invariant coupling λ ijk from complex scalars other than the complex conjugate of the whole coupling; i.e. if λ ijk is permitted for given complex i, j, k then λ ij k is not. This is evidently true -if both are permitted then we can generate a holomorphic mass term at one loop which violates the premise. However, it is important in that we cannot write gauge singlets as complex scalars if they have couplings violating the above condition, no matter how small the couplings -for example for sneutrinos in see-saw models.

B. How to use the new routines
To study a model with SARAH the general procedure is as follows: the user should download and run SARAH with the demanded model. SARAH derives all analytical expressions for mass matrices, vertices, renormalisation group equations as well as loop corrections and exports this information into Fortran source code. The Fortran source code is compiled together with SPheno and all numerical calculations are then performed by the new SPheno module. This includes a calculation of the entire mass spectrum, branching ratios as well as flavour and other precision observables [64]. For the mass spectrum all one-loop corrections to any particle are included in a diagrammatic way [21]. For a supersymmetric model there are now in addition three options to get two-loop corrections in the Higgs sector. The first two are based on the effective potential approach presented in Ref. [61] while the new routines are called by the (now default) third option.
A step-by-step description to obtain a spectrum generator for an arbitrary model $MODEL implemented in SARAH reads as follows: . Also some hard-coded corrections are available which are based on results in literature: 8->8 uses the known α S (α b + α t ) corrections for the MSSM, NMSSM, TMSSM or any variant thereof with up to four neutral CPeven Higgs fields and including models with Dirac gauginos [65]. 8->9 uses the corrections of option 8 and adds the two-loop MSSM (α t + α b + α τ ) 2 results based on Refs. [35-37, 40, 41]. Note that the last two options are not included by default in the SPheno output of SARAH. To include them, the user must make sure to include in the SPheno.m of the considered model We have intensively used the SPheno output to validate our new two-loop functions, in particular: • We found a numerical agreement of more than 10 digits between our code and using public routines for the MSSM based on Refs. [35-37, 40, 41] for the self-energies. In order to perform this validation, it is necessary to use the same assumptions: turn off the first and second generation Yukawa couplings; take the Goldstone boson and light Higgs masses in the loops to be zero, and set the tree-level mixing angle of the neutral CP-even scalars to α = β − π/2. The excellent agreement between all four possibilites to calculate the two-loop Higgs masses in the CMSSM is also shown in Fig. 5.The parameter point used here is the same as in [61], • Similarly, we found full agreement with available results for the α s (α b + α t ) corrections in the NMSSM [20] and for Dirac gauginos [65].
• We compared the full two-loop corrections, i.e. also including corrections not involving the strong interaction, for the NMSSM, models with Dirac gauginos as well as for the B-L-SSM [66] against the results using the other two options based on a completely independent implementation in SARAH presented in Ref. [61]. We usually found very good agreement. Tiny differences were based on numerical artefacts in the routines using the effective potential ansatz. Similarly, we could reproduce the results of Ref. [67] for the two-loop contributions to the Higgs mass stemming from R-parity violating couplings.
The new routines of course provide better stability. For example, in the routines based on numerically taking the derivatives of the potential, it is necessary to take care with the initial step size; if there are neutral scalars which have small expectation values then the results from those methods could become inaccurate -this problem occurs in general for any neutral scalar having expectation value v i M SU SY . Less significantly, the numerical method can suffer (small) errors when there are small couplings present, such that they do not induce a sufficient shift in particle masses or couplings upon variation of the Higgs vevs to accurately take the derivative. Hence, it is very important to have two independent implementations of generic two-loop Higgs mass calculations in SARAH/SPheno: this is the only possibility to cross check results for models beyond the (N)MSSM at the moment. Thus, we highly motivate users to test all options for the model under consideration and to compare the results.

IV. SUMMARY
We have presented the derivation of a new set of expressions for calculation of the tadpoles and self-energies at the two-loop level. These expressions include all generic diagrams which do not vanish in the gaugeless limit and are valid in the limit of zero external momenta. This set of loop functions is simpler than the set of expressions obtained by taking the limit p 2 → 0 in the pole mass functions available in literature so far. This allows for a rapid numerical evaluation of the Higgs mass. We have implemented these functions in Fortran and included them in the new version of SARAH 4.5.0. This provides the possibility to automatically calculate the Higgs mass in a wide range of supersymmetric models with a guaranteed numerical accuracy and stability. The obtained precision for the Higgs mass is comparable with the one dedicated spectrum generators provide so far for the MSSM, and can now be applied to the study of a wide variety of models.
Aside from accuracy and stability, one of the principal advantages of this approach is that it is readily extendable. It would be straightforward to extend the calculation to non-zero momentum by changing the functions in the code and linking with the library TSIL. On the other hand, including the electroweak contributions should be possible by applying these techniques to the full effective potential; we presented the expressions in this case for the tadpoles in the appendix, but the second derivatives are currently unknown -as are the full set of equivalent expressions in the diagrammatic approach. Furthermore, to truly reach the full two-loop precision we would require the two-loop shift in the Z-mass that determines the electroweak expectation value. We hope to return to these issues in future.
Ref. [57] and make often use of results presented there. Note that we use the standard notation where Q is the renormalisation scale.

a. One loop functions
At the one-loop level only two functions are needed: Clearly, we find the following relation to J(x) which is widely used at two loops: where we have introduced the subscript to denote that the external momentum is zero. Let us also denote for future use which is actually symmetric on all three indices. T 0 (x, y, z) = − ∂ ∂x I(x, y, z) (A.8) These functions have to following properties: (i) I(x, y, z) is symmetric in all arguments; (ii) T 0 (x, y, z) is symmetric in the last two arguments; (iii) U 0 (x, y, z, u) is invariant under the exchange z ↔ u and x ← y; (iv) M 0 (x, y, z, u, v) is invariant under the interchanges (x, z) ↔ (y, u), (x, y) ↔ (z, u), and (x, y) ↔ (u, z). An explicit expression for I(x, y, z) is for instance given by [55] I(x, y, z) Note that x, y ≤ z has been assumed here.
c. Relations required for the pole functions The above loop functions are used as a basis for the various loop functions. However, we also find additional combinations such as .
(A. 16) we should compare this to Hence we can write V SSSSS (x, y, y, u, v) = −V (x, y, z, u) and remember that V SSSSS is symmetric on its first three entries.

d. Simplified loop functions
For the amplitudes with one or more massless (and possibly identical) fields we find simplified expressions for the loop integrals, some of which are collected below.
It is also sometimes necessary to consider the case of small mass splittings. The results for I(δ, x, y), I(δ, x, x), I(δ, δ , x) in the limits δ → 0, δ → 0 can be found in Ref. [55] and we do not repeat them here.
Appendix B: Second derivatives of the effective potential In this appendix we present the results for the second derivatives of the effective potential. Largely these are identical to those in [56] with the external momentum set to zero, but for sake of completeness we repeat the full set here. However, certain expressions become much simpler in this limit, notably some complicated functions involving both fermion and scalar propagators, and those involving gauge bosons.
The full contribution is

Diagrams with only scalar propagators
The first contribution, including only scalar propagators, comprises the eight diagrams shown in 3. These are unchanged from [56] and are given by Here the loop integral functions are given by: It should be noted (for example, for the purposes of evaluating the colour factors) that topologies X SSS , Y SSS , Z SSS arise from differentiating V (2) SS (and hence the tadpole T SS ) while the others arise from differentiating V (2) SSS (and hence the tadpoles T SSS and T SSSS ).

Diagrams with scalar and fermion propagators
The contributions from diagrams with the topology W are where we can slightly simplify the loop functions: W SSF F (x, y, z, u) = −2W SSSS (x, y, z, u), (B.13) W SSF F (x, y, z, u) = −(z + u − y)U 0 (x, y, z, u) − I(x, z, u) + B 0 (x, y)(J(z) + J(u)). (B.14) The contributions from diagrams of the topology M with four fermions are The results from diagrams of the topology M with three fermions are The results from diagrams of topology V with four fermions are These represent significant simplifications over the full pole contributions. To recapitulate, y, z, u, v). (B.35) Note that f (2,0,0) F F S is symmetric on its first three indices.

Diagrams with one vector propagator
For self energies of neutral scalars where all gauge groups are unbroken, the diagrams involving one vector propagator are particularly simple.

Effective potential
The full two-loop effective potential in the Landau gauge was given in Ref. [55]: