Automatised matching between two scalar sectors at the one-loop level

Nowadays, one needs to consider seriously the possibility that a large separation between the scale of new physics and the electroweak scale exists. Nevertheless, there are still observables in this scenario, in particular the Higgs mass, which are sensitive to the properties of the UV theory. In order to obtain reliable predictions for a model which involves very heavy degrees of freedom, the precise matching to an effective theory is necessary. While this has been so far only studied for a few selected examples, we present an extension of the Mathematica package SARAH to perform automatically the matching between two scalar sectors at the full one-loop level for general models. We show that we can reproduce all important results for commonly studied models like split- or high-scale supersymmetry. One can now easily go beyond that and study new ideas involving very heavy states, where the effective model can either be just the standard model or an extension of it. Also scenarios with several matching scales can be easily considered. We provide model files for the MSSM with seven different mass hierarchies as well as two high-scale versions of the NMSSM. Moreover, it is explained how new models are implemented.


Introduction
The Standard Model (SM) of particle physics is a very successful theory which has been completed with the discovery of the Higgs boson at the Large Hadron Collider (LHC) [1,2]. On the other side, there are observations like dark matter for which no viable candidate exists within the SM. While it has been expected that solutions to the open problems of the SM, like e.g. supersymmetry (SUSY), exist close to the eleca e-mail: martin.gabelmann@kit.edu b e-mail: milada.muehlleitner@kit.edu c e-mail: florian.staub@kit.edu troweak scale, the LHC has not found any direct signal for new physics so far. Therefore, the possibility of a large gap between the electroweak (EW) and the scale of new physics has been studied more intensively in the recent years. The most prominent idea in this direction is 'split supersymmetry' (split-SUSY) in which the SUSY scalars are much heavier than the SM particles and the SUSY fermions [3][4][5]. In this setup, most of the appealing properties of SUSY like gauge coupling unification and a dark matter candidate are kept, but the coloured particles are too heavy to be produced at the LHC. Mechanisms have been proposed how split-SUSY could arise from string theory [6,7], and also the question of naturalness has been discussed [8]. Moreover, the ansatz of high-scale SUSY, i.e. that all SUSY particles are much heavier than the EW scale, is taken seriously nowadays [9,10]. While it is widely believed that these models suffer from a large fine-tuning, it has pointed out that large SUSY scales can be combined with the relaxion mechanism to solve the big and the small hierarchy problem simultaneously [11]. The idea of SUSY with very large mass scales is not restricted to the Minimal Supersymmetric extension of the SM (MSSM), but has also been applied to other SUSY models like the Nextto-MSSM (NMSSM) [12] or models with Dirac gauginos [13][14][15][16][17].
Even if states beyond the SM (BSM) are too heavy to be produced at current colliders, they often still have an in-print in experimental results, see e.g. Refs. [18,19]. The precise measurement of the Higgs boson mass of m h = 125.09 GeV [20] at the LHC has added another very important constraint in this direction. Consequently, large efforts were put in a precise Higgs boson mass calculation in split-or high-scale SUSY [9,10,[21][22][23]. The reason for this endeavour is that the commonly used fixed order calculations of the Higgs boson mass in SUSY models should only be applied in the case of a small separation between the EW scale and the SUSY scale. Otherwise, the presence of large logarithms introduces a large uncertainty in the prediction of the numerical value of m h [23][24][25][26]. This can be resolved either by the standard ansatz of an effective field theory (EFT) in which the heavy states are integrated out [27][28][29][30][31][32][33][34][35][36], or by a hybrid method in which the fixed-order calculation is combined with the higher-order leading logarithms extracted from an EFT [37][38][39][40]. In both cases, one needs to know how the couplings among the light states depend on the full theory. In terms of the EFT ansatz this means that the full model involving heavy and light states must be matched to an effective theory at the scale at which the heavy degrees of freedom are integrated out. The matching at leading order is straight-forward and the relations often can be read off from the tree-level Lagrangians of both models. However, tree-level relations are usually not sufficient to obtain the necessary precision in the Higgs boson mass prediction. Therefore, higher-order corrections are needed. Of course, the matching procedure at the full one-loop level is already much more time-consuming. Depending on the details of the full and effective model also several subtleties like infra-red divergences can occur as discussed in Ref. [41].
In order to facilitate these studies, we have developed an automatised process to perform the matching between the scalar sectors of two renormalisable theories. This feature has been implemented in the Mathematica package SARAH [42][43][44][45][46] and provides the functionality to obtain analytical expressions for the matching conditions at the oneloop level. Also the interface between SARAH and SPheno [47,48] has been extended to include the matching between an EFT and a UV-complete theory. In that way, one can obtain very quickly numerical predictions for the Higgs boson mass but also for all kind of other observables that concern the Higgs boson. It is worth to stress that this functionality is not restricted to split-or high-scale versions of the MSSM. A large variety of SUSY, but also non-SUSY models, with large BSM scales can be studied with the presented tool-chain. Also the considered EFT need not to be the SM, but could be a singlet extension, a Two-Higgs-Doublet-Model (THDM), or an even more complicated model. Concerning the nature of the heavy states, we restrict our attention to heavy fermions and scalars. The implementation of integrating out heavy vector bosons at the one-loop level is reserved for future work. However, the low-energy EFT can still contain an extended gauge sector which is also matched at the one-loop level. Nevertheless, we will mainly concentrate in the given examples on the established MSSM scenarios because they offer the possibility to compare our generic approach with results available in the literature.
This paper is organised as follows. In Sect. 2 we explain our generic matching procedure. In Sect. 3 the new routines in SARAH are explained while a comparison with the literature is done in Sect. 4. We summarise in Sect. 5.

General ansatz
We consider a general, renormalisable gauge theory with a set of scalars {φ i } and fermions {ψ i } charged under unspecified (sub-)sets of the theories gauge group. Without loss of generality, one can always assume that the scalars are real. The Lagrangian can be written as where all gauge and representation indices have been suppressed. The covariant derivative D μ and the gauge fields are chosen such that the field strength tensors {F μν a } form diagonal kinetic terms (in case of multiple gauged U(1) groups). In the following it is always assumed that all gauge groups are broken near the scale of EW symmetry breaking. If particles with very different masses appear in such a theory, one can categorise the particle content into light fields ( Integrating out all heavy fields leads to an effective theory which contains only light degrees of freedom where the last line contains operators with dimension greater than four. Concerning a precise prediction of Higgs boson masses, only purely scalar operators with ascending dimensionality may be of interest for the matching. However, for d > 4, their influence on the scalar potential is of the order v i/M j , where v i is the vacuum expectation value (VEV) of a light and M j the mass of a heavy field, v i M j . Supposed that the fundamental theory is renormalisable, it follows from the decoupling theorem, that the higher-dimensional operators become unimportant if M j → ∞. The question arises, at which scale the v i/M j terms are no longer relevant for a precise Higgs boson mass calculation. The impact of dimension-six terms, compared to ordinary threshold corrections (d ≤ 4), on the Higgs boson mass in a matching of the SM to the MSSM was studied in Ref. [49]. It was found that for 500 GeV < M j < 1000 GeV, a two-loop matching of these operators yields corrections on m h in the sub-GeV range, which rapidly drop for M j > 1 TeV. Since the focus of this work is on BSM scenarios with M j 1 TeV, we neglect all v i/M j contributions during the matching. Thus, we assume that all VEVs responsible for the breaking of a low-energy gauge theory can be neglected compared to the masses of the heavy states. In particular, this means that all gauge bosons as well as chiral fermions are treated as massless in the computation of the matching conditions. All information about the heavy states is encoded in the effective couplings and massesλ,κ,m 2 ,M andỸ . The purely scalar interactionsλ,κ and mass squaredm 2 contain the crucial information about the scalar sector of the EFT, hence, they have the biggest impact on the properties of the light scalars. We know today, that (at least) one of these light scalars must have couplings comparable to the predictions of an SM-like Higgs boson and the mass must be about 125 GeV. Thus, even if the mass scale of the heavy fields is well above the reach of the LHC, we can test if the fundamental UV theory is consistent with the Higgs boson mass measurements through a precise calculation of the effective couplings at the matching scale and the Higgs properties at the weak scale. In order to determine the effective couplings in terms of parameters of the UV-theory, one assumes the matching condition that the n-loop m-point amplitudes involving the same external (light) states must yield the same result in the infra-red (IR) regime of the UV-theory (i.e. the scale where the heavy fields are integrated) and the EFT, Note, that the external fields in the two theories to be matched must be treated equally. Thus, additional wave-function renormalisations involving internal heavy fields may con-tribute to Eq. (4) by also matching the first derivative of the 2-point function w.r.t the external momentum of the light fields.
In this paper, we are going to calculate M using the Feynman diagrammatic approach neglecting all external momenta. The tree-level matching condition for a quartic coupling, is depicted in Fig. 1. Due to the assumption of vanishing external momenta and vanishing light masses, infra-red divergences appear on both sides of Eq. (5). Since the treelevel matching for cubic couplings is trivial, the divergences cancel exactly. Thus, the effective quartic couplingsλ abcd are given bỹ As already mentioned, the matching at tree level is not sufficient for a precise prediction of the properties of the scalar sector at the low-energy scale. Thus, one needs to include loop corrections changing the matching conditions tõ where δ X (φ L , φ H , ψ L , ψ H ), with X = λ, κ, denote the sum of all one-loop contributions that contain only light fields, mixed heavy and light fields as well as only heavy fields in the loop . Likewise δλ can only arise from diagrams involving light fields in the loop since there are no heavy states present in the EFT. All generic diagrams which can contribute to treelevel and one-loop amplitudes of any renormalisable scalar operator are given in Appendix A. Again, IR divergences caused by light fields are present on both sides which need to cancel in the matching conditions. A detailed discussion on these cancellations is beyond the scope of this paper but was recently discussed in Ref. [41]. In summary, the matching condition can be expressed in terms of IR-finite pieces where the one-loop contributions δ X are computed using modified loop integrals where the IR divergent pieces have been subtracted. For instance, the scalar two-point integral B 0 with vanishing external momentum (for simplicity we omit the vanishing external momentum in the argument of all loop function) and vanishing masses suffers from a logarithmic IR divergence which will necessarily cancel in the matching condition Eq. (4). Thus, the replacement of the B 0 with the modified loop function makes this cancellation manifest without the need to compute the corresponding IR-divergent diagrams in the EFT. Thus, the calculation of the matching conditions can be performed in a straight-forward way by using the IR-safe loop functions B 0 , B 1 , C 0 , D 0 ,Ḃ 0 andḂ 1 defined in "Appendix B".

Renormalisation scheme
A simple renormalisation scheme which is applicable to a wide range of models is the MS/DR scheme. Therefore, we are going to stick mainly to this scheme. The only exception is the treatment of the off-diagonal wave-function renormalisation (WFR) of the scalar fields. It has been proposed in Ref.
[10] that these contributions can be dropped by assuming finite counter-terms for some input parameters. For instance, in the high-scale MSSM one could assume a counter-term for tan β which exactly cancels the off-diagonal WFR contributions. This approach is a more economic calculation and can lead to performance improvements in the runtime. However, it depends on the considered model and the chosen input parameters if such a scheme is possible. Therefore, we provide the possibility to include or exclude the one-loop contributions from off-diagonal WFR constants during the calculation.
For an appropriate choice of the WFR treatment it is worth to mention the equivalence between excluding the offdiagonal WFR constants and the extraction of effective quartic couplings from a pole-mass matching [24,41]. Thus, for the comparison with tools that use a pole-mass matching, the inclusion of off-diagonal WFR constants should be disabled in the calculation.

Parametrisation of the results at the matching scale
Using matching conditions to calculate the effective couplings yields solutions that are functions of the parameters of the UV theory. However, in some cases it might be better to (at least partially) give their dependence on the EFT parameters. This is especially the case for the SM gauge and Yukawa couplings because their values are known very precisely. Therefore, one also needs to match these couplings at a suitable loop-level. Concerning the matching of the scalar sector, the EFT parameters that enter the scalar matching conditions at tree-level need to be matched at the one-loop level (and re-inserted into the scalar tree-level matching). For all other parameters, a tree-level matching is sufficient as long as we stick to a one-loop matching of the scalar couplings. For non-supersymmetric models the scalar parameters which we want to match are free parameters, i.e. in these cases a matching of the SM parameters at tree-level is always sufficient. This is different for supersymmetric models because the scalar couplings are related to the other couplings through F-and D-terms. We concentrate on the D-term contributions, i.e. the matching of the gauge couplings, because this is the part important for the matching of scalars that couldat least in principle -provide a SM-like Higgs boson. The matching of the gauge couplings is parametrised by and receives two different contributions: 1. Thresholds from heavy fields: where g i is the gauge coupling with respect to the gauge group i, I i 2 (x) is the Dynkin index of the field x with respect to the gauge group i, and C i (x) is a multiplicity factor taking into account the charges under non-Abelian gauge groups others than i, i.e. in the case of the SM gauge group, this counts the colour and isospin multiplicity in the loop. 2. MS-DR conversion: required if an MS and a DR renormalised quantity are to be matched. This is e.g. the case if non-SUSY models are matched onto SUSY ones. There are two different contributions which affect the quartic couplings: -The finite shifts of the gauge couplings for an SU (N ) group are [50] -Quartic vertices receive an additional shift from MS-DR conversion from the diagrams shown in Fig. 2.
The amplitude difference of this diagram between the two schemes is where c 1 and c 2 are the two involved vertices between two scalars and two vector bosons.
The calculation of the two different contributions was implemented in SARAH and are automatically included in the matching procedure.
where Y abc is for instance a running SM Yukawa coupling and δY abc contains corrections from diagrams containing heavy fields, which are obtained with IR-safe loop functions, as well as MS-DR conversion if necessary. If the EFT is not the SM but an extension with additional fermions, also new Yukawa-like couplings are present below the matching scale. A good example for such a scenario is for instance split-SUSY with effective gaugino-Higgsino-Higgs couplings. Of course, the one-loop relation to calculate these couplings is just given by inverting Eq. (18), i.e.
Thus, in a generic approach, both types of Yukawa coupling corrections, above (SM-like) and below (BSM-like) the matching scale are obtained simultaneously. Necessary ingredients are the one-loop diagrams depicted in Fig. 3 together with the wave-function corrections of the external states.

Implementation in SARAH and SPheno
In the last section all necessary ingredients for a matching of two arbitrary renormalisable scalar sectors at the one-loop level were introduced. In this section we describe the implementation as well as the usage in the computer programs SARAH and SPheno.

General information about SARAH and SPheno
SARAH 1 is a Mathematica package optimised for an easy, fast and exhaustive study of BSM models. For a given model, which is defined in form of three input files, SARAH derives all tree-level properties, i.e. mass matrices, tadpole equations and vertices. Moreover, the analytical calculations of one-loop self-energies and tadpoles as well as of two-loop renormalisation group equations (RGEs) are fully automatised in SARAH based on generic results given in literature [51][52][53][54][55][56][57][58][59][60][61]. With version 3, SARAH became the first 'spectrumgenerator-generator': all analytical information derived by SARAH can be exported to Fortran code which provides a fully-fledged spectrum generator based on SPheno. A SARAH generated SPheno version calculates all masses at the full one-loop level, and includes the dominant twoloop corrections for neutral scalars [62][63][64]. Beyond that, SPheno makes predictions for two-and three-body decays, flavour and precision observables [65,66], and the EW finetuning. In order to define the properties of the generated SPheno version, SARAH needs an additional input file usually called SPheno.m. This input contains the following information: -The input parameters of the model -The choice for the renormalisation scale -The boundary conditions at the electroweak scale, at the renormalisation scale and at the GUT scale -Optional: a condition to dynamically determine the GUT scale, e.g. g 1 (m GUT ) = g 2 (m GUT ) -A list of particles for which the two-and three body decays should be calculated.
Since the SPheno.m file will be important for the discussion in the following, we give an example in Appendix C how such a file may look like. For more details, we refer to the manual 1 SARAH is available at hepforge: sarah.hepforge.org.
as well as the SARAH wiki page. 2 In the following section, we discuss various aspects that arise in an automatised matching between two models and how they have been considered through the implementation of two independent approaches.

Available options to perform the matching
The matching of two scalar sectors can be motivated by a precise investigation of very different properties of the theories to be matched. The largest contributions to threshold corrections often have their origin in one common sector of the heavy spectrum. It can be of particular interest to track this origin down in order to learn more about which parts of a given UV-theory are essential for the predictions in an EFT framework. For this purpose an analytical evaluation of threshold corrections is preferred. The analytical solutions can also easily be ported to other computer programs which is a key feature of many existing SARAH routines. As already discussed, the matching of an EFT onto a UVcomplete model does not only influence many low-energy observables but also enters the RGE running and other predictions above the matching scale. Considering the whole picture of the matching procedure and its numerical influence in all sectors of the theories to be matched, a numerical calculation of threshold corrections is preferred because it can easily be embedded into existing routines of the generated SPheno code.
With SARAH version 4.14.0 we provide two different possibilities to perform the matching between two arbitrary scalar sectors: 1. An analytical calculation within Mathematica 2. A fully numerical calculation using only the SPheno interface.
It is important to stress that both options are not based on the same routines, but have been implemented independently. Thus, they offer the possibility to double check the obtained results. A schematic comparison of the two approaches is  Fig. 4. A summary of the description given here is also available at the SARAH wiki page. 3 For the analytical calculation it is necessary that all mass matrices in the model can be diagonalised analytically. Thus, it is usually necessary to work with a set of assumptions which simplify the most general mass matrices in a given model. In theories with spontaneous symmetry breaking, a high degree of different mixing patterns is introduced through the presence of VEVs. It has already been argued that, if these VEVs are responsible for the generation of masses in the EFT, i.e. if the low-energy Lagrangian is invariant under the symmetries broken by these VEVs, a common assumption is to neglect all small VEVs. In addition, flavour violating effects are usually negligible. The only exception are scenarios in which large contributions to flavour violation occur in the new physics sector. This could for instance happen in the MSSM with large off-diagonal trilinear soft-terms which can have a big effect on the Higgs boson mass [67]. Thus, if any of these assumptions is not justified, it is necessary to switch to the purely numerical calculation.
Although the focus of the Mathematica interface is on the derivation of analytical expressions for the matching conditions, additional routines have been implemented to make these results easily usable in numerical calculations. This has the advantage that the obtained code for numerical evaluations can be much faster than the fully numerical interface because many simplifications can be performed on the ana-lytical level. In addition, the obtained results can be exported into L A T E X files which makes a evaluation of the expressions in a human readable format possible. On the other hand, the fully numerical implementation has several advantages: (i) the RGE running above the matching scale can be performed, (ii) shifts to fermionic couplings can be included, (iii) several EFTs appearing in models with more than one matching scale can be automatically linked (iv) flavour violating effects can be included.
Before describing the user interface of the new routines, we want to comment on a few subtleties to provide a better understanding on the importance of certain user inputs.
Therefore, the correct matching condition to calculate λ SM becomes where κ i denotes tree-level vertices in the UV theory whileδ are the corresponding one-loop shifts. The relative normalizations between operators in the considered UV and the effective theory, such as for example the factor − 1 /3 in Eq. (22), have to be provided by the user. 3. Superposition of fields: when matching a scalar sector involving multiple (light) scalar fields with identical quantum numbers, often linear combinations of external fields contribute to the matching of different parameters. For instance, consider the couplings λ 4 and λ 5 in a THDM: where the two SU (2) doublets H 1 and H 2 have the same hypercharge. We find that any vertex involving λ 4 receives also contributions from λ 5 and vice versa. For instance consider the couplings after splitting the two doublets into their charged (H ± 1,2 ), CP-even (h 1,2 ) and CP-odd (A 1,2 ) components (note that the gauge eigenstates introduced here also correspond to the mass eigenstates as we assume vanishing VEVs). For simplicity, we assume real parameters . Thus, to obtain the matching conditions for λ 4 and λ 5 separately, it is necessary to calculate the superpositions These conditions are user input as well.

Analytical approach
In order to use the Mathematica interface to obtain analytical expressions for the matching conditions, one needs to initialize a Mathematica kernel, load SARAH and start the considered high-scale model. This can be done by opening a new Mathematica notebook and entering the commands In [1] << SARAH . m where <Model> could be for instance MSSM or NMSSM. In the next step, there are two possibilities to obtain the matching conditions analytically: 1. one can calculate individual effective couplings in an interactive mode or 2. use a batch mode to calculate several matching conditions at once and to optionally obtain L A T E X , Fortran and SPheno outputs.
We are going to give details about both options which are based on the new command where possible options are If the interactive mode is demanded, the option InputFile has to be omitted while values for Parametrisation, Assumptions and SolveTadpoles should be provided to allow for an analytical diagonalization of all mass matrices. The usage of the batch mode requires only the option InputFile and serves a high reproducibility of the obtained results by providing only one single input file.
The provided assumptions and parametrisations are used to calculate analytical expressions for all masses and rotation matrices. If this is not possible, because Mathematica cannot diagonalize the mass matrices analytically (using the build-in functions Eigensystem and SingularValueDecomposition), one can either use the purely numerical interface explained in Sect. 3.4 or choose appropriate simplifying assumptions.

Interactive mode: calculating individual matching conditions
Initializing the matching routines using the InitMatching function with the options described in the previous paragraph while not specifying the option InputFile enables the interactive mode. The necessary vertices of the high-scale theory are calculated or loaded from a previous session and the masses/rotation matrices are derived. However, no further calculations are performed at this point. Example initialization: consider a high-scale MSSM scenario where all SUSY particles have a degenerate mass MSUSY while only the SM Higgs remains light. A possible parametrisation may look like Note that the symbols MSUSY, At, Ab, Ae and TanBeta are not defined in the MSSM model file. Thus, additional information about these symbols must be provided using the Assumptions option, otherwise they are assumed to be arbitrary complex numbers. The initialization is invoked by There are a few important comments concerning the parametrisation which we have used in this example: -The symbol epsUV is used to indicate dimensionful parameters X which are to be neglected in the UV theory.
One should always use this parameter instead of the simpler rule X->0 to avoid problems caused by a division by 0. -It is recommended to express all matching conditions in terms of the running parameters of the effective theory, see Sect. 2.3. Therefore, we express the MSSM gauge and Yukawa couplings by the SM ones using the suffix Q which marks the running parameters (instead of g1 we e.g. specify it to be g1Q). For these parameters, only the tree-level matching conditions are required. The one-loop matching conditions for the gauge couplings, discussed in Sect. 2.3, are automatically derived. -Delta[a,b] is the SARAH internal symbol for the Kronecker delta δ ab . We use it here to include only contributions from third generation Yukawa couplings, and to force diagonal soft masses for the sfermions.
-In order to simplify the analytical calculation, we assume that all parameters are real. This is translated by conj[x_]->x. The object conj is the SARAH internal command for complex conjugation.
When all calculations are finished, it is possible to validate if the obtained mass spectrum at the matching scale is as expected In [5] ? M which yields the result As expected, the spectrum at the matching scale contains one massless CP-even Higgs boson which corresponds to the SM-like Higgs boson. Also all SM-like fermions remain massless while the heavy fields are degenerate in the mass parameter MSUSY.
The rotation matrices are stored in the array ReplacementRotationMatrices and read in our example Let us now continue with the description of the analytical interface. After the successful initialization and calculation of all mass and rotation matrices, one can compute the leading order (LO) and next-to leading order (NLO) corrections to an amplitude with the external fields given in the list fieldslist The fieldslist can contain two, three or four scalar fields including their generation indices to obtain effective mass parameters and cubic or quartic couplings. Note that the matching of effective mass parameters is only demanded if no spontaneous symmetry breaking occurs in this sector of the theory. The possible options for the function EFTcoupNLO are -  [2],hh [2],hh [2]).
-SimplifyResults -> $BOOL -Default: True -Description: whether to simplify the results using the given assumptions or not.
The full expression at the one-loop order is rather lengthy. Therefore, we make a few approximations and include only the terms involving the top quark Yukawa coupling. This can be achieved by setting all other couplings to zero. The command where we have introduced the stop mixing parameter X t = A t −μ tan −1 β, yields where the symbol UVscaleQ is the name for the renormalisation scale used in the loop functions. Note, because of the assumption g i → 0 this corresponds only to the leading oneloop shift but not to the full NLO expression (including the tree-level contributions), i.e. we found which is the well-known leading one-loop shift maximized for X t = √ 6M SUSY . Advanced Examples: the root directory of the new SARAH version includes the file Example_Matching.nb which contains already evaluated Cells that describe the example usage of all possible Options of EFTcoupNLO (e.g. the selection of specific topologies or debugging) within the high-scale MSSM.

Batch mode
The complexity of the calculation requires a high degree of reproducibility of the results. For this purpose it is possible to write input files that contain all necessary information for the matching to a given EFT model. This includes all information already discussed in the interactive mode. In addition, the correspondence between effective couplings in the lowenergy model and amplitudes in the UV model, as it was demonstrated for the THDM matching, have to be defined.
The batch mode is invoked during the initialisation by specifying the input file <FileName> located in the directory of the loaded SARAH model The mandatory content of the input file is Up to $NameUV and $MatchingConditions this is the same information which is otherwise passed to InitMatching and EFTcoupLO/EFTcoupNLO in the interactive mode. In addition, one can define options to control the generation of L A T E X or SPheno output. This is described in more detail below. First, consider an input file example which defines a high-scale SUSY scenario Matching_SimpleHighScaleSUSY.m ¤ 1 $NameUV="SimpleHighScaleSUSY"; .hh [1].hh [1].hh [1]  Here, we skipped most of the lines for $ParametrisationUV because they are similar to the definition of MyParametrisation in the last subsection. For simplicity, we set here all trilinear sfermion couplings as well as the matching scale UVscaleQ equal to MSUSY.
If the option InputFile->"Matching_SimpleHighScaleSUSY.m" is given to InitMatching, SARAH will calculate all matching conditions defined in $MatchingConditions. The information is stored in the arrays The meaning of these lines is -$EFTcouplingsToTeX: if set to True, all information obtained during the matching is exported into a L A T E X file ready to be compiled by standard L A T E X compilers. -$AdditionalTeXsymbols: a list containing replacement rules that define the correspondence between L A T E X and Mathematica expressions which are for instance used in the defined parametrisation. This will improve the readability of the L A T E X document significantly.
Thus, for our chosen example, the entries might read where the additional backslash is a necessary escape character. The L A T E X files are saved in the same directory $SARAH_Directory/Output/$Model/EWSB/Matching/ $NameUV as the other outputs.
SPhenooutput With little effort, it is also possible to generate a SPheno version which includes the analytical matching conditions to be used within an iterative running between the matching and the EW scale. In order to do so, two steps are necessary: 1. Export the Mathematica expressions into Fortran code and write a corresponding SPheno.m file 2. Run the EFT model using this SPheno.m The first step is again steered through the input file of InitMatching by adding the following information The export into SPheno routines is enabled with the first line. This option is sufficient to obtain Fortran routines for all matching conditions at the one-loop level. All other information must be given to automatically generate a suitable SPheno.m for the EFT model. Most variables have a 1:1 correspondence to the standard variables (without the $SPheno prefix) used in SPheno.m files discussed in "Appendix C". The new option is $SPhenoMatchingScale which defines at which scale the matching should be performed.
Running InitMatching with an input file containing these lines, produces two outputs: HighScaleSUSY.f90, located in the output directory of the MSSM model, which contains the matching conditions in Fortran format where most of the terms in the sum have been omitted as they are not important for the discussion. -a Mathematica file named SPhenoEFT_MSSM_ SimpleHighScaleSUSY.m which is located in the model directory of the SM. This file may look like One can see that this file contains the information given to InitMatching (line 1-12). In addition, the information about the matching and the corresponding Fortran routines (using parameter without the Q prefix) have been automatically added by SARAH (line 14-16 and 21-26).
The second step to generate a numerical code that includes the computed matching conditions is to run a new Mathematica kernel and call the SARAH routine MakeSPheno using the generated SPheno.m, i.e. This generates all necessary Fortran routines for the high-scale SUSY implementation. The code is compiled in the same way as other SARAH generated SPheno modules: The actual behaviour of the compiled SPheno code is described and compared with the fully numerical approach in the next section.

Matching at two scales
The analytical matching procedure discussed so far supports the derivation of effective scalar couplings from a high-scale theory at a single matching scale. Thus, towers of effective theories where the different sets of RGEs are needed between the different matching scales are not a priori possible in this approach. On the other side, SPheno always provides the possibility to perform a pole-mass matching between a given BSM model and the SM as described in detail in Ref. [25]. Thus, the functionality can be used to obtain precise prediction for scenarios like where large scale separations between the two BSM models as well as the SM exist. This is for instance the case for split-SUSY where the electroweakinos are in the multi-TeV range. Thus, such scenarios are already fully covered. An even more general implementation to allow for an arbitrary number of matching scales and an RGE running in-between is only possible with the numerical approach which we discuss next. A schematic overview about the numerical evaluation 4 SPheno can be downloaded from spheno.hepforge.org. of a parameter point when using the analytical calculation of matching conditions is shown in Fig. 5.

Numerical approach
The second option to generate a SPheno version for an effective model including the matching conditions to a UV theory is to set up a suitable SPheno.m for the EFT from the very beginning. This file must include the following information in addition to the standard information which is usually defined in the SPheno.m files, see "Appendix C": SPheno.m  where F1 and F2 are the involved fermions and S is the involved scalar. Yukawa-like interactions are chiral couplings. Therefore, the main difference to purely scalar couplings is the second argument containing PL/PR (for P L ,R = 1 2 (1 ± γ 5 )) to define which part of the coupling is meant. Moreover, the keyword ShiftCoupNLO can be used just to obtain the one-loop shift to a coupling, e.g. Several examples for the usage of these options are given below.

One matching scale without RGE running above
We start again with the simplest example of high-scale SUSY without any RGE running above the matching scale. Thus, the produced SPheno code will generate the same results as the one with the analytical approach in the last section. In order to set up a high-scale SUSY version with degenerate SUSY masses at the matching scale, the corresponding lines in the SPheno.m located in the model directory of the Note that for simple high-scale theories, without additional light fields, the SPheno.m could also be stored in the SM model directory as the two models are technically the same. The newly introduced models are described in Sect. 3.5.
The definitions are very similar to the analytical approach: the symbol epsUV again has been used to neglect specific parameters at the matching scale. An important difference is that we have not singled out the contributions from only third generation Yukawas because this would not give any performance improvement for the numerical calculation. Note, that it is also not necessary to define the matching for the Yukawas when running down. Moreover, we have used the option to define the scale where the RGE running should stop as function of an input parameter (UseParameterAsGUTscale={m0}). 5 Thus, SPheno will run the RGEs only to that scale and evaluate the SUSY boundary conditions.
The process to generate the SPheno output and to compile the Fortran code is identical to the final steps for the analytical approach: 1. Run MakeSPheno of SARAH with the new input file

Copy the files and compile SPheno
For a SPheno version generated in that way, two additional flags are available in the Les Houches input file to have some control over the calculations:

LesHouches.in
Thus, these flags can be used to: 201 Turn on/off all one-loop contributions to the matching.
By default, they are turned on. This might be helpful to check the size and importance of the one-loop corrections.
202 Turn on/off the contributions from the off-diagonal wave-function renormalisation. By default, they are turned off. See Sect. 2.2 for more details.

Running above the matching scale
We can modify the last example easily to include also the running above the matching scale. This might be for instance necessary if one wants to apply the SUSY boundary conditions at the scale where the gauge couplings do unify but not at the matching scale. In order to do so, one needs to remove UseParameterAsGUTscale = {m0} from the last example and put instead SPheno.m.GUT Thus, SPheno stops the running once the condition g 1 (Q) = g 2 (Q) is fulfilled.
In addition, the matching conditions for the Yukawas are changed to

SPheno.m.GUT
The need for the normalization onto the tree-level rotation matrix elements ZH is described in the next section. In that way, we can include the one-loop shifts to all Yukawa couplings which are necessary to have a consistent RGE running with two-loop SUSY RGEs between the matching and GUT scale, see also the discussion in Sect. 2.4. Note, we did not consider any generation indices for the involved fermions, i.e. the result of ShiftCoupNLO is a 3 × 3 matrix. If one wants to safe program run-time it is possible to consider the one-loop shifts to the top Yukawa couplings only.

SPheno.m.GUT
Moreover, the shifts for the gauge couplings are applied automatically.

Several matching scales
With the above settings one can now implement an arbitrary number of matching scales. However, as we have noted already in Sect. 3.3.3, the pole-mass matching to the SM is automatically included in the SPheno output. Thus, if a second matching scale, which is not too far away from the EW scale, is needed, one can simply rely on that. However, if more than two matching scales are needed, or if the matching to the SM should take place at such a high scale where the pole-mass matching might suffer from numerical problems, 6 one can now start to build up towers of EFTs by defining more matching scales in SPheno.m. For instance, the full input to define the tower SM → THDM → THDM + electroweakinos → MSSM is given in "Appendix D". In this example we also make use of the functionality to calculate new fermionic couplings at the one-loop level below a matching scale: Here, g u,d 1,2 are the split-SUSY couplings between the Higgs boson and a Higgsino-Gaugino pair, see e.g. Ref. [10]. We include these corrections by considering the one-loop amplitude between the Higgs boson and a pair of neutralinos. In this example we have also used another feature: we have not explicitly defined the generation indices of the involved neutralinos. The reason for this is: even if the neutralino mass matrix contains only zero's under the given approximations (μ, M i M SUSY ), it is not clear how the mass eigenstates are ordered in the numerical run. Therefore, we have used the name of the gauge eigenstates. By doing that, SPheno checks during the numerical evaluation which of the mass eigenstates has the biggest contribution of the given gauge eigenstate. Of course, if the rotation matrix for the neutralinos is not equivalent to the unit matrix, i.e. if some mixing appears for instance because of effects of non-vanishing μ, one needs to define Thus, the rotation to mass eigenstates, which should take place just at the weak scale, is divided out.

Summary
A summary of the numerical evaluation of a parameter point with SPheno which includes several matching scales M M n and optionally also the running to the GUT scale M G is given in Fig. 6.

Included models and input files in SARAH
Several models which make use of the new functionality have already been implemented and are part of the publicly available SARAH version. All hierarchies considered for the MSSM so far are summarised in Fig. 7. Also for the NMSSM with very heavy particles two models exist: the high-scale NMSSM, where all SUSY fields are integrated out and a split-NMSSM, where the singlet and the SUSY fermions are kept.
The names of the new models that make use of the numerical approach are listed in Table 1. Also for the analytical approach several input files are now included in SARAH. Those are summarised in Table 2. Based on these examples and by the explanations in this section, it is now straightforward for the users to implement their own scenarios.

Examples, self-consistency checks and comparisons with other codes
The following section describes realistic examples of practical applications of the presented framework. We consider different high-scale SUSY scenarios which were already studied intensively in literature. In particular comparisons between predictions for the SM Higgs boson mass derived with our generic setup against dedicated tools and calculations are made. In this context, we demonstrate also the perfect agreement between the two available options to use SARAH/SPheno for numerical studies. Finally, we also show that one can easily obtain precise results for other highscale extensions for which no other tool existed so far.

Low-energy limits of the MSSM
In the introduction it was already mentioned that SUSY models with a SUSY breaking scale well above the electroweak scale became more popular in the recent years. While in these scenarios the direct observation of SUSY states is difficult or even impossible, these models are severely constrained by the Higgs boson mass measurements. For instance, if the masses of all superpartners are degenerate, the highest possible SUSY breaking scale in the MSSM is about 10 10 GeV [10]. For higher scales, the predicted m h always becomes too large. Since the Higgs boson mass in these models is the crucial observable, a precise calculation is mandatory and specialised codes have been developed to get reliable predictions. We are going to consider three different cases: (i) split-SUSY in which all SUSY scalars are very heavy, but electroweakinos might stay moderately light, (ii) high-scale SUSY in which all SUSY masses and the additional Higgs with Here, m 2f are the soft masses squared for all chiral superfields, M A is the mass of the heavy Higgs doublet, M i are the soft gaugino masses, μ is the Higgsino mass term in the superpotential, and B μ , T i are the soft-breaking equivalents of the μ-term and the Yukawa couplings in the superpotential.

Split-SUSY: MSSM → SM & electroweakinos & gluinos
Split-SUSY with very heavy SUSY scalars but significantly lighter SUSY fermions keeps most of the nice SUSY properties like gauge coupling unification and provides a viable dark matter candidate. In this setup, the full MSSM is matched to the SM extended by additional fermions. The Lagrangian of  Table 1 The names of the new models which are part of SARAH 4.14.0. The hierarchy in the last column refers to Fig. 7. For the split-NMSSM also a light singlet is present, i.e. the hierarchy is similar to (c), but not identical where the Yukawa couplings g u,d 1,2 are as in the example of Sect. 3.4.3, σ a are the pauli matrices and α = 1, ..., 8. In order to calculate the Higgs boson mass in this model, the common approach is to (i) decouple the SUSY scalars at the scale M SUSY and calculate λ SM (M SUSY ) including important higher-order corrections, (ii) run the split-SUSY RGEs to the scale M χ of the remaining SUSY states and calculate the shift in λ SM (M χ ), (iii) run the SM RGEs to m t and calculate m h (m t ) at the two-loop level. The full results for the oneloop matching conditions at M SUSY and M χ were given in Ref. [10]. Also the dominant two-loop corrections to λ SM of For the case of electroweakino masses of 1 TeV we show also the SPheno result when using a two-loop fixed-order calculation in the EFT. We see, in agreement with a previous study in Ref. [64], that the two-loop contributions of the additional fermions have only a mild effect on the SM-like Higgs boson mass.

High-scale SUSY: MSSM → SM
An even more extreme setup than split-SUSY is high-scale SUSY in which all SUSY partners are very heavy. Thus, the effective model is just the SM, i.e. Footnote 7 continued 21 2g 2 2dg 2u (g 2 1d +g 2 1u ) misses one power ofg 2u . We fixed that and in all following results the patched version of SusyHD is used.
and the only visible impact of SUSY is the prediction of λ SM at the matching scale M SUSY M χ M A . The matching conditions at the SUSY scale are just the combination of the two matching conditions for split-SUSY applied at a single matching. Thus, it is obvious that also for this case a full agreement between our analytical results and those of Ref. [10] exists. However, Ref. [10] includes also the dominant two-loop corrections in the case of high-scale SUSY which also entered the code SusyHD. Therefore, it's worth to discuss also the numerical differences between SPheno and SusyHD for the case of high-scale SUSY. The model implemenation in SARAH is called

HighScaleSUSY/MSSM
The results are summarised in Fig. 9. In addition to the comparison to SusyHD we also compare the results to two other calculations: a standard fixed-order calculation as well as an EFT calculation based on the pole-mass matching [25]. In the pole-mass matching, the quartic coupling λ SM is calculated from the condition which can be translated into where Π SM h are the loop corrections to m h known from the SM. The pole-mass matching has the advantage that also terms v M SUSY are included and that only two-point functions need to be calculated instead of four-point functions, see Ref. [24] for more details. On the other side, this approach has also some drawbacks. It is mainly restricted to the SM as EFT, but it is not straightforward to be used in models with several light scalars. Also a consistent matching at the two-loop level needs some fiddling with the running parameters which enter the different parts of Eq. (38), see Ref. [70]. While SPheno by default used MS parameters to calculate Π SM  Fig. 9. The difference between both results is a two-loop effect and could be taken as estimate of the remaining uncertainty in the one-loop pole-mass matching.
Moreover, we find that the pole-mass matching becomes also numerically unstable -at least in SPheno-once M SUSY v is used because the loop functions used for the pole-mass calculations are not optimised for these cases: we see in Fig. 9 that the pole-mass matching breaks down at M SUSY 5 · 10 5 GeV. Nevertheless, we find that the agreement between the pole-mass matching and the direct matching procedure presented here is very good for SUSY scales up to 100 TeV. One finds also that the fixed-order calculation agrees perfectly with the pole-mass matching for M SUSY below 1 TeV. Of course, for larger SUSY scales, the discrepancy between the fixed-order calculation and all EFT calculations grows very rapidly.
We come back to the comparison with SusyHD: we see that the agreement between SPheno and SusyHD is also very good and the differences are always of the level of 1 GeV or below. The 1 GeV differences appear only for the choice A 0 = 2M SUSY and M SUSY around the TeV scale. In that case, the two-loop corrections missing in SPheno play some role. However, for larger M SUSY or smaller trilinear terms, these two-loop corrections cause only a moderate shift -or become even completely negligible. Thus, we think that it is not a substantial drawback of our setup that 'only' one-loop corrections are included so far.

High-scale SUSY with intermediate M A : M SSM → THDM
In the last example we have assumed that all BSM particles are very heavy and degenerate. An important deviation of this ansatz is the possibility that the second Higgs doublet remains light, i.e. only fields with negative R-parity are very heavy. In this case, the low-energy theory of the MSSM is a Two-Higgs-Doublet-Model type-II. 8 The Lagrangian of the EFT is 8 Strictly speaking, one obtains a THDM type-III when integrating out all SUSY fields in the MSSM because the 'wrong' Yukawa couplings ∼ H * dq u are loop-induced. However, this becomes mainly important for flavour violating observables and has no visible impact on our discussion of the Higgs boson mass prediction here.
One can make the following association between fields at the SUSY scale to calculate the matching conditions However, this choice is not unique as there is no preferred basis of Higgs doublets in a general THDM, i.e. one could also interchange H 1 and H 2 or take any linear combination of them. With the common choice made in Eq. (40), one can simultaneously apply a rotation into the mass basis on (H 1 , H 2 ) and (H u , −iσ 2 H * d ) so that the tree-level mixing angle tan β of the MSSM coincides with the effective THDM.
The dominant threshold corrections to λ 1 -λ 7 involving third generation Yukawa couplings are available in literature [29]. We have double checked the analytical expressions derived by SARAH and found full agreement.
The importance of the proper matching to the THDM for the case M A M SUSY has been pointed out in Ref. [36]. It was found that in particular for small tan β very large difference to a one-scale matching appear. In order to demonstrate that, we compare in Fig. 10  Since we have demonstrated the importance of performing the matching to the THDM properly for the case of a light second Higgs doublet, it is clear that codes were developed to include these effects. The first tool in this direction was MhEFT which uses a purely EFT ansatz [36]. In a recent update of FeynHiggs a hybrid ansatz combining the fixedorder calculation with higher-order terms was implemented [40].The overall agreement between both codes turned out to be good once a careful translation between the parameters in both renormalisation schemes was done. Since MhEFT is much closer to the ansatz of SARAH/SPheno we are going to compare our results with this tool. 9 For this purpose, we have set up the model  One can also see that the impact of the additional twoloop corrections implemented in MhEFT is very moderate. Thus the main source of the difference is the determination of the running top Yukawa coupling. In contrast, the additional one-loop corrections due to gauginos, which were presented very recently also in Ref. [40], can be easily included in SPheno using the numerical matching interface. For the considered choice of parameters these have numerically a bigger effect than the two-loop corrections and cause a shift of 1-1.5 GeV.

High-scale NMSSM
Up to now we have only discussed examples of models involving very heavy BSM particles which could already be studied with public tools like SusyHD, MhEFT, FlexibleSUSY or FeynHiggs. These are just different low-energy limits of the MSSM. However, our framework is not restricted to this case and in principle any SUSY or non-SUSY model could be considered as high-scale theory. We show that crucial differences compared to the MSSM show already up in the case of the NMSSM. The NMSSM involves an additional gauge singlet superfieldŜ which leads to the following superpotential after imposing a Z 3 symmetry to forbid all dimensionful parameters where W Y represents the terms involving Yukawa couplings that are identical to the MSSM. The NMSSM-specific soft-SUSY breaking parameters are The scalar singlet S can receive a VEV even without EWSB which causes an effective Higgsino mass term We can now study what the impact of the additional gauge singlet in a high-scale SUSY scenario is. For this purpose, we impose the following relation among the parameters The full high-scale model has three free parameters The MSSM limit is obtained for λ S → 0. We have implemented this model in SARAH as

HighScaleSUSY/NMSSM
The predicted mass for the SM-like Higgs boson as function of the SUSY scale M SUSY is shown in Fig. 12  Of course, one could now start to consider also other lowenergy limits of the NMSSM. However, this is beyond the scope of this paper here and interesting applications are given elsewhere [71].

Summary
We have presented an extension of the Mathematica package SARAH which derives the one-loop matching conditions for effective scalar couplings based on a UV theory. Two different approaches exists, which are based on either an analytical or fully numerical calculation. The full agreement between both calculations and analytical results available in literature has been pointed out. Furthermore, good agreement with specialised codes to study Split-or High-scale SUSY like SusyHD or MhEFT was shown. Since our approach is completely general, it can be used to study UV completions of a large variety of BSM models with and without an extended Higgs sector.

A: Generic diagrams
In this appendix we provide a complete list of all possible one-loop diagrams with 2, 3 and 4 external scalars and internal fermions or scalars. The results were obtained in the limit of vanishing external momenta using the computer programs FeynArts and FormCalc [72,73].
We distinguish between topologies, where neither the statistical nature (spin=0,1, 1 /2) nor the mass (light or heavy i.e. zero or non-zero) of internal fields is specified and generic diagrams, where the spin of all fields is specified but not their mass hierarchies. This is still a model-independent graph but gives the possibility to write down a model-independent generic amplitude. Couplings and masses appearing in the expression of a generic amplitude are seen as generic couplings that do not have any relation to parameters used in other generic amplitudes.
As already mentioned in Sect. 3, SARAH's analytical matching interface does not only compute the full one-loop contribution but is also able to compute a subset of diagrams based on e.g. a choice of topologies. For this rea-son, a notation with a successive structure is introduced. A topology is described by a string consisting of maximum four characters. It starts with the specification of the diagram type which can be tree-level (T), self-energy (S), WFR (W) or ordinary one-loop diagram (blank) followed by a letter specifying the type of the involved loop integral defined in Appendix B. In the next digit, the number of external fields is specified. If the diagram contains an internal single propagator additional to a loop (2) or a loop only (1) is denoted by the next digit (blank means that no diagrams with additional internal lines exist). The last digit is a counting index (blank means only one diagram of that type exists). An example explanation of the notation is given in Fig. 13 All possible topologies and groups are stored in the list TopoNotation. In the following, we list the analytical expressions for all generic amplitudes as well as the topologies they belong to in the format used in the SARAH matching routines.

A.1.1: Quartic couplings
There are two tree-level topologies with four external scalars. The first one is a local quartic coupling which could for example be given by supersymmetric D-terms and/or F-terms while the second one has one internal propagator, necessarily heavy and of bosonic nature. Thus, fermions can only enter one-loop and higher-order corrections.

A.1.2: Cubic couplings
Since there is only one three-point topology, the matching of trilinear couplings at tree level is trivial.

A.1.3: Bilinear parameters
The two-point function is necessary for the matching of scalar sectors that involve non-Higgs scalar fields i.e. scalars that do not develop a VEV. In this case the scalar masses and couplings are independent parameters and have to be matched separately. (A.4)

A.2: One-loop contributions
At the one-loop order, we distinguish between irreducible diagrams and reducible diagrams which contain an additional internal propagator line.

A.2.1: Irreducible diagrams
Quartic couplings Since we consider renormalisable theories, fermions can only enter the one-loop corrections to quartic couplings in box diagrams (topology D, note that D [4] is not a valid topology as it is a reserved Mathematica symbol). The generic diagrams and amplitudes are

A.2.3: Wave-function corrections
Contributions on external legs are divided into diagonal and off-diagonal WFR topologies named DiagonalWFRs and OffdiagonalWFRs. They consist of a Z-factor times the treelevel amplitude times a combinatorial factor. In the following, we give the analytical expressions for the Z factors. with A 0 (m 2 1 ) + A 0 (m 2 2 ) + (m 2 1 + m 2 2 )B 0 (m 2 1 , m 2 2 ) + m 1 m 2 (c L 1 c L 2 + c R 1 c R 2 )B 0 (m 2 1 , m 2 2 ) , (A.31) where the dotted notation is introduced in "Appendix B".

B: IR-safe loop-functions
In this appendix we give analytical expressions for all loop functions used in the matching routines. In particular, we list the limits for all possible combinations of vanishing and equal masses as they are needed to provide numerical stability of the matching routines. The common prefactor (B.34) The integrand I n is symmetric w.r.t to the masses and thus also the loop functions are symmetric w.r.t. their arguments.
One can reduce all integrals to the A 0 integral by using partial fractioning and integration by parts. We define the following abbreviations for finite logarithmic terms The content of the SPheno specific input file for SARAH is the following: 1. Input parameter (MINPAR, EXTPAR): a list of parameters which should be read by SPheno from the block MINPAR or EXTPAR in a LesHouches file. Note that there are no hard coded entries for MINPAR or EXTPAR. This makes it necessary to define these blocks also for models with already existing SLHA conventions. However, this also provides more freedom in varying the model and the free parameters. 2. RealParameters: By default, all parameters defined in MINPAR or EXTPAR are assumed to be complex, i.e. it is possible to use also the block IMMINPAR to define the imaginary part. However, some Fortran functions like sin cannot be used with complex numbers, therefore is is necessary to define parameters like tan β explicitly as real.