N-loop running should be combined with N-loop matching

We investigate the high-scale behaviour of Higgs sectors beyond the Standard Model, pointing out that the proper matching of the quartic couplings before applying the renormalisation group equations (RGEs) is of crucial importance for reliable predictions at larger energy scales. In particular, the common practice of leading-order parameters in the RGE evolution is insufficient to make precise statements on a given model's UV behaviour, typically resulting in uncertainties of many orders of magnitude. We argue that, before applying N-loop RGEs, a matching should even be performed at N-loop order in contrast to common lore. We show both analytical and numerical results where the impact is sizeable for three minimal extensions of the Standard Model: a singlet extension, a second Higgs doublet and finally vector-like quarks. We highlight that the known two-loop RGEs tend to moderate the running of their one-loop counterparts, typically delaying the appearance of Landau poles. For the addition of vector-like quarks we show that the complete two-loop matching and RGE evolution hints at a stabilisation of the electroweak vacuum at high energies, in contrast to results in the literature.


I. INTRODUCTION
Minimal extensions of the Standard Model (SM) are invaluable tools in the pursuit of physics beyond the SM (BSM). They offer the possibility of studying different effects at energy scales testable by the Large Hadron Collider (LHC) in a comparably clean environmenti.e. the models typically contain the minimal numbers of new fields to exhibit novel phenomenology. All additional states also alter the high-scale behaviour of the model compared to the SM expectations. For instance, it is known today that the SM becomes metastable if it is extrapolated to very high energies [1][2][3][4]: at a scale of 10 9−11 GeV the quartic coupling runs negative. The scale at which the potential becomes unstable could be significantly affected by the presence of new states -it could even completely disappear. This would then indicate that the BSM model is valid up to the Planck scale. An opposite effect can occur if large couplings are present. In this case, a Landau pole might be present which points towards the breakdown of the theory. Both effects, the presence of Landau poles or deeper vacua, can also be used to directly constrain the parameters of a new physics model. A parameter point can be discarded if the model becomes strongly interacting at energies already * braathen@lpthe.jussieu.fr † goodsell@lpthe.jussieu.fr ‡ mkrauss@th.physik.uni-bonn.de § toopferk@uni-mainz.de ¶ florian.staub@kit.edu probed by the LHC, or if the life-time of the electroweak breaking vacuum is too short on cosmological time scales. Many of these effects have already been studied in the literature for plethora of different models such as singlet extensions [5][6][7][8][9], triplet extensions [10,11], Two-Higgs-Doublet-Models (THDMs) [12][13][14][15][16][17][18][19][20][21] or models with vectorlike fermions [22]. These studies utilise the one-and sometimes even two-loop renormalisation group equations (RGEs). However, less care was was taken in the determination of the parameters which enter the RGE running. Often, two-loop RGEs were combined with a tree-level matching.
A proper determination including higher-order corrections of the quartic coupling, which enters the RGE running, was so far only performed for the SM [1]. 1 It was shown that even the next-to-next-to-leading-order (NNLO) shifts to λ are important for determining the fate of the model. This is remarkable, because it is well known that the loop corrections to the Higgs mass are small if they are calculated at Q = m t : the corresponding shifts in λ are only 2.5%. While the corrections from top quarks are of a similar order in many BSM models, other corrections like the ones from Higgs self-interactions can be much larger. This has been recently pointed out in Ref. [27] in the example of THDMs where the two-loop 1 Loop corrections in the scalar sector were taken into account in Ref. [23] for a singlet extension and in Refs. [24,25] for a THDM. These studies did not however investigate the impact on the high-scale behaviour of the model. In Ref. [26] in turn, a one-loop matching has been performed for that purpose in the context of a seesaw-II as well as a left-right symmetric model.
corrections to the Higgs masses were calculated for the first time. 2 In this context it is important to note that the common lore "N -loop running needs N − 1-loop matching" is applicable only in certain scenarios. This degree of matching catches the leading logarithms correctly but misses finite contributions which can be relevant in many BSM applications. We show in this work that higher order corrections can be very important for the study of the UV behaviour of a theory leading to four main conclusions: 1. The threshold corrections at low energy can lead to substantial (finite) shifts in the running parameters of a model: therefore N -loop RGEs with N -loop matching is required for consistency.
2. The change from one-loop to two-loop running can flatten the running at large values of the coupling, preventing the onset of a Landau pole at high energies -leading to a form of asymptotic safety.
3. Alternatively, in the case where the running drives some quartic coupling negative, higher-order corrections can lead to significant changes to the predicted scale of metastability. 4. As a by-product of the above, we find that new fermionic fields at low energies can stabilise the SM potential.
We illustrate the above with a detailed examination of three examples: a singlet extension, the SM extended by vector-like quarks and the THDM. The paper is organised at follows. In Sec. II we give a step-by-step prescription for the general matching procedure including effects at the N − 1and N -loop levels. In Sec. III we give details into the procedure used to obtain higher order corrections to the quartic couplings in the different models considered, before we discuss in Sec. IV the numerical results, providing insights including approximate formulae.

II. MATCHING AND RUNNING
To extrapolate a theory from the electroweak scale to high energies, we require two ingredients: 1. The value of the couplings at the "low scale" where the running starts.
2. The RGE running of all parameters. 2 In this context, it was pointed out that using "on-shell" masses and couplings as input can be quite dangerous because it hides the presence of large couplings which could even spoil perturbation theory. A similar observation was made for another model, the Georgi-Machacek model, in Ref. [28].

A. Renormalisation Group Equations
We shall always work in the MS scheme. In this scheme, the β-functions, which describe the energy dependence of the parameters Θ, are defined as Here, µ is an arbitrary mass scale. β i can be expanded in a perturbative series: are the one-and two-loop contributions to the running which we are mainly interested in. The expressions for the two-loop running of the parameters appearing in a given model can be obtained from the generic expressions valid for a general quantum field theory, given in Refs. [29][30][31][32].

B. Matching
The renormalised coupling constants Θ i in d = 4 − 2 dimensions, which enter the running, are related to the corresponding bare couplings Θ 0 i by Here, C i are constant factors depending on the character of Θ i 3 . The coefficients a i are the result of a perturbative expansion. In general, two approaches are possible to determine the MS parameters as function of physical observables such as masses.
1. In an on-shell calculation the physical observables are identical at each loop-level, but all finite and infinite corrections are absorbed into the counterterms of the Lagrangian parameters (δΘ OS i ).
2. In an MS calculation the counter-terms of the Lagrangian parameters (δΘ MS i ) contain only the divergences. Therefore, the calculated masses depend on the loop-level at which the calculation is performed.
The bare Lagrangian parameters are identical in both cases In an on-shell calculation, however, there is no generic set of renormalisation group equations known, and therefore to explore a theory at high energies it is necessary to use MS equations -i.e. to extract the underlying MS parameters of the theory, and then run them. On the other hand, in an MS calculation, the physical parameters are functions of the Lagrangian parameters: so if we are given the physical quantities, we must invert these functions to extract the MS ones. This is where complications appear, and why many studies in the literature resort to simply using tree-level matching.
For example, suppose that we want to extract the quartic coupling of the SM from the Higgs mass. The Higgs mass M h is, however, calculated in terms of the underlying Lagrangian parameters as a loop expansion via the on-shell condition This is in general a highly non-linear equation in λ; but fortunately since the series is perturbative we can solve it through expanding to find , which is simple enough for the Standard Model and extensions without scalar mixing -so we shall give analytical expressions in Secs. IV B 1 and IV C. On the other hand, for more complicated models, we need to solve Eq. (5) through iteration, and we shall adopt this approach in general for the numerical studies. Whichever way we solve for λ, as we shall argue in the next subsection, when we are using N -loop RGEs, to obtain a consistent expansion of our MS parameters in terms of the physical ones we should use N -loop matching.

C. Loop Level of Matching and Running
It is commonly accepted in matching a high-energy theory onto an effective field theory that, if the running is performed using N -loop RGEs, it is only necessary to use N − 1-loop threshold corrections. However, the rationale for this criterion is less well known: it corresponds to matching the logarithmic terms of the calculations, and neglects the finite parts. To take perhaps the simplest and best-known example, suppose that we want to integrate out two heavy Dirac fermions with masses M 2 > M 1 that couple with charges Q 2 , Q 1 to some U (1) gauge theory. Suppose that the contribution of other fields to the one-loop beta function is b 0 , so at high energies the beta function is b 0 + 4 3 (Q 2 1 + Q 2 2 ), and we take the gauge coupling at high energies to be g(Λ). Now, to determine the gauge coupling at a low energy µ the classic prescription is to run first to M 2 , match at tree level, then run to M 1 , match again, and then run down to µ. This gives the one-loop value of the gauge coupling to be This can be rewritten as which is also the result of simply matching the two theories at the scale µ and including the threshold corrections: for corrections to the gauge coupling, the thresholds contain only logarithmic corrections. However, suppose that instead we wanted to match the two theories at a scale M that we have chosen to be different to M 1 and/or M 2for example, because other sectors of the theory contain particles of that mass. Now, if we just match the two theories at tree level, we would find after first running to M and then matching: In this case, we have a difference between the two values of the gauge coupling at one loop: in other words, we conclude that for single-scale matching, even in this simple case there is a discrepancy of one-loop order between the two calculations. What we have captured is, as mentioned above, the logarithmic terms only. The reader might complain that this is a slightly strange example since the threshold corrections themselves contain only logarithmic terms. We can instead take another example, more relevant for this work, of the quartic coupling of a real scalar φ coupled to a Majorana (Weyl) fermion ψ: If y and m are real, then the one-loop threshold correction to λ from integrating out the fermion ψ at matching scale M is The finite, non-logarithmic correction is highlighted in bold. On the other hand, Hence the tree-level matching clearly finds the correct logarithmic terms, but we miss the finite correction of − 2y 2 π 2 , even at one loop. One argument against using N -loop matching has been that the logarithmic corrections should be the most important: if the couplings are all small, then e.g. 2y 2 π 2 should be a very small correction and we will only have significant contributions to the couplings when we run down from high energies from the terms enhanced by logarithms. However, this preconception is biased from the idea of the SM where the strong gauge coupling and top Yukawa run to smaller values at high energies, and the quartic runs to zero; this also applies for models with no new scalars such as Split Supersymmetry, so for example we see that the supersymmetric threshold corrections become smaller as we take the supersymmetry scale higher -when the theory in the ultraviolet is the Minimal Supersymmetric SM.
For a generic BSM theory -especially ones with new quartic scalar couplings -this expectation is no longer true. Indeed, if we neglect other iterations then it is well known that a quartic coupling will tend to increase as we increase the energy -see e.g. Eq. (12) with y = 0. For a generic theory with new quartics, as we increase their values relative to the Yukawa or gauge couplings we expect for some given value there will be a transition from decreasing couplings in the UV (leading to metastability or instability of the potential) to the appearance of a Landau pole. In this latter case, since the quartic couplings will be even larger at high energies, including the finite parts of threshold corrections becomes vital for a consistent matching.
So far we have discussed threshold corrections when integrating out a heavy theory at high energies. This should not be confused with the matching that we need to do when we are running up to investigate the appearance of a Landau pole or the vacuum stability of a model. When we are applying threshold corrections at low energies (around the electroweak scale), then if we want to investigate whether a theory is stable, metastable or has a cut-off, and if so at what scale, then clearly it is important that the starting point for our calculation is determined accurately. For example, taking our toy model above and neglecting the coupling y, then integrating up from M to a scale µ we find a Landau pole at approximately but if we have incorrectly determined λ(M ) by an amount δλ, then the ratio of the correct cut-off scale to the erro- In the Standard Model, the difference between the treelevel value for λ and the two-loop value is tiny when the extraction is performed at the top mass. However, as we shall see, in other models a shift of 10% in the quartic coupling is conservative (and we should not forget the famous example of the MSSM where δλ λ).

III. NUMERICAL SET-UP
For the numerical calculations we make use of the Mathematica package SARAH [33][34][35][36][37][38] to produce a spectrum generator based on SPheno [39][40][41]. SPheno includes routines to obtain the full one-loop corrections to all masses as well as the two-loop corrections to real scalars. The two-loop calculations are done in the gauge-less limit and based on the generic results of Refs. [42,43]. In non-SUSY BSM models, these results suffer in general from the so called "Goldstone boson catastrophe" even in the gauge-less limit because the couplings of the Goldstone do not disappear in this limit, but are proportional to the cubic and quartic potential parameters. Therefore, we also make use of the results of Ref. [44] which provides a general solution to this problem.
In practice, we perform the following steps to calculate the mass spectrum based on a set of MS parameters: 1. The running couplings Θ(Q) at the scale Q = m t are taken as input, while the SM parameters are evolved to this scale including all known SM corrections, i.e. three-loop running and two-loop matching for strong coupling g 3 and top Yukawa Y t .

2.
The tree-level tadpoles T i are solved to fix the remaining free parameters, which in what follows are typically the mass parameters µ 2 i |φ| 2 . 3. The tree-level mass are calculated by diagonalising the tree-level mass matrices 4. The n-loop corrections to the tadpoles δ (n) t i are calculated. The imposed renormalisation conditions are which cause shifts in µ 2 i : 5. The one-and two-loop self-energies for real scalars are calculated for external gauge eigenstates. The pole masses are the eigenvalues of the loopcorrected mass matrix calculated as Here,M φ is the tree-level mass matrix including the shifts (16).
The calculation of the one-loop self-energies in both cases is done iteratively for each eigenvalue i until the on-shell condition is fulfilled. The renormalised rotation matrix is taken to be the one calculated for p 2 = m 2 φ1 .
If a chosen set of input parameters Θ(Q) results in the desired physical masses and mixing angles when using a N -loop calculation, we refer to them as N -loop couplings. Thus, with tree-level relations we have leadingorder (LO) parameters, while the one-and two-loop mass corrections result in the next-to-leading-order (NLO) and NNLO couplings, respectively.
Finding the correct set of MS couplings corresponding to the desired physical parameters at loop level is non-trivial. In what follows we use a simple fitting algorithm which varies the input parameters until the desired masses and mixing angles are obtained.

A. Singlet Extension
We start with the SM extended by a real gauge singlet S. The potential reads After electroweak symmetry breaking the CP-even scalar components mix to two physical states h, H via a rotation angle α. At tree level we can use m h , m H and t α ≡ tan α as input to calculate the quartic couplings Aside from the typical requirement that the quartic couplings remain perturbative, the constraints from perturbative unitarity need to be taken into account. For that, we can evaluate the scalar 2 → 2 scattering amplitude in the limit of high energies and demand that the eigenvalues stay below 8 π. The conditions for the model under consideration read which leads to λ S < 4π/3 for small λ and λ SH . Al- (1, 1) 0.33 0.24 -0.97 6.4 · 10 8 3.2 · 10 8 (1, 2) 1.3 · 10 12 2.5 · 10 9 (2, 1) 0.32 0.18 -0.94 5.1 · 10 10 2.5 · 10 10 (2, 2) 1.0 · 10 14 2.0 · 10 11    though a seemingly weak constraint at first sight, this can become a severe constraint particularly in the case of small v S , cf. Eq. (21). In the following, we compare different approaches for the matching of the quartic couplings. In Fig. 1, we show the values of the quartic couplings at the scale 10 6 GeV as a function of the scale Q where the matching is performed. This is done by first running the SM RGEs to the scale Q where we then match the quartic couplings to the spectrum at tree level, one loop as well as two loop. The final step is the running of the singlet extended SM RGEs (both at both one-or two-loop order) up to 10 6 GeV. In the left-hand plots, we use one-loop RGEs for the cases of tree-level matching (dotted) and one-loop matching (dashed). Although there is no dependence of the quartics on the matching scale when using tree-level matching, the scale dependence induced by the RGEs is larger than for the case of one-loop matching and oneloop RGEs. In the right-hand panels, we show the quartics at 10 6 GeV, evaluated with two-loop RGEs, when using one-loop (dashed) and two-loop matching (solid lines). Once again, the scale dependence is decreased when using N -loop matching with N -loop RGEs, where in this case N = 2.
Also visible in Fig. 1 are the large differences between the eventual coupling values when using the "traditional" approach of tree-level matching with one-loop RGEs and the approach we advertise, two-loop matching and twoloop RGEs. This also means that large differences are expected when evaluating the cut-off scale of a theory, i.e. the scale at which the model becomes non-perturbative or violates unitarity. In Tab. I, we show the cut-off scale of a particular parameter point when using n-loop matching in conjunction with m-loop RGEs. We show the scale at which the quartics become non-perturbative (Λ 4π ) separately to the case where either perturbativity or unitarity is violated (Λ unit. 4π ). The corresponding running of the individual couplings is displayed in Fig. 2. Note that we only display the cases (n, m) = (T, 1), (1, 1), (1, 2) and (2, 2) in this figure as dotted, dot-dashed, dashed and full lines corresponding to the comparison of N versus N − 1 matching. The impact of the two-loop RGEs is a moderation of the one-loop RGEs: while the one-loop β-function of λ S is given by β , so that λ S tends to grow very rapidly, there is a moderating term from the two-loop RGEs which goes with . Therefore, using the oneloop RGEs only, λ S grows large very quickly -whereas the unitarity limit is reached at a much later scale when using two-loop RGEs. Nevertheless, a complete stalling of the evolution is typically only reached at λ S values which already violate the unitarity limit according to Eq. (25), see the black dashed (full) line between 10 9 and 10 13 GeV (10 12 and 10 15 GeV). The moderation of the evolution of λ SH and λ is not as pronounced. For λ SH , the corresponding β-function grows with 1 16π 2 12λ S λ SH with only a small moderating effect from the two-loop RGEs. As a consequence, it becomes larger than 4 π before λ S and then drags the latter with it. In total, in particular because of the large two-loop contributions to β λ S , there can be several orders of magnitude between the eventual cut-off scales when using one-or two-loop RGEs.
The effect of using a two-loop matching instead of a tree matching, in turn, is a reduction of the quartic couplings. The reason is the positive mass corrections to m H , leading to smaller MS couplings when doing the proper loop-level matching. As shown here, the impact can be large and we observe positive shifts in the eventual cut-off scale by several orders of magnitude when including the matching. Fig. 3 we show in the m H -v S plane the differences between using N − 1-loop and N -loop matching when applying N -loop RGE running. The cut-off scale here and in what follows is defined as the scale at which either one of the couplings grows larger than 4π or any of the conditions for perturbative unitarity are violated, each evaluated with the running MS quartic couplings. The grey contours in the left-hand pane of Fig. 3 display the ratio of the evaluated cut-off scales for N = 1.

Finally in
In particular for small v S , which leads to large quartic couplings, the effects are quite drastic as loop effects become very important. The differences between one-and two-loop matching (shown as blue coloured contours) are significantly milder in this region, the maximum difference is just a factor of three. For large v S , instead, the quartic couplings are comparably small, leading to large cut-off scales in general. This also means, however, that during the long RGE running, small shifts in couplings can lead to more drastic effects as is seen in the upper region in the plot with v S 350 GeV. However, the cutoff scale differences stay below an order of magnitude for N = 2.
On the right-hand side of Fig. 3 we present the difference in cut-off scales between the most extreme cases, tree-level matching using one-loop RGEs versus two-loop matching using two-loop RGEs. In particular for small values of the singlet VEV, the eventual cut-off scale can be many orders of magnitude larger than the cut-off scale evaluated with tree-level matching. Grey contours show the ratios of the singlet quartic couplings at Q = m t between the two matching approaches, λ

Analytical approximation
We may now make a further simplification to the singlet extension studied in Sec. IV A, namely adding an additional Z 2 symmetry under which the singlet scalar is charged -for clarity we will call this model the Z 2 SSM to distinguish it from the SSM. This symmetry forbids nonzero values for the couplings κ 1 , κ 2 and for the singlet VEV v S , and furthermore eliminates mixing in the Higgs sector. Therefore, the derivation of analytic expressions for the radiative corrections to the matching of the Higgs quartic coupling λ, and their comparison to numerical studies, are significantly simpler, and follow the procedure outlined in Sec. II B. Here we will be interested in the part of the corrections that come on top of the purely SM corrections due to the singlet scalar and shall give expressions including two-loop contributions.
The one-and two-loop corrections to the Higgs mass in the SM are well-known and small; however, in our model there may be large corrections from the singlet scalar. In order to extract the two-loop contributions via Eq. (8) we require the two-loop mass correction, and also the derivatives of the one-loop part. However, our two-loop calculation is performed in the gaugless limit in Feynman gauge, so we require the full one-loop Higgs mass correction in this limit: Here we have defined m 2 S ≡ M 2 S + 1 2 λ SH v 2 , m 2 h ≡ λv 2 which are the tree-level squared masses of the singlet and Higgs respectively, while a complete list of the definitions for our loop functions can be found in Refs. [27,44] which are based upon the basis defined in Refs. [42,45]. This gives us where we use the shorthand notation B(x , y) ≡ ∂ ∂x B(x, y). We note that the infra-red-divergent piece B(0, 0) will cancel against an equivalent piece from ∆ Z2SSM M 2 h , similarly to the effect noted in Ref. [46]. For the Z 2 SSM , we obtain This expression is valid for the gaugeless limit but with generic external momentum (so we can take the momentum in the loop integrals on-shell as the procedure demands, if we wish). However, if we take the "generalised effective potential limit" introduced in Refs. [27,44] and employed in SARAH, then the penultimate line vanishes and the loop functions simplify considerably. We can then obtain a further simplified version of this expression by replacing m 2 h by its tree-level value λv 2 and by performing an expansion in powers of v 2 /m 2 S and keeping only the leading and sub-leading terms, giving

Numerical study
Because the Z 2 symmetry forbids some couplings, the corrections to the matching conditions can be understood in terms of only three parameters added to the SM ones: λ SH , M S , and (to a lesser extent) λ S . The effects of using loop-corrected matching and RGEs in the Z 2 SSM are similar to those observed in Sec. IV A for the SSM, although for most values of λ SH and M S the shift to the quartic coupling has only a very small effect on the value of the cut-off scale. We give in Tab. II our results for λ obtained for the three different orders of matching, for both small and large λ SH and for two choices of M S . For small λ SH the one-loop shift to λ is small, because of a cancellation between the purely-SM part -dominated by the effect of the top quark -and the singlet part of δ (1) λ.  If one then considers larger values of λ SH , the term from the singlet becomes dominant over the SM one, and δ (1) λ is a large negative shift -the evolution of λ, extracted at different orders, as a function of λ SH is also shown in Fig. 5, discussed below. At two loops however, there is no cancellation between SM and singlet contributions, and δ (2) λ is always a negative shift to the Higgs quartic, as was observed previously for the general SSM. On the other hand, it is always small, showing -importantlythat perturbativity of the model is preserved. Having fewer parameters allows for a more detailed study of the different phases of the theory. Indeed, there are two transitions that occur respectively: • Between a metastable and a stable vacuum of the theory; for the physically relevant values of λ around 0.25-0.26, this happens for λ SH ∼ 0.3 and depends very little on M S or λ S .
• Between a UV-complete model and a UVincomplete one -in other words the cut-off scale of the model becomes smaller than the Planck scale for sufficiently large couplings. Fig. 4 shows an example in which the order of the matching performed to extract λ causes differences in the stability of the vacuum of the theory. Indeed, while the curve with two-loop matching and two-loop RGE running (solid line) crosses to negative values of λ -for 10 10 GeV Q 10 16 GeV -the curve with tree-level matching (dashed line) does not, because of the negative shift to the initial value of λ at scale Q = m t at two-loop order. Two-loop corrections to the matching of λ may exclude some parameter points that appear viable when only using a tree-level matching and are therefore important in the discussion of allowed regions of parameter space. Comparing the dashed and dotted lines, we also observe the stabilising effect of the use of the two-loop RGEs, as discussion in Sec. IV A.  are lost below the Planck scale. The left panel of Fig. 5 presents the whole range of couplings that we considered, while the right panel shows an enlargement of the region in which the transition between stable and metastable phases occurs.
We observe that the UV-complete phase of the model corresponds to smaller values of the inputs at scale Q = m t -which can easily be understood as large values of the couplings at m t naturally lead to even larger values at higher scales. Furthermore, we can see that the phase of the model with stable vacua is associated with larger values of λ SH , and that when λ decreases, the value of λ SH needed to ensure a stable vacuum increases. While the SM part of the β-function of λ is negative and tends to drive it to negative values, the additional piece in β λ in the Z 2 SSM is positive and is of the form β When lowering λ(m t ) a higher value of λ SH is needed so that the β-function of λ changes sign earlier, and that λ does not run negative at some scale.
The blue lines in Fig. 5 give λ(m t ) obtained from requiring that m h = 125.1 GeV as a function of λ SH . The different curves correspond to the different orders at which the matching can be done: dotted for tree-level matching, dashed for one-loop and solid for two-loop order. The most important point to notice is that, as for the vacuum stability (see Fig. 4), there is a value of λ SH -here around 0.65 -for which the UV-completeness -in other words whether perturbativity or unitarity is broken at some scale below M P l -of a given parameter point depends greatly on the order at which λ(m t ) has been extracted from the Higgs mass. Moreover, this is not only a matter of using a loop-corrected matching instead of a tree-level one, but the loop order at which it is performed does also matter.
C. Vector-like quarks and stability of the SM From the SM, it is known that the quartic coupling λ runs negative at a scale Q 10 9 − 10 11 GeV, leading to a metastable but long-lived vacuum [1,2]. While extensions with a heavy singlet similar to the previous subsections can have a stabilising effect on the potential [6,7], fermionic extensions typically have the opposite effect through the negative impact of the vector-like (VL) fermions Yukawa coupling on the running of λ, see e.g. [47,48]. A model where the latter is compensated by the effect of the former is discussed in Ref. [22].
Here, we shall extend the SM by one generation of a VL quark doublet Q as well as an up-type quark singlet t with their corresponding counterpartsQ ,t , with quantum numbers under the SM gauge group of t : (3, 1, − 2 3 ), t : (3, 1, 2 3 ), Q : (3, 2, 1 6 ),Q : (3, 2, − 1 6 ). The Lagrangian of the model reads (in terms of two-component spinors) For simplicity we take m Q = m T ≡ M Q ; we then find at one loop that, with the normalisation of the Higgs quartic Let us first consider the impact of the new vector-like states on the running quartic Higgs coupling. For simplicity, we consider here and in the following examples only one extra non-zero Yukawa interaction Y t and consequently setỸ t = 0 as it does not play a role in the following discussion. 4 Then for matching at µ = m t with M Q < TeV, the shifts to λ are less than 10% for Y t 0.7, but grow rapidly to ∼ 50% for Y t ∼ 1. On the other hand, the direct impact of Y t on the running of λ at one loop is given by which contributes significantly to the negative slope of λ for large values but plays a negligible role when Y t is small. In the latter case, the impact of the new fermions on the running of the gauge couplings may outweigh their direct impact on λ. Consider the potential of Eq. (30). Due to the additional coloured fermions, the running of g 3 changes at one-loop to i.e. it decreases more slowly when increasing the scale compared to the SM. In addition, we also obtain a shift in α MS S (m t ) of with respect to the SM. In total, both effects increase the influence of the strong force on the running of λ, adding positively to the slope. The impact on λ is shown in Fig. 6 where the running λ is computed using two-loop RGEs when assuming the pure SM (blue) and the VL extension (purple and black). No matching was applied yet here (i.e. also the shift in α S was neglected) -the changes in the VL case therefore entirely stem from the altered running of the gauge couplings, most importantly Eq. (33). As starting value for λ we used the best-fit value from Ref. [1]. The increased g 3 throughout the energy scales leads to a positive contribution to the slope of λ. As is seen, it can even lead to a stabilisation of the potential at high energies as long as the direct impact of Y t is kept under control by taking it small. This is seen in the purple curve where we have chosen Y t = 0.3. For larger values, the known destabilising effect can overcome the stabilisation from g 3 . As a consequence, the scale of metastability would coincide with the SM for values of Y t ∼ 0.5 and decreases quickly with larger values. This is also shown in the figure for Y t = 0.7 (black line) where λ enters the metastable region already at energies of ∼ 10 5 GeV. We remark that the inclusion of the shift in α S according to Eq. (34) would lead to an even milder running of λ. In fact, just using the one-loop RGEs for the case Y t = 0.3 the quartic coupling stays positive over the entire energy range. We will discuss the effects of the proper matching, including the shifts in λ, in what follows. Solving Eq. (31) for the matched λ VLQ at µ = m t and keepingỸ t = 0, we see that the shifts are slightly negative for small Y t 0.45 and positive for larger Yukawa couplings. This has as a consequence that for low Y t where the VL quarks help increasing the scale of metastability, loop corrections have the opposite effect. However, the size of the shifts stays below 1 %. In Fig. 7, the contour lines represent the predictions for the scale of metastability as a function of Y t and M Q = m Q = m T using tree-level (dotted), one-loop (dashed) as well as two-loop matching (full lines) while applying two-loop RGEs. The colour code in the background quantifies the relative twoloop shift in λ(m t ), (λ (2) − λ (1) )/λ (1) , which stays below roughly half a percent. Nevertheless, the impact of these small shifts is non-negligible: the corresponding Y t values at which λ crosses zero at a given scale typically change by more than 10 % between tree-and two-loop matching. That is, the scale at which metastability occurs is very sensitive to the starting value of λ -meaning that matching is absolutely crucial to make reliable statements. After including the correct shifts in λ, the picture nevertheless remains the same as that for small The situation is reversed if Y t is large. In that case, we enter the known scenario in which the additional impact of Y t on the RGEs of λ drives it negative faster when compared to the SM, further destabilising the vacuum. We show this in Fig. 8. The background shading indicates the scale Λ 0 at which λ crosses zero as a function of Y t and M Q , using two-loop matching, whereas the contour lines represent the relative changes with respect to using one-loop matching, Λ . As expected, Λ 0 is well below the pure SM prediction, and becomes smaller for larger Y t and smaller M Q . The differences in Λ 0 between one-and two-loop matching are quite mild here -"only" O(100 %) -since the two-loop corrections to λ are small. We remark however, that going one order lower and comparing tree-level with one-loop matching using one-loop RGEs, we would see up to an order of magnitude differences in the eventual scale Λ 0 .
Summarising, we have shown that vector-like quarks can have both a destabilising but also a stabilising effect on the SM Higgs potential, and that the correct inclusion of threshold effects is crucial for obtaining precise predictions about the fate of the electroweak vacuum.

D. Two-Higgs Doublet model
Finally, as a last example we study the impact of looplevel matching on Two-Higgs Doublet models. 5 Here we 5 For an overview, see for instance Ref. [49].
will restrict ourselves to the CP-conserving version with a softly broken Z 2 symmetry. The corresponding scalar potential can be written as Note that our sign convention for M 2 12 differs from most definitions in the literature. After electroweak symmetry breaking, we decompose the scalar fields according to where v 2 1 + v 2 2 = v 2 and we define t β = tan β = v 2 /v 1 . The charged (neutral CP-odd) fields mix to one physical charged Higgs H ± (pseudoscalar A) and the corresponding would-be Goldstone bosons. At LO, the angle β coincides with the mixing angle in the pseudoscalar and charged Higgs sector. In the CP-even sector, there are two fields which mix to one light and one heavy eigenstate, with masses m h and m H .
In the same fashion as for the models used above, we can relate the masses and mixing angles to the quartic couplings, leading to the following relations: Analogously to the singlet extension of the SM (Sec IV A), we define the cut-off scale of a particular scenario as the scale at which either one of the λ i becomes larger than 4 π or the unitarity constraints using the running couplings are violated. The latter are too long to show here but can e.g. easily be computed using the SARAH implementation of the model in conjunction with Appendix D of Ref. [28].
First we are going to look at the matching at the top mass scale. It has already been pointed out in Ref. [27] that the loop corrections to the mass spectrum of THDMs can be significant. In Fig. 9, we show on the left-hand side the size of the individual couplings λ i for the three matching orders as a function of the charged Higgs mass.
i |. Note that for λi where i = 3, 4, 5 the NLO differences have been multiplied by an additional factor of 10 to increase visibility. mass according to Eqs. (37) to (41), whereas the λ i evaluated with higher-order matching contain the shifts due to self-energy and tadpole corrections. We see that large differences of O(100 %) or even larger can appear between leading and next-to-leading order. The size of the relative shifts is displayed on the right-hand side of each panel. As expected for a converging perturbative series, the differences between one-and two-loop matching are much less pronounced; however they can still range around tens of percent. Obviously these large differences necessarily have a significant effect on the validity of the theory at higher scales. In the following we will therefore investigate the changes in cut-off scales between the different approaches. As mentioned in the Introduction, the two-loop RGEs are well known but often neglected in the literaturealthough it is known that large differences can appear; see e.g. Ref. [17]. Similar to the singlet-extended SM, the two-loop RGEs tend to moderate the one-loop running. As a result, Landau poles typically appear at much higher scales when including the two-loop effects. For instance, for both i = 1 and 2, 16π 2 β λi ⊃ 24 λ 2 i (1− 13 16π 2 λ i ). The two-loop contribution thereby counteracts the large one-loop slope, stalling the evolution for λ i just below 4 π. In contrast, the two-loop RGEs to λ 3 for instance do have a mitigating effect on the evolution; however a complete stalling only occurs for values much larger than 4π. In Fig. 10 we show for a particular parameter point the running of the couplings λ 1 (black) and λ 3 (blue) for n-loop matching and m-loop RGEs. The influence of the two-loop RGEs can be best gauged for the cases (n, m) = (1, 1) versus (1, 2), displayed as dot-dashed and dashed lines, respectively. It is clearly seen that, while both quartic couplings run into a Landau pole close to Q = 3 TeV when using one-loop RGEs, the inclusion of the two-loop terms leads to a significant flattening and therefore splitting between the two cases.
Let us look at the impact of the threshold corrections next. As shown in the example of Fig. 9, the threshold corrections for the λ i can be significant. This can also be seen in the starting values of the couplings at m t FIG. 11. Comparison of the cut-off scales between using tree matching and one-loop RGEs and both two-loop matching and RGEs for a THDM region with large M 2 12 = −(750 GeV) 2 . The loop-level spectrum was evaluated taking the tree-level values m H ± = 1.14 TeV and tα = −0.95 as inputs. The ratio of the loop-corrected charged Higgs mass to its tree-level input is shown as grey contour lines. The quartic couplings for the case of tree matching were obtained using the leading-order relations Eqs. (37) to (41), taking the spectrum of the two-loop calculation as input. We further fixed tan β = 1.14 and applied the Yukawa scheme of type I.
in Fig. 10: the values for λ 1 using tree-level (one-loop) [two-loop] matching are 1.88 (1.45) [1.14] whereas for λ 3 , the values are 5.7 (4.5) [4.3]. The decrease in value at higher loop orders comes from the fact that in this particular scenario, the average loop corrections to the scalar masses are positive. As a result, one obtains a negative shift in the λ i at the matching scale when demanding an on-shell renormalisation scheme. Consequently, the cut-off scale increases with every additional loop order. Finally, the purple lines show the maximal eigenvalue of the scalar 2 → 2 scattering matrix as well as the corresponding upper bound of 8 π from perturbative unitarity. We see that, while the parameter point would seem to violate perturbative unitarity already at scales below 400 GeV when using the conventional approach, i.e. treelevel matching with one-loop RGEs, it actually only does so just below Q = 2 TeV. Existing studies would have discarded such a parameter point, due to the breakdown of unitarity at energies probeable by the LHC. To that end including matching can result in an increase in the experimental viability of large regions of parameter space.
However, it need not be the case that the cut-off scale is raised by higher loop effects. Indeed, for large values of |M 12 | and therefore large heavy scalar masses, the mass corrections can be large and negative -leading to the opposite effects, i.e. a decreased cut-off scale due to larger quartic couplings after the inclusion of the proper matching. An example where this happens is presented in Fig. 11. For this figure we have evaluated the spectrum at the two-loop level while fixing the tree-level input values of t α as well as m H ± which enter the spectrum calculation at the loop level. Therefore, the loop-corrected m H ± varies over this plane. The grey contours show the ratio of the loop-corrected charged Higgs mass over the treelevel input, m (1) H ± /m (T ) H ± . To obtain the LO couplings, i.e. the case of tree-level matching, we take the loop-corrected spectrum and calculate λ i according to Eqs. (37) to (41). Finally, we run the couplings up in scale using two-loop RGEs for the two-loop-and one-loop RGEs for the treelevel-matched couplings in order to evaluate the cut-off scale. The coloured contours show the ratio of the cutoff scales, Λ (T,1) /Λ (2,2) , obtained with tree-level matching and one-loop RGEs and with two-loop matching and two-loop RGEs, respectively. In particular in the region where all heavy scalar masses are approximately equal, we observe large differences in cut-off scales. In fact, while the tree-level matching approach suggests a cut-off at O(10 7 GeV), the full two-loop matching procedure demands new physics restoring unitarity and perturbativity already at the TeV scale.
Concluding, the conventional approach of tree-level matching and one-loop RGEs can both over-but also underestimate the cut-off scale by many orders of magnitude. It is therefore of crucial importance to (i) take into account the -known -RGEs beyond one-loop and to (ii) consistently match the couplings before running.

V. SUMMARY AND CONCLUSIONS
In this paper, we have investigated the impact that matching plays in the high-scale validity of minimal extensions of the Standard Model (SM). We argued that the usual approach of using N − 1 matching when utilising N -loop RGEs neglects important contributions in the presence of large couplings. In fact, for most nonsupersymmetric models, studies beyond tree-level matching and one-loop RGEs are rare or even absent. We analysed in different scenarios the impact of both matching at two-loop order as well as the two-loop RGEs, highlighting the differences with respect to previous approaches. For simple models, we provided an analytical computation of the matching conditions. We pointed out how sensitive the cut-off scale of the real-singlet-extended SM is to the loop order of both matching and RGE running and showed that the scale dependence decreases for Nwith respect to N − 1-loop matching. Imposing an additional Z 2 symmetry on this model furthermore enabled us to study the fate of the electroweak vacuum as well as the UV completion analytically. We highlighted regions of parameter space where the model can in principle be valid up to the Planck scale -a statement which crucially depends on the proper matching of the quartic couplings at the low scale.
In a scenario where the SM is extended by vector-like quarks, we showed that the impact of the latter can ac-tually increase the Higgs quartic interaction such that it does not become negative at higher scales -an observation which we have not encountered before in the literature. The reason is that, despite the negative impact of the additional Yukawa coupling on the running of λ, the presence of additional coloured states modifies the running strong coupling in such a way that it adds positively to the β-function of λ. Also in this scenario, the matching of λ before the RGE evolution has a significant impact on the predicted high-scale behaviour of the model.
As a final example we showed in a Two-Higgs Doublet model that the loop-level matching of the quartic couplings can lead to significant changes in both the MS values and subsequently the cut-off scale of the theory.
To conclude, we observe that robust statements about the UV behaviour of non-supersymmetric, weakly coupled BSM models can only be made when including, at the very least, loop-level matching. We stress that the required loop-level corrections, as well as two-loop RGEs, are readily accessible with the computer tool SARAH for any general renormalisable field theory. In light of our results, we strongly encourage its use when accurate high-scale predictions are required.