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

Johannes Braathen, Mark D. Goodsell, Manuel E Krauss, Toby Opferkuch,

Florian Staub

► To cite this version:

Johannes Braathen, Mark D. Goodsell, Manuel E Krauss, Toby Opferkuch, Florian Staub. N-loop running should be combined with N-loop matching. Physical Review D, 2018, 97 (1), pp.015011. 10.1103/PhysRevD.97.015011. hal-01701025

HAL Id: hal-01701025

https://hal.sorbonne-universite.fr/hal-01701025

Submitted on 5 Feb 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
N-loop running should be combined with N-loop matching

Johannes Braathen,¹,²,* Mark D. Goodsell,¹,²,† Manuel E. Krauss,³,‡ Toby Opferkuch,³,⁴,§ and Florian Staub⁵,⁶,||

¹LPTHE, UPMC Univ. Paris 6, Sorbonne Universités, 4 Place Jussieu, F-75252 Paris, France

²LPTHE, CNRS, 4 Place Jussieu, F-75252 Paris, France

³Bethe Center for Theoretical Physics & Physikalisches Institut der Universität Bonn,
Nußallee 12, 53115 Bonn, Germany

⁴PRISMA Cluster of Excellence & Mainz Institute for Theoretical Physics,
Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany

⁵Institute for Theoretical Physics (ITP), Karlsruhe Institute of Technology,
Engesserstraße 7, D-76128 Karlsruhe, Germany

⁶Institute for Nuclear Physics (IKP), Karlsruhe Institute of Technology,
Hermann-von-Helmholtz-Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany

<figure><img src="image_2.png" /></figure>

(Received 1 December 2017; published 22 January 2018)

We investigate the high-scale behavior of Higgs sectors beyond the Standard Model, pointing out that the proper matching of the quartic couplings before applying the renormalization 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 behavior, 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 sizable 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 stabilization of the electroweak vacuum at high energies, in contrast to results in the literature.

DOI: 10.1103/PhysRevD.97.015011

# 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 environment—i.e. the models typically contain the minimal numbers of new fields to exhibit novel phenomenology. All additional states also alter the high-scale behavior 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–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 probed by the LHC, or if the lifetime of the electroweak breaking vacuum is too short on cosmological time scales.

Many of these effects have already been studied in the literature for a plethora of different models such as singlet extensions [5–9], triplet extensions [10,11], two-Higgs-doublet models (THDMs) [12–20] or models with vector-like fermions [21]. These studies utilize the one- and sometimes even two-loop renormalization group equations (RGEs). However, less care was taken in the determination of the parameters which enter the RGE running. Often, two-loop RGEs were combined with a tree-level matching.

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.

<footnote>*braathen@lpthe.jussieu.fr</footnote>

<footnote>†goodsell@lpthe.jussieu.fr</footnote>

<footnote>‡mkrauss@th.physik.uni-bonn.de</footnote>

<footnote>§toopferk@uni-mainz.de</footnote>

<footnote>||florian.staub@kit.edu</footnote>
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].<fnref n="1" /> It was shown that even the next-to-next-to-leading-order (NNLO) shifts to $\lambda$ 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_i$: the corresponding shifts in $\lambda$ 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. [26] in the example of THDMs where the two-loop corrections to the Higgs masses were calculated for the first time.<fnref n="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 behavior 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 stabilize 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 organized at follows. In Sec. II we give a step-by-step prescription for the general matching procedure including effects at the N-1- and 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 formulas.

## 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.

### A. Renormalization group equations

We shall always work in the $\overline{\text{MS}}$ scheme. In this scheme, the $\beta$ functions, which describe the energy dependence of the parameters $\Theta$, are defined as

$$ \beta_i = \mu \frac{d\Theta_i}{d\mu}. \qquad (1) $$

Here, $\mu$ is an arbitrary mass scale. $\beta_i$ can be expanded in a perturbative series:

$$ \beta_i = \sum_n \frac{1}{(16\pi^2)^n} \beta_i^{(n)}. \qquad (2) $$

$\beta_i^{(1)}$ and $\beta_i^{(2)}$ 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. [28–31].

### B. Matching

The renormalized coupling constants $\Theta_i$ in $d = 4 - 2\epsilon$ dimensions, which enter the running, are related to the corresponding bare couplings $\Theta_i^0$ by

$$ \Theta_i^0 \mu^{-C_i \epsilon} = \Theta_i + \sum_n \frac{a_i^{(n)}}{\epsilon^n}. \qquad (3) $$

Here, $C_i$ are constant factors depending on the character of $\Theta_i$.<fnref n="3" /> The coefficients $a_i$ are the result of a perturbative expansion. In general, two approaches are possible to determine the $\overline{\text{MS}}$ parameters as functions 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 ($\delta\Theta_i^{\text{OS}}$).

<footnote n="3">³Gauge and Yukawa couplings have $C_i = 1$, quartic couplings $C_i = 2$.</footnote>

<footnote n="1">¹Loop corrections in the scalar sector were taken into account in Ref. [22] for a singlet extension and in Refs. [23,24] for a THDM. These studies did not however investigate the impact on the high-scale behavior of the model. In Ref. [25] 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.</footnote>

<footnote n="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. [27].</footnote>
(2) In an $\overline{\text{MS}}$ calculation the counterterms of the Lagrangian parameters $(\delta\Theta_i^{\overline{\text{MS}}})$ 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

$$ \Theta_i^0 = \Theta_i^{\text{OS}} - \delta\Theta_i^{\text{OS}} = \Theta_i(\mu) - \delta\Theta_i^{\overline{\text{MS}}}. \quad (4) $$

In an on-shell calculation, however, there is no generic set of renormalization group equations known, and therefore to explore a theory at high energies it is necessary to use $\overline{\text{MS}}$ equations—i.e. to extract the underlying $\overline{\text{MS}}$ parameters of the theory, and then run them. On the other hand, in an $\overline{\text{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 $\overline{\text{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

$$ M_h^2 = \lambda v^2 + \sum_{n=1}^{\infty} \frac{1}{(16\pi^2)^n} \Delta^{(n)} M_h^2(\lambda). \quad (5) $$

This is in general a highly nonlinear equation in $\lambda$; but fortunately since the series is perturbative we can solve it through expanding

$$ \lambda = \lambda^{(0)} + \frac{1}{16\pi^2} \delta^{(1)}\lambda + \frac{1}{(16\pi^2)^2} \delta^{(2)}\lambda + \dots \quad (6) $$

to find

$$ \lambda^{(0)} = \frac{M_h^2}{v^2}, $$

$$ \delta^{(1)}\lambda = -\frac{1}{v^2} \Delta^{(1)} M_h^2|_{\lambda=\lambda^{(0)}}, $$

$$ \delta^{(2)}\lambda = -\frac{1}{v^2} \left[ \delta^{(1)}\lambda \frac{\partial}{\partial\lambda} \Delta^{(1)} M_h^2 + \Delta^{(2)} M_h^2 \right]_{\lambda=\lambda^{(0)}}, \quad (7) $$

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 $\lambda$, as we shall argue in the next subsection, when we are using N-loop RGEs, to obtain a consistent

expansion of our $\overline{\text{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 + \frac{4}{3}(Q_1^2 + Q_2^2)$, and we take the gauge coupling at high energies to be $g(\Lambda)$. Now, to determine the gauge coupling at a low-energy $\mu$ 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 $\mu$. This gives the one-loop value of the gauge coupling to be

$$ \frac{8\pi^2}{g^2(\mu)} = \frac{8\pi^2}{g^2(\Lambda)} - \left( b_0 + \frac{4}{3}Q_1^2 + \frac{4}{3}Q_2^2 \right) \log \frac{M_2}{\Lambda} - \left( b_0 + \frac{4}{3}Q_1^2 \right) \log \frac{M_1}{M_2} - b_0 \log \frac{\mu}{M_1}. \quad (8) $$

This can be rewritten as

$$ \frac{8\pi^2}{g^2(\mu)} = \frac{8\pi^2}{g^2(\Lambda)} - b_0 \log \frac{\mu}{\Lambda} - \frac{4}{3} Q_1^2 \log \frac{M_1}{\Lambda} - \frac{4}{3} Q_2^2 \log \frac{M_2}{\Lambda}, $$

which is also the result of simply matching the two theories at the scale $\mu$ 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_2$—for 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

$$ \begin{aligned} \frac{8\pi^2}{\tilde{g}^2(\mu)} &= \frac{8\pi^2}{g^2(\Lambda)} - \left( b_0 + \frac{4}{3}Q_1^2 + \frac{4}{3}Q_2^2 \right) \log \frac{M}{\Lambda} - b_0 \log \frac{\mu}{M} \\ &= \frac{8\pi^2}{g^2(\mu)} - \frac{4}{3}Q_1^2 \log \frac{M_1}{M} - \frac{4}{3}Q_2^2 \log \frac{M_2}{M}. \end{aligned} \quad (9) $$

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. Another example, more relevant for this work, is the quartic coupling of a real scalar $\phi$ coupled to a Majorana (Weyl) fermion $\psi$:

$$ \mathcal{L} \supset -\frac{1}{24} \lambda \phi^4 - \left( \frac{1}{2} (m + y\phi)\psi\psi + \text{H.c.} \right). \quad (10) $$

If $y$ and $m$ are real, then the one-loop threshold correction to $\lambda$ from integrating out the fermion $\psi$ at matching scale $M$ is

$$ \delta\lambda = -\frac{y^2}{16\pi^2} \left[ \mathbf{32y^2} + (24y^2 - 4\lambda) \log \frac{m}{M} \right]. \quad (11) $$

The finite, nonlogarithmic correction is highlighted in bold. On the other hand,

$$ \frac{d\lambda}{d\log \mu} = -\frac{1}{16\pi^2} [y^2(24y^2 - 4\lambda) - 3\lambda^2]. \quad (12) $$

Hence the tree-level matching clearly finds the correct logarithmic terms, but we miss the finite correction of $-\frac{2y^2}{\pi^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. $\frac{2y^2}{\pi^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 cutoff, 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 $\mu$ we find a Landau pole at approximately

$$ \Lambda \approx M \exp \left( -\frac{16\pi^2}{3\lambda(M)} \right), \quad (13) $$

but if we have incorrectly determined $\lambda(M)$ by an amount $\delta\lambda$, then the ratio of the correct cutoff scale to the erroneous one $\Lambda'$ is

$$ \begin{aligned} \frac{\Lambda'}{\Lambda} &\approx \exp \left( \frac{16\pi^2}{3} \frac{\delta\lambda}{\lambda^2(M)} \right) \\ &\approx 200 \exp \left( \frac{16\pi^2}{3} \left( \frac{\delta\lambda}{\lambda^2(M)} - 0.1 \right) \right). \end{aligned} \quad (14) $$

In the Standard Model, the difference between the tree-level value for $\lambda$ 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 minimal supersymmetric SM where $\delta\lambda \gtrsim \lambda$).

### III. NUMERICAL SETUP

For the numerical calculations we make use of the MATHEMATICA package SARAH [32–37] to produce a spectrum generator based on SPHENO [38–40]. 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 gaugeless limit and based on the generic results of Refs. [41,42]. In nonsupersymmetric BSM models, these results suffer in general from the so-called “Goldstone boson catastrophe” even in the gaugeless 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. [43] 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 $\overline{\text{MS}}$ parameters:

1. The running couplings $\Theta(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 $\mu_i^2|\phi|^2$.

(3) The tree-level mass are calculated by diagonalizing the tree-level mass matrices

(4) The n-loop corrections to the tadpoles $\delta^{(n)}t_i$ are calculated. The imposed renormalization conditions are

$$T_i + \sum_j^{n} \delta^{(j)} t_i = 0, \qquad (15)$$

which cause shifts in $\mu_i^2$:

$$\mu_i^2 \rightarrow \mu_i^2 + \sum_j^n \delta^{(j)} \mu_i^2. \qquad (16)$$

(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 loop-corrected mass matrix calculated as

$$M_{\phi}^{(n)}(p^2) = \tilde{M}_{\phi}^{(2L)} - \sum_{j}^{n} \Pi_{\phi}^{(j)}(p^2). \qquad (17)$$

Here, $\tilde{M}_{\phi}$ 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

$$[\mathrm{eig} M_{\phi}^{(n)}(p^2 = m_{\phi_i}^2)]_i \equiv m_{\phi_i}^2 \qquad (18)$$

is fulfilled. The renormalized rotation matrix is taken to be the one calculated for $p^2 = m_{\phi_1}^2$.

If a chosen set of input parameters $\Theta(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 leading-order (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 $\overline{\text{MS}}$ couplings corresponding to the desired physical parameters at loop level is nontrivial. In what follows we use a simple fitting algorithm which varies the input parameters until the desired masses and mixing angles are obtained.

# IV. MODELS AND RESULTS

## A. Singlet extension

We start with the SM extended by a real gauge singlet $S$. The potential reads

$$V = \mu^2|H|^2 + \frac{1}{2}M_{\tilde{S}}^2S^2 + \kappa_1|H|^2S + \frac{1}{3}\kappa_2S^3 \\ + \frac{1}{2}\lambda|H|^4 + \frac{1}{2}\lambda_{SH}S^2|H|^2 + \frac{1}{2}\lambda_S S^4. \qquad (19)$$

After electroweak symmetry breaking the CP-even scalar components mix to two physical states $h$, $H$ via a rotation angle $\alpha$. At tree level we can use $m_h$, $m_H$ and $t_\alpha \equiv \tan \alpha$ as input to calculate the quartic couplings

$$\lambda_H = \frac{m_h^2 + m_H^2 t_\alpha^2}{v^2(1 + t_\alpha^2)}, \qquad (20)$$

$$\lambda_S = \frac{\kappa_1 v^2}{8 v_S^3} - \frac{\kappa_2}{4 v_S} + \frac{(m_H^2 + m_h^2 t_\alpha^2)}{4(1 + t_\alpha^2) v_S^2}, \qquad (21)$$

$$\lambda_{SH} = - \frac{m_H^2 t_\alpha - m_h^2 t_\alpha + \kappa_1 v + \kappa_1 t_\alpha^2 v}{v v_S (1 + t_\alpha^2)}, \qquad (22)$$

where $v$ and $v_S$ are the vacuum expectation values (VEVs) of the Higgs doublet and singlet scalar respectively. Interpreting $m_h$, $m_H$ and $t_\alpha$ as physical on-shell parameters, the quartic calculated via Eqs. (20)–(22) are therefore the LO values.

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 \rightarrow 2$ scattering amplitude in the limit of high energies and demand that the eigenvalues stay below $8\pi$. The conditions for the model under consideration read

$$|\lambda_{SH}| < 8\pi, \qquad (23)$$

$$|\lambda| < 8\pi, \qquad (24)$$

$$|-6\lambda_S - 3\lambda \pm \sqrt{4\lambda_{SH}^2 + 9(\lambda - 2\lambda_S)^2}| < 16\pi, \qquad (25)$$

which leads to $\lambda_S < 4\pi/3$ for small $\lambda$ and $\lambda_{SH}$. Although 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 (at both one- and 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 one-loop 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
<figure><img src="image_3.png" /><figcaption>FIG. 1. Values of the running quartic couplings at the scale $Q = 10^6$ GeV using one-loop (left) and two-loop RGEs (right) as functions of the matching scale at which the quartic couplings were calculated. The labels $(n, m)$ refer to $n$-loop-level matching of the quartics and $m$-loop RGEs. We choose the parameters of the singlet extended SM at the matching scale to be $m_h = 125$ GeV, $m_H = 400$ GeV, $\tan\alpha = 0.3$ and $v_S = 300$ GeV. Cubic terms were set to zero to ensure a scale-invariant input. Note that the y axes' ranges are different in each panel.</figcaption></figure>

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 two-loop RGEs. This also means that large differences are expected when evaluating the cutoff scale of a theory, i.e. the scale at which the model becomes nonperturbative or violates unitarity. In Table I, we show the cutoff 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 nonperturbative ($\Lambda_{4\pi}$) separately to the case where either perturbativity or unitarity is violated ($\Lambda_{4\pi}^{\text{unit}}$). 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

<table><caption>TABLE I. Values of the quartic couplings and the cutoff for different combinations of parameters, $\lambda^{(n)}$, and RGEs, $\beta^{(m)}$ at the loop orders $n$ and $m$ respectively. $\Lambda_{4\pi}$ is the scale at which the quartic couplings first exceed $4\pi$ while $\Lambda_{4\pi}^{\text{unit}}$ is the naive $4\pi$ cutoff augmented by the unitarity constraint of Eq. (25). The observables are fixed at the $n$-loop order to be $m_h = 125$ GeV, $m_H = 700$ GeV, $\tan\alpha = 0.1$ while the remaining input parameters are chosen as $\kappa_1 = 0$ GeV, $\kappa_2 = 2000$ GeV, $v_S = 175$ GeV.</caption><thead><tr><th>(n, m)</th><th>λ</th><th>λ<sub>S</sub></th><th>λ<sub>SH</sub></th><th>Λ<sub>4π</sub> [GeV]</th><th>Λ<sub>4π</sub><sup>unit</sup> [GeV]</th></tr></thead><tbody><tr><td>(T, 1)</td><td rowspan="2">0.34</td><td rowspan="2">1.1</td><td rowspan="2">-1.1</td><td>6.4 × 10<sup>3</sup></td><td>3.2 × 10<sup>3</sup></td></tr><tr><td>(T, 2)</td><td>8.0 × 10<sup>6</sup></td><td>1.3 × 10<sup>4</sup></td></tr><tr><td>(1, 1)</td><td rowspan="2">0.33</td><td rowspan="2">0.24</td><td rowspan="2">-0.97</td><td>6.4 × 10<sup>8</sup></td><td>3.2 × 10<sup>8</sup></td></tr><tr><td>(1, 2)</td><td>1.3 × 10<sup>12</sup></td><td>2.5 × 10<sup>9</sup></td></tr><tr><td>(2, 1)</td><td rowspan="2">0.32</td><td rowspan="2">0.18</td><td rowspan="2">-0.94</td><td>5.1 × 10<sup>10</sup></td><td>2.5 × 10<sup>10</sup></td></tr><tr><td>(2, 2)</td><td>1.0 × 10<sup>14</sup></td><td>2.0 × 10<sup>11</sup></td></tr></tbody></table>

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 $\beta$ function of $\lambda_S$ is given by $\beta_{\lambda_S}^{(1)} = \frac{1}{16\pi^2}(36\lambda_S^2 + \lambda_{SH}^2)$, so that $\lambda_S$ tends to grow very rapidly, there is a moderating term from the two-loop RGEs which goes with $\frac{1}{(16\pi^2)^2}(-816\lambda_S^3 - 20\lambda_S\lambda_{SH}^2)$. Therefore, using the one-loop RGEs only, $\lambda_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 $\lambda_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

<figure><img src="image_4.png" /><figcaption>FIG. 2. The running of the quartic couplings for the point given in Table I. The line styles refer to the loop order of the matching and RGE running as described in Table I, namely $(n, m)$ refers to the matching at $n$-loop order with $m$-loop RGEs. The solid red line is the $4\pi$ perturbativity limit, while the dashed red line is the unitarity constraint of $4\pi/3$ obtained from Eq. (25) in the limit $\lambda_S \gg \lambda_{SH}, \lambda$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Difference in the predicted cutoff scale depending on the matching performed as a function of the singlet VEV $v_S$ and the heavy CP-even Higgs mass $m_H$. Left: We show the ratio of the obtained cutoff given matching at $N$ versus $N-1$ order using the RGEs at $N$-loop order. The colored (grey) contours use the two (one)-loop RGEs, therefore showing the ratio of the matching at two (one) loop versus one loop (tree level), respectively. Right: Ratio of the most precise calculation performed (both matching and RGEs at two-loop order) versus the leading order (tree-level matching and one-loop RGEs). The grey contours correspond to the ratios of the quartic coupling $\lambda_S$ for these two scenarios. Here we have fixed the physical parameters such that $m_h = 125$ GeV, $\tan \alpha = 0.2$, while the remaining parameters are chosen as $\kappa_1 = 0$ GeV and $\kappa_2 = 1000$ GeV.</figcaption></figure>

$\lambda_{SH}$ and $\lambda$ is not as pronounced. For $\lambda_{SH}$, the corresponding $\beta$ function grows with $\frac{1}{16\pi^2}$, $12\lambda_S\lambda_{SH}$ with only a small moderating effect from the two-loop RGEs. As a consequence, it becomes larger than $4\pi$ before $\lambda_S$ and then drags the latter with it. In total, in particular because of the large two-loop contributions to $\beta_{\lambda_S}$, there can be several orders of magnitude between the eventual cutoff 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 $\overline{\text{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 cutoff scale by several orders of magnitude when including the matching.

Finally in 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 cutoff scale here and in what follows is defined as the scale at which either one of the couplings grows larger than $4\pi$ or any of the conditions for perturbative unitarity are violated, each evaluated with the running $\overline{\text{MS}}$ quartic couplings. The grey contours in the left-hand panel of Fig. 3 display the ratio of the evaluated cutoff scales for $N=1$. 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 colored contours) are significantly milder in this region: the maximum difference is just a factor of 3. For large $v_S$, instead, the quartic couplings are comparably small, leading to large cutoff 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 \gtrsim 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 cutoff 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 cutoff scale can be many orders of magnitude larger than the cutoff 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, $\lambda_S^{(2)}/\lambda_S^{(T)}$. Already at the matching scale, differences of an order of magnitude between tree and two-loop matching can appear, emphasizing the requirement for proper matching and running when analyzing the high-scale behavior of a given model.

## B. Singlet extension with an additional $\mathbb{Z}_2$ symmetry

### 1. Analytical approximation

We may now make a further simplification to the singlet extension studied in Sec. IV A, namely adding an additional $\mathbb{Z}_2$ symmetry under which the singlet scalar is charged—for clarity we will call this model the $\mathbb{Z}_2$ singlet-extended Standard Model (SSM) to distinguish it from the SSM. This symmetry forbids nonzero values for the couplings $\kappa_1, \kappa_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 $\lambda$, 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 gaugeless limit in Feynman gauge, so we require the full one-loop Higgs mass correction in this limit:

$$ \Delta^{(1)} M_h^2 (p^2) = 3y_t^2 (4m_t^2 - p^2) B(m_t^2, m_t^2) - \frac{3}{2} \lambda^2 v^2 B(0, 0) \\ - \frac{9}{2} \lambda^2 v^2 B(m_h^2, m_h^2) - \frac{1}{2} v^2 \lambda_{SH}^2 B(m_S^2, m_S^2). \qquad (26) $$

Here we have defined $m_S^2 = M_S^2 + \frac{1}{2} \lambda_{SH} v^2$, $m_h^2 = \lambda 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. [26,43] which are based upon the basis defined in Refs. [41,44]. This gives us

$$ \begin{aligned} \lambda^{(0)} &= \frac{M_h^2}{v^2} \\ \delta^{(1)} \lambda &= \delta_{SM}^{(1)} \lambda + \frac{1}{2} \lambda_{SH}^2 B(m_S^2, m_S^2) \\ \delta^{(2)} \lambda &= \delta_{SM}^{(2)} \lambda - \frac{1}{v^2} [\Delta_{Z_2SSM}^{(2)} M_h^2]_{\lambda=\lambda^{(0)}} \\ & \quad + \frac{3}{2} \lambda_{SH}^2 \lambda B(m_S^2, m_S^2) (3B(m_h^2, m_h^2) \\ & \quad + 3m_h^2 B((m_h^2)', m_h^2) + B(0, 0)) \end{aligned} \qquad (27) $$

where we use the shorthand notation $B(x', y) \equiv \frac{\partial}{\partial x} B(x, y)$. We note that the infrared-divergent piece $B(0, 0)$ will cancel against an equivalent piece from $\Delta_{Z_2SSM}^{(2)} M_h^2$, similarly to the effect noted in Ref. [45].

For the $Z_2$SSM, we obtain

$$ \begin{aligned} \delta^{(2)} \lambda = \delta_{SM}^{(2)} \lambda & - \frac{1}{2} \lambda_{SH} (\lambda_{SH}^3 v^2 M_{SSSSS} (m_S^2, m_S^2, m_S^2, m_S^2, m_h^2) + 6 \lambda \lambda_{SH}^2 v^2 M_{SSSSS} (m_h^2, m_S^2, m_h^2, m_S^2, m_S^2) \\ & - 6 \lambda \lambda_{SH} U_{SSSS} (m_h^2, m_h^2, m_S^2, m_S^2) - 4 \lambda_{SH}^2 U_{SSSS} (m_S^2, m_S^2, m_h^2, m_S^2) + 9 \lambda^2 \lambda_{SH} v^2 V_{SSSSS} (m_h^2, m_h^2, m_S^2, m_S^2) \\ & + 2 \lambda_{SH}^3 v^2 V_{SSSSS} (m_S^2, m_S^2, m_h^2, m_S^2) - \lambda_{SH}^2 Y_{SSSS} (m_S^2, m_S^2, m_h^2, m_S^2) - 9 \lambda^2 Y_{SSSS} (m_h^2, m_h^2, m_S^2, m_S^2) \\ & - 12 \lambda_{SH} \lambda_S Y_{SSSS} (m_S^2, m_S^2, m_h^2, m_S^2) - 6 \lambda_{SH} \lambda_S Z_{SSSS} (m_S^2, m_S^2, m_h^2, m_S^2) - 12 \lambda \lambda_{SH} Z_{SSSS} (m_h^2, m_h^2, m_S^2, m_S^2) \\ & - \frac{1}{2v^2} \lambda_{SH}^2 (S_{SSS} (m_h^2, m_S^2, m_S^2) - I(m_h^2, m_S^2, m_S^2)) + \frac{9}{2} \lambda_{SH}^2 \lambda m_h^2 B(m_S^2, m_S^2) B((m_h^2)', m_h^2). \end{aligned} \qquad (28) $$

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 “generalized effective potential limit” introduced in Refs. [26,43] 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_h^2$ by its tree-level value $\lambda v^2$ and by performing an expansion in powers of $v^2/m_S^2$ and keeping only the leading and subleading terms, giving

$$ \begin{aligned} \delta^{(2)} \lambda \approx \delta_{SM}^{(2)} \lambda & - \frac{9}{4v^2} \lambda_{SH} \lambda A(m_S^2) + \lambda_{SH}^3 [1 - 2\overline{\log} m_S^2 + \overline{\log}^2 m_S^2] \\ & + \frac{1}{4} \lambda_{SH}^2 \lambda [-18 - 6\overline{\log}^2 m_S^2 + (36\overline{\log} m_h^2 - 12)\overline{\log} m_S^2] \\ & + 3\lambda_{SH}^2 \lambda_S [-1 + \overline{\log} m_S^2 + \overline{\log}^2 m_S^2]. \end{aligned} \qquad (29) $$

## 2. 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: $\lambda_{SH}$, $M_S$, and (to a lesser extent) $\lambda_S$. The effects of using loop-corrected matching and RGEs in the $Z_2$SSM are similar to those observed in Sec. IVA for the SSM, although for most values of $\lambda_{SH}$ and $M_S$ the shift to the quartic coupling has only a very small effect on the value of the cutoff scale. We give in Table II our results for $\lambda$ obtained for the three different orders of matching, for both small and large $\lambda_{SH}$ and for two choices of $M_S$. For small $\lambda_{SH}$ the one-loop shift to $\lambda$ is small, because of a cancellation between the purely-SM part—dominated by the effect of the top quark—and the singlet part of $\delta^{(1)}\lambda$. If one then considers larger values of $\lambda_{SH}$, the term from the singlet becomes dominant over the SM one, and $\delta^{(1)}\lambda$ is a large negative shift—the evolution of $\lambda$, extracted at different orders, as a function of $\lambda_{SH}$ is also shown in Fig. 5, discussed below. At two loops however, there is no cancellation between SM and singlet contributions, and $\delta^{(2)}\lambda$ 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—importantly—that perturbativity of the model is preserved.
<table><caption>TABLE II. Values of the Higgs quartic $\lambda$ obtained from matchings at tree-level, one-loop and two-loop orders, for different choices of $M_S$ and $\lambda_{SH}$. The singlet quartic coupling $\lambda_S$ is set to be 0.1.</caption><thead><tr><th>$M_S$ [GeV]</th><th>$\lambda_{SH}$</th><th>$\lambda_{tree}$</th><th>$\lambda_{1\ell}$</th><th>$\lambda_{2\ell}$</th></tr></thead><tbody><tr><td rowspan="2">500</td><td>0.25</td><td>0.2610</td><td>0.2623</td><td>0.2551</td></tr><tr><td>3</td><td>0.2610</td><td>0.1885</td><td>0.1794</td></tr><tr><td rowspan="2">1000</td><td>0.25</td><td>0.2610</td><td>0.2651</td><td>0.2546</td></tr><tr><td>3</td><td>0.2610</td><td>0.1548</td><td>0.1385</td></tr></tbody></table>

Having fewer parameters allows for a more detailed study of the different phases of the theory. Indeed, there are two transitions that occur respectively:

(1) Between a metastable and a stable vacuum of the theory; for the physically relevant values of $\lambda$ around 0.25–0.26, this happens for $\lambda_{SH} \sim 0.3$ and depends very little on $M_S$ or $\lambda_S$.

(2) Between a UV-complete model and a UV-incomplete one—in other words the cutoff scale of the model becomes smaller than the Planck scale for sufficiently large couplings.

Figure 4 shows an example in which the order of the matching performed to extract $\lambda$ 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 $\lambda$—for $10^{10}$ GeV $\lesssim Q \lesssim 10^{16}$ GeV—the curve with tree-level

<figure><img src="image_6.png" /><figcaption>FIG. 4. Running of the Higgs quartic coupling as a function of the renormalization scale $Q$, having taken $\lambda_{SH} = 0.28$, $\lambda_S = 0.1$ and $M_S = 500$ GeV. The value of $\lambda(m_t)$ is obtained by requiring $m_h = 125.15$ GeV, with different orders of matchings depending on the curve. The solid line corresponds to the use of two-loop matching and two-loop RGE running, the dashed line to tree-level matching and two-loop RGEs, and the dotted line to tree-level matching and one-loop RGEs. Note that because of the cancellation that occurs in the one-loop correction for small $\lambda_{SH}$ (discussed in the main text), the curves we would have obtained using one-loop matching would have been very similar to those with tree-level matching.</figcaption></figure>

matching (dashed line) does not, because of the negative shift to the initial value of $\lambda$ at scale $Q = m_t$ at two-loop order. Two-loop corrections to the matching of $\lambda$ 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 stabilizing effect of the use of the two-loop RGEs, as discussion in Sec. IVA.

Figure 5 shows how both types of transitions occur in this model. The different domains in this figure were obtained as follows: we start with values of the couplings, at the scale $Q = m_t$, in the range $\lambda \in [0, 0.35]$ and $\lambda_{SH} \in [0, 4]$, and take $\lambda_S = 0.1$ and $M_S = 500$ GeV. We then use two-loop RGEs to run the couplings up to the Planck scale, and we verify whether the Higgs quartic $\lambda$ becomes negative at any point, and whether perturbativity or unitarity 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 the 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 $\lambda_{SH}$, and that when $\lambda$ decreases, the value of $\lambda_{SH}$ needed to ensure a stable vacuum increases. While the SM part of the $\beta$ function of $\lambda$ is negative and tends to drive it to negative values, the additional piece in $\beta_\lambda$ in the $Z_2$SSM is positive and is of the form $\beta_\lambda^{(1)} \supset \frac{1}{16\pi^2}\lambda_{SH}^2$. When lowering $\lambda(m_t)$ a higher value of $\lambda_{SH}$ is needed so that the $\beta$ function of $\lambda$ changes sign earlier, and that $\lambda$ does not run negative at some scale.

The blue lines in Fig. 5 give $\lambda(m_t)$ obtained from requiring that $m_h = 125.1$ GeV as a function of $\lambda_{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 $\lambda_{SH}$—here around 0.65—for which the UV completeness—in other words whether perturbativity or unitarity is broken at some scale below $M_{Pl}$—of a given parameter point depends greatly on the order at which $\lambda(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 $\lambda$ runs negative at a scale $Q \approx 10^9-10^{11}$ GeV, leading to a
<figure><img src="image_7.png" /><figcaption>FIG. 5. Different phases in the $Z_2$SSM shown in the $\lambda_{SH}-\lambda$ plane, where the couplings are taken at the scale $Q = m_t$. The orange shaded region of parameter space corresponds to UV-complete theories, i.e. none of the three quartic couplings ($\lambda, \lambda_{SH}, \lambda_S$) become nonperturbative and the constraints from unitarity are not violated before the Planck scale; the black shaded region corresponds to theories with stable vacua. The thin blue lines give $\lambda$ as a function of $\lambda_{SH}$ when imposing $m_h = 125.15$ GeV with a matching condition at respectively tree level (dotted curve), one-loop order (dashed curve) and two-loop order (solid line). The other parameters of the scalar sector are $\lambda_S = 0.1$, and $M_S = 500$ GeV.</figcaption></figure>

metastable but long-lived vacuum [1,2]. While extensions with a heavy singlet similar to the previous subsections can have a stabilizing 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 $\lambda$; see e.g. Refs. [46,47]. A model where the latter is compensated by the effect of the former was discussed in Ref. [21].

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 counterparts $\tilde{Q}', \tilde{t}',$ with quantum numbers under the SM gauge group of $t': (\tilde{3}, 1, -\frac{2}{3})$, $\tilde{t}': (\tilde{3}, 1, \frac{2}{3})$, $Q': (\tilde{3}, 2, \frac{1}{6})$, $\tilde{Q}': (\tilde{3}, 2, -\frac{1}{6})$. The Lagrangian of the model reads (in terms of two-component spinors)

$$ \mathcal{L} = \mathcal{L}_{\text{SM}} - (Y'_t Q' \cdot H t' + \tilde{Y}'_t \tilde{Q}' \cdot \tilde{H} t' + m_T \tilde{t}' t' + m_Q \tilde{Q}' Q' + \text{H.c.}). \quad (30) $$

For simplicity we take $m_Q = m_T \equiv M_Q$; we then find at one loop that, with the normalization of the Higgs quartic coupling $\mathcal{L} \supset -\frac{1}{2}\lambda|H|^4$

$$ \lambda_{\text{SM}} = \lambda_{\text{VLQ}} - \frac{1}{16\pi^2} \left[ (Y'_t + \tilde{Y}'_t)^2 (5Y'^2_t - 2Y'_t \tilde{Y}'_t + 5\tilde{Y}'^2_t) \right. \\ \left. \qquad + 6(Y'^4_t + \tilde{Y}'^4_t) \log \frac{M_Q^2}{\mu^2} \right] \\ + \frac{2\lambda_{\text{VLQ}}}{16\pi^2} \left[ (Y'_t - \tilde{Y}'_t)^2 + 3(Y'^2_t + \tilde{Y}'^2_t) \log \frac{M_Q^2}{\mu^2} \right]. \quad (31) $$

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 nonzero Yukawa interaction $Y'_t$ and consequently set $\tilde{Y}'_t = 0$ as it does not play a role in the following discussion.<fnref n="4" /> Then for matching at $\mu = m_t$ with $M_Q < \text{TeV}$, the shifts to $\lambda$ are less than 10% for $Y'_t \lesssim 0.7$, but grow rapidly to ~50% for $Y'_t \sim 1$. On the other hand, the direct impact of $Y'_t$ on the running of $\lambda$ at one loop is given by

$$ 16\pi^2 \beta_{\lambda}^{(1)} \supset 12Y_{t}^{2}(\lambda - Y_{t}^{2}), \quad (32) $$

which contributes significantly to the negative slope of $\lambda$ 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 $\lambda$. Consider the potential of Eq. (30). Due to the additional colored fermions, the running of $g_3$ changes at one loop to

$$ 16\pi^2 \beta_{g_3}^{(1)} = \left( -7 + \frac{4}{3}n_T + \frac{2}{3}n_Q \right) g_3^3 \rightarrow -5g_3^3, \quad (33) $$

i.e. it decreases more slowly when increasing the scale compared to the SM. In addition, we also obtain a shift in $\alpha_S^{\overline{\text{MS}}}(m_t)$ of

$$ \alpha_S^{\overline{\text{MS}}} \rightarrow \frac{\alpha_S^{\overline{\text{MS}}}}{1 - \frac{\alpha_S^{\overline{\text{MS}}}}{\pi} \log(M_Q/m_t)} \quad (34) $$

with respect to the SM. In total, both effects increase the influence of the strong force on the running of $\lambda$, adding

<footnote n="4">⁴Although this leads to a stable lightest VL quark, there could for instance be couplings to a hidden sector, leading to a relaxation of the direct collider constraints.</footnote>
<figure><img src="image_8.png" /><figcaption>FIG. 6. Simplified comparison between the running of $\lambda$ in the SM with and without vector-like states. Here, we used full one-loop (dashed lines) and two-loop RGEs (full lines) in both models and as a starting point the SM best-fit values from Ref. [1]. For the purple (black) lines we use $Y'_t = 0.3(0.7)$.</figcaption></figure>

positively to the slope. The impact on $\lambda$ is shown in Fig. 6 where the running $\lambda$ 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 $\alpha_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 a starting value for $\lambda$ 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 $\lambda$. As is seen, it can even lead to a stabilization 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 destabilizing effect can overcome the stabilization from $g_3$. As a consequence, the scale of metastability would coincide with the SM for values of $Y'_t \sim 0.5$ and decreases quickly with larger values. This is also shown in the figure for $Y'_t = 0.7$ (black line) where $\lambda$ enters the metastable region already at energies of $\sim 10^5$ GeV. We remark that the inclusion of the shift in $\alpha_S$ according to Eq. (34) would lead to an even milder running of $\lambda$. 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 $\lambda$, in what follows.

Solving Eq. (31) for the matched $\lambda_{VLQ}$ at $\mu = m_t$ and keeping $\tilde{Y}'_t = 0$, we see that the shifts are slightly negative for small $Y'_t \lesssim 0.45$ and positive for larger Yukawa couplings. This has as a consequence that for low $Y'_t$ where the VL quarks help to increase 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

<figure><img src="image_9.png" /><figcaption>FIG. 7. Contours of the scale $\Lambda_0$ at which $\lambda$ runs negative with full two-loop RGEs and one-loop (dashed) or two-loop (full lines) matching, and for comparison when using the best-fit value $\lambda^{SM} = 0.25208$ at the top mass scale (dotted). Black lines correspond to $\Lambda_0 = 10^{18}$ GeV, and blue lines to $10^{11}$ GeV. The background shows the two-loop shift in $\lambda(m_t)$ in percent, defined as $(\lambda^{(2)} - \lambda^{(1)})/\lambda^{(1)}$.</figcaption></figure>

lines) while applying two-loop RGEs. The color code in the background quantifies the relative two-loop shift in $\lambda(m_t)$, $(\lambda^{(2)} - \lambda^{(1)})/\lambda^{(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 $\lambda$ 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 $\lambda$—meaning that matching is absolutely crucial to make reliable statements. After including the correct shifts in $\lambda$, the picture nevertheless remains the same as that for small $Y'_t$; the impact of the VL quarks on $\alpha_S$ can be such that the scale of metastability is increased with respect to

<figure><img src="image_10.png" /><figcaption>FIG. 8. The scale of metastability $\Lambda_0$ in the case of large $Y'_t$ using two-loop running with two-loop matching for $\lambda$. The black contours show the size of $\Lambda_0$ with respect to using one-loop matching.</figcaption></figure>
the SM, leading to the possibility of absolute stability all the way up to the Planck scale.

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 $\lambda$ drives it negative faster when compared to the SM, further destabilizing the vacuum. We show this in Fig. 8. The background shading indicates the scale $\Lambda_0$ at which $\lambda$ 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, $\Lambda_0^{(2,2)}/\Lambda_0^{(1,2)}$. As expected, $\Lambda_0$ is well below the pure SM prediction, and becomes smaller for larger $Y'_t$ and smaller $M_Q$. The differences in $\Lambda_0$ between one- and two-loop matching are quite mild here—"only" $\mathcal{O}(100\%)$— since the two-loop corrections to $\lambda$ 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 $\Lambda_0$.

Summarizing, we have shown that vector-like quarks can have both a destabilizing but also a stabilizing 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 loop-level matching on two-Higgs-doublet models.<fnref n="5" /> Here we will restrict ourselves to the CP-conserving version with a softly broken $Z_2$ symmetry. The corresponding scalar potential can be written as

$$V = m_1^2 \Phi_1^\dagger \Phi_1 + m_2^2 \Phi_2^\dagger \Phi_2 + \lambda_1 (\Phi_1^\dagger \Phi_1)^2 + \lambda_2 (\Phi_2^\dagger \Phi_2)^2 \\ + \lambda_3 (\Phi_1^\dagger \Phi_1)(\Phi_2^\dagger \Phi_2) + \lambda_4 (\Phi_2^\dagger \Phi_1)(\Phi_1^\dagger \Phi_2) \\ + M_{12}^2 (\Phi_1^\dagger \Phi_2 + \Phi_2^\dagger \Phi_1) + \frac{\lambda_5}{2} ((\Phi_2^\dagger \Phi_1)^2 + (\Phi_1^\dagger \Phi_2)^2). \quad (35)$$

Note that our sign convention for $M_{12}^2$ differs from most definitions in the literature. After electroweak symmetry breaking, we decompose the scalar fields according to

$$\Phi_k = \begin{pmatrix} \phi_k^+ \\ \frac{1}{\sqrt{2}}(v_k + \phi_k^0 + i\sigma_k) \end{pmatrix}, \quad i = 1, 2, \quad (36)$$

where $v_1^2 + v_2^2 = v^2$ and we define $t_\beta = \tan\beta = v_2/v_1$. The charged (neutral CP-odd) fields mix to one physical charged Higgs $H^\pm$ (pseudoscalar A) and the corresponding would-be Goldstone bosons. At LO, the angle $\beta$ 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:

$$\lambda_1 = \frac{1 + t_{\beta}^2}{2(1 + t_{\alpha}^2)v^2} (m_h^2 t_{\alpha}^2 + m_H^2 + M_{12}^2 t_{\beta} (1 + t_{\alpha}^2)), \quad (37)$$

$$\lambda_2 = \frac{M_{12}^2 (1 + t_{\beta}^2)}{2 t_{\beta}^3 v^2} + \frac{(1 + t_{\beta}^2) (m_h^2 + m_H^2 t_{\alpha}^2)}{2 t_{\beta}^2 (1 + t_{\alpha}^2) v^2}, \quad (38)$$

$$\lambda_3 = \frac{1}{(1 + t_{\alpha}^2) t_{\beta} v^2} [(m_H^2 - m_h^2) t_{\alpha} (1 + t_{\beta}^2) \\ + 2 m_{H^{\pm}}^2 (1 + t_{\alpha}^2) t_{\beta} + M_{12}^2 (1 + t_{\alpha}^2) (1 + t_{\beta}^2)], \quad (39)$$

$$\lambda_4 = \frac{1}{t_{\beta} v^2} (-M_{12}^2 (1 + t_{\beta}^2) + m_A^2 t_{\beta} - 2 m_{H^{\pm}}^2 t_{\beta}), \quad (40)$$

$$\lambda_5 = \frac{1}{t_{\beta} v^2} (-M_{12}^2 (1 + t_{\beta}^2) - m_A^2 t_{\beta}). \quad (41)$$

Analogously to the singlet extension of the SM (Sec. IVA), we define the cutoff scale of a particular scenario as the scale at which either one of the $\lambda_i$ becomes larger than $4\pi$ 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. [27].

First we are going to look at the matching at the top mass scale. It has already been pointed out in Ref. [26] 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 $\lambda_i$ for the three matching orders as a function of the charged Higgs mass. The leading-order $\lambda_i$ are simple linear functions of this mass according to Eqs. (37)–(41), whereas the $\lambda_i$ evaluated with higher-order matching contain the shifts due to self-energy and tadpole corrections. We see that large differences of $\mathcal{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 cutoff scales between the different approaches.

As mentioned in the Introduction, the two-loop RGEs are well known but often neglected in the literature—although 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

<footnote n="5">⁵For an overview, see for instance Ref. [48].</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 9. The size of the quartic couplings at $m_t$ at LO (blue dotted), NLO (grey dashed) and NNLO (blue solid) as a function of $m_{H^±}$. Here, we set the other physical masses and mixing angle to $m_h=125$ GeV, $m_H=750$ GeV, $m_A=730$ GeV and $\tan\alpha = -0.71$. The other model parameters are $\tan\beta = 1.4$ and $M_{12}^2 = -500^2$ GeV$^2$. On the right, we show the relative differences $\delta\lambda_i^{\text{LO}} = |(\lambda_i^{(T)} - \lambda_i^{(1)})/\lambda_i^{(T)}|$ and $\delta\lambda_i^{\text{NLO}} = |(\lambda_i^{(1)} - \lambda_i^{(2)})/\lambda_i^{(1)}|$. Note that for $\lambda_i$ where $i = 3, 4, 5$ the NLO differences have been multiplied by an additional factor of 10 to increase visibility.</figcaption></figure>

scales when including the two-loop effects. For instance, for both $i = 1$ and 2, $16\pi^2\beta_{\lambda_i} \supset 24\lambda_i^2(1 - \frac{13}{16\pi^2}\lambda_i)$. The two-loop contribution thereby counteracts the large one-loop slope, stalling the evolution for $\lambda_i$ just below $4\pi$. In contrast, the two-loop RGEs to $\lambda_3$ for instance do have a mitigating effect on the evolution; however a complete stalling only occurs for values much larger than $4\pi$. In Fig. 10 we show for a particular parameter point the running of the couplings $\lambda_1$ (black) and $\lambda_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

<figure><img src="image_12.png" /><figcaption>FIG. 10. RGE running of the individual quartic couplings $\lambda_1$ (black) and $\lambda_3$ (blue) for the parameter point defined by $m_H = 511$ GeV, $m_A = 607$ GeV, $m_{H^±} = 605$ GeV, $t_\beta = 1.45$, $t_\alpha = -0.82$, $M_{12}^2 = -(250 \text{ GeV})^2$, using n-loop-level matching and m-loop RGEs. The dotted lines stand for (T, 1), dot-dashed for (1,1), dashed for (1,2) and full lines for (2,2). Coupling values of $\pm 4\pi$ are indicated by a red solid line. In addition, we display the largest eigenvalue of the scalar $2 \to 2$ scattering amplitude in purple. The upper bound of $8\pi$ is indicated by a dashed red line. We used the Yukawa scheme of type II.</figcaption></figure>

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 $\lambda_i$ can be significant. This can also be seen in the starting values of the couplings at $m_t$ in Fig. 10: the values for $\lambda_1$ using tree-level (one-loop) [two-loop] matching are 1.88 (1.45) [1.14] whereas for $\lambda_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 $\lambda_i$ at the matching scale when demanding an on-shell renormalization scheme. Consequently, the cutoff scale increases with every additional loop order. Finally, the purple lines show the maximal eigenvalue of the scalar $2 \to 2$ scattering matrix as well as the corresponding upper bound of $8\pi$ 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. tree-level 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 cutoff 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
<figure><img src="image_13.png" /><figcaption>FIG. 11. Comparison of the cutoff scales between using tree matching and one-loop RGEs and both two-loop matching and RGEs for a THDM region with large $M_{12}^2 = -(750 \text{ 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 (37)–(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.</figcaption></figure>

opposite effects, i.e. a decreased cutoff 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 tree-level input, $m_{H^±}^{(1)}/m_{H^±}^{(T)}$. To obtain the LO couplings, i.e. the case of tree-level matching, we take the loop-corrected spectrum and calculate $\lambda_i$ according to Eqs. (37)–(41). Finally, we run the couplings up in scale using two-loop RGEs for the two-loop and one-loop RGEs for the tree-level-matched couplings in order to evaluate the cutoff scale. The colored contours show the ratio of the cutoff scales, $\Lambda^{(T,1)}/\Lambda^{(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 cutoff scales. In fact, while the tree-level matching approach suggests a cutoff at $O(10^7 \text{ 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 cutoff 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 SM. We argued that the usual approach of using $N-1$ matching when utilizing $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 analyzed 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 cutoff 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 $N$- with 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 actually 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 $\lambda$, the presence of additional colored states modifies the running strong coupling in such a way that it adds positively to the $\beta$ function of $\lambda$. Also in this scenario, the matching of $\lambda$ before the RGE evolution has a significant impact on the predicted high-scale behavior 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 $\overline{\text{MS}}$ values and subsequently the cutoff scale of the theory.

To conclude, we observe that robust statements about the UV behavior of nonsupersymmetric, 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 renormalizable field theory. In light of our results, we strongly encourage its use when accurate high-scale predictions are required.

ACKNOWLEDGMENTS

M. E. K., T. O. and F. S. thank the LPTHE in Paris for their hospitality while part of this work was completed. J. B. and M. D. G. acknowledge support from French state funds managed by the Agence Nationale de la Recherche
(ANR), in the context of the LABEX ILP (ANR-11-IDEX-0004-02, ANR-10-LABX-63), and M. D. G. acknowledges support from the ANR grant "HiggsAutomator" (ANR-15-CE31-0002). J. B. was supported by a scholarship from the Fondation CFM. M. E. K. is supported by the DFG Research Unit 2239 "New Physics at the LHC." The work of T. O. has been supported by the German Research

Foundation (DFG) under grant Nos. EXC-1098, FOR 2239 and SFB/Transregio 33 "The Dark Universe" as well as by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 637506, "*Directions*"). F. S. is supported by the ERC Recognition Award ERC-RA-0008 of the Helmholtz Association.

[1] D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio, and A. Strumia, J. High Energy Phys. 12 (2013) 089.

[2] G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori, and A. Strumia, J. High Energy Phys. 08 (2012) 098.

[3] A. Andreassen, W. Frost, and M. D. Schwartz, arXiv: 1707.08124.

[4] A. Spencer-Smith, arXiv:1405.1975.

[5] R. N. Lerner and J. McDonald, Phys. Rev. D 80, 123507 (2009).

[6] O. Lebedev, Eur. Phys. J. C 72, 2058 (2012).

[7] J. Elias-Miro, J. R. Espinosa, G. F. Giudice, H. M. Lee, and A. Strumia, J. High Energy Phys. 06 (2012) 031.

[8] G. M. Pruna and T. Robens, Phys. Rev. D 88, 115012 (2013).

[9] R. Costa, A. P. Morais, M. O. P. Sampaio, and R. Santos, Phys. Rev. D 92, 025024 (2015).

[10] Y. Hamada, K. Kawana, and K. Tsumura, Phys. Lett. B 747, 238 (2015).

[11] N. Khan, arXiv:1610.03178.

[12] H. S. Cheon and S. K. Kang, J. High Energy Phys. 09 (2013) 085.

[13] N. Chakrabarty, U. K. Dey, and B. Mukhopadhyaya, J. High Energy Phys. 12 (2014) 166.

[14] N. Chakrabarty and B. Mukhopadhyaya, Eur. Phys. J. C 77, 153 (2017).

[15] P. Ferreira, H. E. Haber, and E. Santos, Phys. Rev. D 92, 033003 (2015); 94, 059903(E) (2016).

[16] N. Chakrabarty and B. Mukhopadhyaya, Phys. Rev. D 96, 035028 (2017).

[17] D. Chowdhury and O. Eberhardt, J. High Energy Phys. 11 (2015) 052.

[18] S. Gori, H. E. Haber, and E. Santos, J. High Energy Phys. 06 (2017) 110.

[19] P. Basler, P. M. Ferreira, M. Muhlleitner, and R. Santos, arXiv:1710.10410.

[20] P. S. Bhupal Dev and A. Pilaftsis, J. High Energy Phys. 12 (2014) 024; 11 (2015) 147(E).

[21] M.-L. Xiao and J.-H. Yu, Phys. Rev. D 90, 014007 (2014); 90, 019901(E) (2014).

[22] S. Kanemura, M. Kikuchi, and K. Yagyu, Nucl. Phys. B907, 286 (2016).

[23] P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt, and A. Wlotzka, J. High Energy Phys. 02 (2017) 121.

[24] P. Basler, M. Muhlleitner, and J. Wittbrodt, arXiv: 1711.04097.

[25] A. Kobakhidze and A. Spencer-Smith, J. High Energy Phys. 08 (2013) 036.

[26] J. Braathen, M. D. Goodsell, and F. Staub, Eur. Phys. J. C 77, 757 (2017).

[27] M. E. Krauss and F. Staub, arXiv:1709.03501.

[28] M. E. Machacek and M. T. Vaughn, Nucl. Phys. B222, 83 (1983).

[29] M. E. Machacek and M. T. Vaughn, Nucl. Phys. B236, 221 (1984).

[30] M. E. Machacek and M. T. Vaughn, Nucl. Phys. B249, 70 (1985).

[31] M.-x. Luo, H.-w. Wang, and Y. Xiao, Phys. Rev. D 67, 065019 (2003).

[32] F. Staub, arXiv:0806.0538.

[33] F. Staub, Comput. Phys. Commun. 181, 1077 (2010).

[34] F. Staub, Comput. Phys. Commun. 182, 808 (2011).

[35] F. Staub, Comput. Phys. Commun. 184, 1792 (2013).

[36] F. Staub, Comput. Phys. Commun. 185, 1773 (2014).

[37] F. Staub, Adv. High Energy Phys. 2015, 840780 (2015).

[38] W. Porod, Comput. Phys. Commun. 153, 275 (2003).

[39] W. Porod and F. Staub, Comput. Phys. Commun. 183, 2458 (2012).

[40] F. Staub and W. Porod, Eur. Phys. J. C 77, 338 (2017).

[41] S. P. Martin, Phys. Rev. D 70, 016005 (2004).

[42] M. Goodsell, K. Nickel, and F. Staub, Eur. Phys. J. C 75, 290 (2015).

[43] J. Braathen and M. D. Goodsell, J. High Energy Phys. 12 (2016) 056.

[44] S. P. Martin, Phys. Rev. D 68, 075002 (2003).

[45] S. P. Martin and D. G. Robertson, Phys. Rev. D 90, 073010 (2014).

[46] J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori, A. Riotto, and A. Strumia, Phys. Lett. B 709, 222 (2012).

[47] G. Bambhaniya, P. Bhupal Dev, S. Goswami, S. Khan, and W. Rodejohann, Phys. Rev. D 95, 095016 (2017).

[48] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, Phys. Rep. 516, 1 (2012).