Leading two-loop corrections to the Higgs boson self-couplings in models with extended scalar sectors

Abstract

We compute the dominant two-loop corrections to the Higgs trilinear coupling $\lambda _{hhh}$ and to the Higgs quartic coupling $\lambda _{hhhh}$ in models with extended Higgs sectors, using the effective-potential approximation. We provide in this paper all necessary details about our calculations, and present general $\overline{{\mathrm{MS}}}$ expressions for derivatives of the integrals appearing in the effective potential at two loops. We also consider three particular Beyond-the-Standard-Model (BSM) scenarios – namely a typical scenario of an Inert Doublet Model (IDM), and scenarios of a Two-Higgs-Doublet Model (2HDM) and of a Higgs Singlet Model (HSM) without scalar mixing – and we include all the necessary finite counterterms to obtain (in addition to $\overline{{\mathrm{MS}}}$ results) on-shell scheme expressions for the corrections to the Higgs self-couplings. With these analytic results, we investigate the possible magnitude of two-loop BSM contributions to the Higgs self-couplings and the fate of the non-decoupling effects that are known to appear at one loop. We find that, at least as long as pertubative unitarity conditions are fulfilled, the size of two-loop corrections remains well below that of one-loop corrections. Typically, two-loop contributions to $\lambda _{hhh}$ amount to approximately 20% of those at one loop, implying that the non-decoupling effects observed at one loop are not significantly modified, but also meaning that higher-order corrections need to be taken into account for the future perspective of precise measurements of the Higgs trilinear coupling.

1 Introduction

The discovery of a 125-GeV Higgs boson at the CERN LHC [1, 2] in 2012 was an astonishing success for particle physics, establishing the mechanism of Electroweak Symmetry Breaking (EWSB) and completing the particle content of the Standard Model (SM). Nevertheless, there is no doubt that new physics is needed to address deficiencies of the SM, both because of theoretical considerations and of a number of experimental results. At the same time, there has so far been no clear evidence of what this new physics would be, and instead the many ongoing experiments are only setting increasingly stringent constraints on the parameter space of possible Beyond-the-Standard-Model (BSM) theories, leaving us without any guidance about how to address the shortcomings of the SM. One of the most pressing and timely questions in this respect is to understand the structure of the Higgs sector. Indeed, many – if not most – of the best motivated extensions of the SM come with enlarged Higgs sectors – for example supersymmetric models, models with additional gauge symmetries, or bottom-up models to realise e.g. baryogenesis, scalar dark matter, etc. – and the Higgs sector is thus expected to play a special, central role in BSM searches. However, to this point, all measured Higgs-boson properties are in agreement with SM predictions within experimental and theoretical uncertainties (see for example Ref. [3]). This seems to suggest that the BSM scalar states are somehow made difficult to find, being either heavy and beyond the reach of current experiments (and possibly even decoupled if very heavy), or hidden by some symmetry or mechanism.

One example of the latter is alignment [4], which occurs in extended scalar sectors with multiple Higgs doublets when one of the CP-even Higgs mass eigenstates is collinear in field space with the total electroweak vacuum expectation value (VEV). The state aligned with the VEV then obtains SM-like couplings at tree level, while the other scalars are difficult to detect – in particular their couplings to weak gauge bosons vanish in this limit. The alignment can in principle happen in two distinct cases: (i) as a consequence of the decoupling of the additional scalars, or (ii) without decoupling. This second case can happen, for example, because of some symmetry – e.g. in the Maximally Symmetric Two-Higgs-Doublet Model (2HDM) [5] or in inert scalar models with an unbroken $\mathbb {Z}_2$ symmetry [6,7,8] – or of possible dynamics in the ultra-violet (UV); for an example of how an aligned 2HDM can appear as a low-energy limit of a supersymmetric theory with Dirac gauginos see for instance Refs. [9, 10].

However, even in aligned BSM scenarios, properties of the 125-GeV Higgs boson can deviate from their SM predictions in various observables, because of radiative corrections involving the new BSM states. For the Higgs-boson couplings, as found first in Refs. [11, 12], these loop corrections can actually become very significant in some regions of the parameter space of the BSM theories, because of non-decoupling effects. Among the Higgs properties, those that can exhibit the largest such non-decoupling effects are its self-couplings, i.e. its trilinear coupling $\lambda _{hhh}$ and its quartic coupling $\lambda _{hhhh}$.¹

Indeed, both couplings are directly related to the shape of the Higgs potential, of which currently only very little is known with the exception of the existence of the electroweak (EW) minimum and the curvature of the potential around this minimum (determined by the Higgs mass of 125 GeV). While the EW minimum and Higgs mass are common for all BSM models, the Higgs self-couplings provide information about the differences between models beyond the SM. Additionally, the Higgs self-couplings determine the strength of the electroweak phase transition (EWPT). In particular, successful scenarios of electroweak baryogenesis (EWBG) [14,15,16] require the EWPT to be of strong first order, and it has been shown in Refs. [17, 18] that for this to be the case, $\lambda _{hhh}$ must deviate from its SM value by at least $20-30\%$. $\lambda _{hhhh}$ is also important when considering the behaviour of the Higgs potential at large field values because it relates to the Lagrangian scalar quartic couplings, whose running to high scales controls the stability of the Higgs potential – this has been studied at loop level in the SM in Refs. [19,20,21], and in BSM extensions in e.g. Refs. [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35].

Contrary to most of the couplings of the 125-GeV Higgs boson that are now known to a precision of at least a few percent, the Higgs self-couplings are currently not well constrained experimentally and deviations of several hundred percent from their SM values are still allowed at the LHC. For $\lambda _{hhh}$, some limits are already available: using single Higgs production data from LHC Run 2, the ATLAS collaboration set a limit on the ratio

$$\begin{aligned} \kappa _\lambda \equiv \frac{\lambda ^\text {exp.}_{hhh}}{\lambda _{hhh}^\text {SM}} \end{aligned}$$

(1.1)

as $-3.2<\kappa _\lambda <11.9$ at 95% confidence level (CL) [36], while with double Higgs production, the best intervals obtained at 95% CL are respectively $-5.0<\kappa _\lambda <12.1$ from ATLAS [37] (see also Ref. [38]) and $-11<\kappa _\lambda <17$ from CMS [39] (see also Ref. [40]). These measurements are expected to be significantly improved at future colliders. State-of-the-art values for the expected accuracies for the ratio $\kappa _\lambda$ at almost all envisioned future colliders can be found in Ref. [41]. We only recall here some of the main results, considering moreover sensitivities obtained through exclusive analyses (that are typically weakened when considering more complete global fits). First of all, the high-luminosity upgrade of the LHC (the HL-LHC) could reach $0.5<\kappa _\lambda <1.6$ (at 68% CL) with an integrated luminosity of 3$\text { ab}^{-1}$ [42] (see also Ref. [43]). A possible high-energy version of the LHC (the HE-LHC), with centre-of-mass energy of $\sqrt{s}=$ 27 TeV, was found in Ref. [44] to be able to attain $0.54<\lambda _{hhh}/\lambda _{hhh}^\text {SM}<1.46$ (at 68% CL) using 15$\text { ab}^{-1}$ of data.² Turning next to the case of possible future lepton colliders, the initial stage of the International Linear Collider (ILC) running at $\sqrt{s}=$ 250 GeV cannot access $\lambda _{hhh}$ directly via double-Higgs production [46], but could obtain a measurement to 49% accuracy, at 68% CL, in a single-Higgs production analysis using 2$\text { ab}^{-1}$ of data [41]. With the data from further ILC extensions to 500 GeV (4$\text { ab}^{-1}$) or even 1 TeV (8$\text { ab}^{-1}$), it could be possible to reach a precision of 27% or 10% respectively [47] (once again at 68% CL). The CLIC project could, using the combination of 1$\text { ab}^{-1}$ of data at 380 GeV, 2.5$\text { ab}^{-1}$ at 1.5 TeV, and 5$\text { ab}^{-1}$ at 3 TeV, obtain a final accuracy of $0.93<\kappa _\lambda <1.11$ at 68% confidence level [48] (see also Refs. [49, 50]). Further in the future, a 100-TeV FCC-hh hadron collider with 30$\text { ab}^{-1}$ of data could allow reaching a $5\%$ accuracy (at 68% CL) on the measurement of the Higgs trilinear coupling [42, 45] (see also Ref. [51]). Finally, it will most probably not be possible to probe the Higgs quartic coupling experimentally in a foreseeable future, the involved cross sections being far too small – for instance, with 30$\text { ab}^{-1}$ of data, the FCC-hh would only be able to constrain the ratio $\lambda _{hhhh}^\text {exp.}/\lambda _{hhhh}^\text {SM}$ to be approximately between $-4$ and 16 [52] (at 95% confidence level).

On the theoretical side, the first one-loop calculations of $\lambda _{hhh}$ were performed in the SM and the minimal supersymmetric SM (MSSM) in Refs. [53,54,55]. One-loop corrections to the Higgs trilinear coupling have also been studied in several non-supersymmetric extensions of the SM, with singlets [56,57,58,59], additional doublets [11, 12, 59,60,61,62], or triplets [63]. Expressions for $\lambda _{hhh}$, as well as loop calculations of all other Higgs couplings and of Higgs decays, are implemented in the public program H-COUP [64, 65] for the (CP-conserving) 2HDM, the Inert Doublet Model (IDM) and the singlet extension of the SM. Various other public tools also include, or allow obtaining, results for loop corrections to Higgs couplings and decays in these three BSM models: namely 2HDMC [66], PROPHECY4F [67], and 2HDECAY [68] for the 2HDM; sHDECAY [69] and PROPHECY4F (since version 3.0 [70]) for the singlet extension; and finally SPheno [71, 72] together with SARAH [73,74,75,76,77,78] which implement expressions for generic BSM theories [79] that can be applied automatically to a desired model. Since the early works [11, 12] in the context of the 2HDM, it has been known that the radiative corrections involving additional BSM scalars can, in the non-decoupling regime, cause a significant enhancement of $\lambda _{hhh}$ – with deviations from its SM prediction of several tens of percent or even a (few) hundred percent. One should emphasise that there is in principle no problem in finding one-loop corrections larger than the tree-level result, as the loop corrections here do not stem from a perturbation of the tree-level formula, but instead arise from new parameters that enter the calculation only at the one-loop level. Furthermore, the large effects found in Refs. [11, 12] were obtained for parameter points satisfying the criterion of tree-level perturbative unitarity [80] (expressed for the 2HDM in Refs. [81, 82]). Nevertheless, one is quite naturally led to ask what would happen once corrections beyond the one-loop level are included, and in particular whether new huge effects can appear at two loops or not.

Two-loop corrections to $\lambda _{hhh}$ have so far only been considered in a limited number of works in the literature, with different motivations. The earliest calculations were performed in the context of supersymmetric models, namely in the MSSM [83] and Next-to-MSSM (NMSSM) [84], in which Higgs boson masses can be and are calculated to high precision, making it necessary to compute also $\lambda _{hhh}$ to a similar level of accuracy.³ In Refs. [83, 84], the leading $\mathscr {O}(\alpha _s\alpha _t)$ corrections to $\lambda _{hhh}$ at two loops, computed using the effective-potential approximation, were found to be approximately $5-10\%$ of the size of the one-loop corrections, and allowed a significant reduction of the scale dependence of the total results. In addition to this, a third reference, Ref. [85], studied (part of) the leading corrections from the additional scalars in the Inert Doublet Model and how these affect the strength of the EWPT. This calculation found an enhancement of $\lambda _{hhh}$ by a few percent at two loops even if effects of 30–40% appear at one loop; in turn these two-loop contributions slightly weaken the strength of the first-order EWPT.

In Ref. [86], we also computed two-loop corrections to $\lambda _{hhh}$ in an aligned 2HDM and in the IDM, but in that respect, we took a slightly different point of view compared to previous works. Indeed, what we wanted to investigate was the maximal possible size of the two-loop corrections and how non-decoupling effects can be affected by them. We found that two-loop corrections amount typically to 10–20% of the one-loop corrections, and hence while they do not alter dramatically the non-decoupling effect that might appear, they are not entirely negligible.

In the present paper, we therefore build on our previous work, and we provide all needed details about our calculations and the involved technical aspects. We moreover extend both our computations, by including also the case of a singlet extension of the SM – which we will refer to as Higgs Singlet Model (HSM) – and our numerical investigations. In addition to $\lambda _{hhh}$, we also give new formulae for the two-loop corrections to $\lambda _{hhhh}$ in these models. Once again, our main interest is to determine the maximal possible size of the BSM deviations, so we consider scenarios without mixing throughout this work.

This paper is organised as follows: we start by defining our notations for the models that we study in Sect. 2, before describing the set-up of our calculations in Sect. 3. In Sect. 4, we give general results for derivatives of the effective potential, expressed in the $\overline{{\mathrm{MS}}}$ scheme. Then we present our analytical results for the SM and three BSM scenarios, in both $\overline{{\mathrm{MS}}}$ and on-shell schemes, in Sect. 5 and consider numerical examples in Sect. 6. Finally, we discuss implications of our calculations in Sect. 7, before concluding in Sect. 8. Additional details are presented in appendices, with our conventions and definitions of loop functions in Appendix A, full expressions for 2HDMs in Appendix B, and definitions of the intermediate functions used in Sect. 4 in Appendix C.

2 Models

We here recall our conventions for the 2HDM, the IDM, and the HSM, and describe the scenarios that we will be considering. For more complete reviews of these models, see for example Refs. [87,88,89,90] for the 2HDM, Refs. [6, 8, 91] for the IDM, and Ref. [92] for the HSM.

2.1 Two-Higgs-Doublet-Models

We consider first a CP-conserving Two-Higgs-Doublet Model [93], in which the Higgs sector is composed of two $SU(2)_L$–doublets of hypercharge $Y=1/2$. This type of model is in principle plagued by possible large Higgs-mediated flavour-changing neutral currents (FCNCs) at tree level, which would be incompatible with experimental results. Such tree-level FCNCs can however be avoided by requiring that each type of fermion only couples to one of the two Higgs doublets [94, 95], and this can be achieved by imposing a $\mathbb {Z}_2$ symmetry under which the scalar doublets transform respectively as $\varPhi _1\rightarrow \varPhi _1$ and $\varPhi _2\rightarrow -\varPhi _2$, and the different families of fermions have charges $\pm 1$. Several charge assignments are possible for the fermion families, corresponding to distinct types of 2HDM [96,97,98] – but as we only consider effects from top quarks in the following we do not need to specify a type here.⁴ The $\mathbb {Z}_2$ symmetry can be broken softly – i.e. without reintroducing dangerous FCNCs – via an off-diagonal mass term $m_3^2\big (\varPhi _2^\dagger \varPhi _1+\text {h.c.}\big )$. We follow the conventions of Ref. [12] and write the tree-level scalar potential as

$$\begin{aligned} V^{(0)}_\text {2HDM}=&\ m_{1}^2 |\varPhi _1|^2 + m_{2}^2 |\varPhi _2|^2 - m_3^2 \big (\varPhi _2^\dagger \varPhi _1 + \text {h.c.}\big )\nonumber \\&+\frac{\lambda _1}{2}|\varPhi _1|^4 +\frac{\lambda _2}{2}|\varPhi _2|^4 + \lambda _3 |\varPhi _1|^2|\varPhi _2|^2\nonumber \\&+\lambda _4 |\varPhi _2^\dagger \varPhi _1|^2 +\frac{\lambda _5}{2} \Big ((\varPhi _2^\dagger \varPhi _1)^2 + \text {h.c.}\Big )\,. \end{aligned}$$

(2.1)

Our assumption that CP is conserved in the Higgs sector has also allowed us to take all parameters in the above equation – as well as the VEVs of both doublets – to be real. Requiring that this potential is bounded from below implies the following conditions [6, 22, 23, 100, 101]

$$\begin{aligned} \lambda _1>0\,,\quad \lambda _2>0\,,\quad \sqrt{\lambda _1\lambda _2} +\lambda _3+\text {min}\{0,\lambda _4\pm \lambda _5\}>0\,. \end{aligned}$$

(2.2)

We expand each of the scalar doublets as [12]

$$\begin{aligned} \varPhi _i=\begin{pmatrix} w_i^+\\ \frac{1}{\sqrt{2}}(v_i+h_i+i z_i) \end{pmatrix},\quad \text {for }i=1,2. \end{aligned}$$

(2.3)

Here $v_1$ and $v_2$ denote the VEVs of the neutral components of the scalar doublets, and satisfy the relation $v_1^2+v_2^2=v^2\approx (246\text { GeV})^2$. We will further assume that the values of the parameters in the potential ensure that both $v_1\ne 0$ and $v_2\ne 0$ – see case (D) in Ref. [6] for the precise conditions. When this is the case, we can eliminate two parameters from the potential – typically $m_1^2$ and $m_2^2$ – using the minimisation conditions of the potential (i.e. the tadpole equations), which read at tree level

$$\begin{aligned}&\frac{1}{v_1}\frac{\partial V^{(0)}}{\partial h_1}\bigg |_\text {min.}=0=m_1^2-m_3^2\tan \beta +\frac{1}{2} \bigg (\lambda _1c^2_\beta +\lambda _{345}s^2_\beta \bigg )v^2\,, \end{aligned}$$

(2.4)

$$\begin{aligned}&\frac{1}{v_2}\frac{\partial V^{(0)}}{\partial h_2}\bigg |_\text {min.}=0=m_2^2-m_3^2\cot \beta +\frac{1}{2} \bigg (\lambda _2s^2_\beta +\lambda _{345}c^2_\beta \bigg )v^2\,. \end{aligned}$$

(2.5)

The angle $\beta$ in the above two equations is defined from the ratio of VEVs $v_2/v_1\equiv \tan \beta$, and we make use of the following shorthand notations $\lambda _{345}\equiv \lambda _3+\lambda _4+\lambda _5$, $c_x\equiv \cos x$, and $s_x\equiv \sin x$. Once the tadpole equations have been applied, six free parameters remain in the 2HDM Higgs sector, namely

$$\begin{aligned} m_3^2,\,\lambda _{i}\,(i=2,\cdots ,5),\,\tan \beta \,, \end{aligned}$$

(2.6)

where we note that one of the quartic couplings – here we choose $\lambda _1$ – cannot be free as it must be tuned to ensure that one of the CP-even mass eigenstates has a mass of 125 GeV. We will also follow the common choice of trading the off-diagonal mass parameter $m_3$ for a soft $\mathbb {Z}_2$-breaking scale M defined as $M^2\equiv 2m_3^2/s_{2\beta }$. Moreover, while the angle $\beta$ is by definition the angle that rotates away the VEV of one of the two doublets, it is also the angle that diagonalises the CP-odd and charged Higgs mass matrices at tree level. Indeed, applying the rotation matrix $R_\beta$, with

$$\begin{aligned} R_x\equiv \begin{pmatrix} \cos x &{} -\sin x\\ \sin x &{} \cos x \end{pmatrix}\,, \end{aligned}$$

(2.7)

to the component fields of the two doublets by the angle $\beta$, we obtain new states as

$$\begin{aligned} \begin{pmatrix} h_1\\ h_2 \end{pmatrix}=&R_\beta \begin{pmatrix} \phi _1\\ \phi _2 \end{pmatrix}\,,\quad \begin{pmatrix} z_1\\ z_2 \end{pmatrix}=R_\beta \begin{pmatrix} z\\ A \end{pmatrix}\,,\nonumber \\ \begin{pmatrix} w_1^+\\ w_2^+ \end{pmatrix}=&R_\beta \begin{pmatrix} w^+\\ H^+ \end{pmatrix}\,. \end{aligned}$$

(2.8)

In this new basis – often refered to as the Higgs basis – only $\phi _1$ carries a VEV, i.e. $\langle \phi _1\rangle =v$ and $\langle \phi _2\rangle =0$. The fields z, A, $w^+$, and $H^+$ are tree-level mass eigenstates: z and $w^+$ are respectively the neutral and charged would-be Goldstone bosons, while A and $H^+$ are pseudoscalar (i.e. CP-odd) and charged Higgs bosons, with masses (at tree level)

$$\begin{aligned} m_A^2&=M^2-\lambda _5 v^2\,,\nonumber \\ m_{H^\pm }^2&=M^2-\frac{1}{2}(\lambda _4+\lambda _5)v^2\,. \end{aligned}$$

(2.9)

$\phi _1$ and $\phi _2$ are, however, not mass eigenstates in this basis, and their mass matrix reads

$$\begin{aligned} \begin{pmatrix} m_{\phi _1\phi _1}^2 &{} m_{\phi _1\phi _2}^2\\ m_{\phi _1\phi _2}^2 &{} m_{\phi _2\phi _2}^2 \end{pmatrix}\equiv \begin{pmatrix} [\lambda _1c^4_\beta +\lambda _2s^4_\beta +\frac{1}{2}\lambda _{345}s^2_{2\beta }]v^2 &{} \quad -\frac{1}{2}[\lambda _1c^2_\beta -\lambda _2s^2_\beta -\lambda _{345} c_{2\beta }]s_{2\beta }v^2\\ -\frac{1}{2}[\lambda _1c^2_\beta -\lambda _2s^2_\beta -\lambda _{345}c_{2\beta }]s_{2\beta }v^2 &{} \quad M^2 + \frac{1}{4}[\lambda _1+\lambda _2-2\lambda _{345}]s^2_{2\beta }v^2 \end{pmatrix}. \end{aligned}$$

(2.10)

An additional rotation of angle $(\alpha -\beta )$ is necessary to obtain tree-level CP-even mass eigenstates, which will be denoted h and H

$$\begin{aligned} \begin{pmatrix} \phi _1\\ \phi _2 \end{pmatrix}=R_{\alpha -\beta }\begin{pmatrix} H\\ h \end{pmatrix}\,,\quad \text { or equivalently }\begin{pmatrix} h_1\\ h_2 \end{pmatrix}=R_\alpha \begin{pmatrix} H\\ h \end{pmatrix}. \end{aligned}$$

(2.11)

The latter of the two relations shows that $\alpha$ is defined as the CP-even Higgs mixing angle. In turn the tree-level mass eigenvalues for h and H can be found as

$$\begin{aligned} m_H^2=c^2_{\alpha -\beta }m_{\phi _1\phi _1}^2+s_{2(\alpha -\beta )} m_{\phi _1\phi _2}^2+s^2_{\alpha -\beta }m_{\phi _2\phi _2}^2\,,\nonumber \\ m_h^2=s^2_{\alpha -\beta }m_{\phi _1\phi _1}^2-s_{2(\alpha -\beta )} m_{\phi _1\phi _2}^2+c^2_{\alpha -\beta }m_{\phi _2\phi _2}^2\,. \end{aligned}$$

(2.12)

Throughout this paper we will assume that the lightest of the two eigenstates, h, corresponds to the discovered 125-GeV Higgs boson. A limit of particular interest is when the second rotation is not needed to diagonalise the CP-even mass matrix: this is the so-called alignment limit [4], in which the Higgs VEV is aligned in field space with one of the two CP-even mass eigenstates. In terms of mixing angles, two choices are possible to realise this limit, either $s_{\beta -\alpha }=1$ or $c_{\beta -\alpha }=1$ depending on whether h or H is assumed to be the 125-GeV Higgs boson. As we identify the discovered Higgs particle with h, we must require the former condition.

In this limit, the heavy CP-even Higgs mass simplifies to

$$\begin{aligned} m_H^2=M^2+\frac{1}{4}[\lambda _1+\lambda _2-2\lambda _{345}]s^2_{2\beta }v^2\,, \end{aligned}$$

(2.13)

and therefore we find that all the masses of the additional Higgs bosons $\varPhi =H,\,A,\,H^{\pm }$ take the form

$$\begin{aligned} m_\varPhi ^2=M^2+\tilde{\lambda } v^2\,, \end{aligned}$$

(2.14)

where $\tilde{\lambda }$ denotes some simple function of Lagrangian quartic couplings (and $\tan \beta$) – as given in Eqs. (2.9) and (2.13). Moreover, in the alignment limit we can obtain the h-field dependent masses of the additional scalars with the simple replacement $v\rightarrow v+h$, as h is aligned in field space with the VEV v. Similarly, in this limit, the field-dependent mass of the top quark also takes the simple form

$$\begin{aligned} m_t(h)=\frac{y_t}{\sqrt{2}}s_\beta (v+h)\,. \end{aligned}$$

(2.15)

Finally, we should mention that we follow the common choice of trading the five quartic couplings for the four mass eigenvalues $m_h$, $m_H$, $m_A$, and $m_{H^\pm }$, and the CP-even mixing angle $\alpha$ (fixed in our case because we work in the alignment limit). General expressions for this translation are given at tree level for example in Ref. [12], and we only reproduce them here in the limiting case $\alpha =\beta -\pi /2$

$$\begin{aligned} \lambda _1&=\frac{1}{v^2}\big (m_h^2+(m_H^2-M^2)\tan ^2\beta \big )\,,\nonumber \\ \lambda _2&=\frac{1}{v^2}\big (m_h^2+(m_H^2-M^2)\cot ^2\beta \big )\,,\nonumber \\ \lambda _3&=\frac{1}{v^2}\big (m_h^2+2m_{H^\pm }^2-m_H^2-M^2\big )\,,\nonumber \\ \lambda _4&=-\frac{1}{v^2}\big (2m_{H^\pm }^2-m_A^2-M^2\big )\,,\nonumber \\ \lambda _5&=-\frac{1}{v^2}\big (m_A^2-M^2\big )\,. \end{aligned}$$

(2.16)

We should emphasise here that as these are tree-level relations, they can only be used if the masses $m_\varPhi$ are tree-level $\overline{{\mathrm{MS}}}$ mass parameters (the relation between Lagrangian parameters and scalar masses computed at the loop level has been investigated for instance in Refs. [35, 59, 60, 102,103,104,105,106]). Anticipating slightly on the next section’s discussion of our effective-potential calculation, we note that we employ the above relations to express the one- and two-loop contributions to the effective potential in terms of (tree-level) $\overline{{\mathrm{MS}}}$ scalar masses, and once we have taken derivatives of the potential, we add the necessary finite counterterms in order to express our results in terms of physical (i.e. pole) masses.

2.2 The Inert-Doublet Model

The next model we turn to is the Inert-Doublet Model [6, 8] that corresponds to a simple limit of the above 2HDM in which the $\mathbb {Z}_2$ symmetry acts only on one of the two Higgs doublet – say $\varPhi _2$ to fix the notation – and remains unbroken after EWSB. This condition forbids the presence of a mass term $(\varPhi _2^\dagger \varPhi _1+\text {h.c.})$, as well as the appearance of a non-zero VEV for the neutral component of $\varPhi _2$. In this case, the $\mathbb {Z}_2$-odd doublet $\varPhi _2$ cannot mix with the SM-like doublet $\varPhi _1$, nor can it couple to the fermion sector.

We follow the conventions of Ref. [107] and we expand the two scalar doublets as

$$\begin{aligned} \varPhi _1=\begin{pmatrix} G^+\\ \frac{1}{\sqrt{2}}(v+h+iG) \end{pmatrix}\,,\text { and}\quad \varPhi _2=\begin{pmatrix} H^+\\ \frac{1}{\sqrt{2}}(H+iA) \end{pmatrix}\,, \end{aligned}$$

(2.17)

where the notations for the component fields are common with the 2HDM. Because they do not couple to fermions, the components of $\varPhi _2$ – i.e. H, A, $H^\pm$ – are referred to as inert scalars. The lightest of the two neutral of these $\mathbb {Z}_2$-odd states, which we will assume to be H in the following, constitutes a candidate for dark matter (DM) [8, 91].

With the requirements of gauge invariance, the $\mathbb {Z}_2$ symmetry, and assuming once again that there is no new source of CP violation in the Higgs sector, the tree-level scalar potential of the IDM can be written as

$$\begin{aligned} V^{(0)}_\text {IDM}=&\ \mu _1^2|\varPhi _1|^2+\mu _2^2|\varPhi _2|^2+\frac{\lambda _1}{2}|\varPhi _1^2|^4 +\frac{\lambda _2}{2}|\varPhi _2^2|^4\nonumber \\&+\lambda _3|\varPhi _1|^2|\varPhi _2|^2+\lambda _4| \varPhi _1^\dagger \varPhi _2|^2\nonumber \\&+\frac{\lambda _5}{2}\left( (\varPhi _1^\dagger \varPhi _2)^2+\text {h.c.}\right) \,, \end{aligned}$$

(2.18)

where all parameters are real. As only one of the two Higgs doublets acquires a VEV, we have a single tadpole equation

$$\begin{aligned} \frac{1}{v}\frac{\partial V^{(0)}}{\partial h}\bigg |_\text {min.}=0=\mu _1^2+\frac{1}{2}\lambda _1v^2\,, \end{aligned}$$

(2.19)

which we use to eliminate the mass parameter $\mu _1$. We are then left with five free parameters in the Higgs sector, namely

$$\begin{aligned} \mu _2,\,\lambda _i\ (i=2,\cdots ,5)\,, \end{aligned}$$

(2.20)

while v is related to the Fermi constant and $\lambda _1$ is constrained by $m_h=125\text { GeV}$. We note that for the $\mathbb {Z}_2$ symmetry to remain exact after EWSB and the minimum of the potential to correspond to the correct EW minimum, the Lagrangian parameters are constrained, as shown e.g. in case (C) in Ref. [6]. Concurrently, the condition of the potential being bounded from below also gives conditions on the parameters, which are the same as those for the 2HDM given in Eq. (2.2).

We can obtain the tree-level, field-dependent, masses of the inert scalars as

$$\begin{aligned} m_H^2(h)&=\mu _2^2+\lambda _H (v+h)^2\,,\nonumber \\ m_A^2(h)&=\mu _2^2+\lambda _A (v+h)^2\,,\nonumber \\ m_{H^\pm }^2(h)&=\mu _2^2+\lambda _3 (v+h)^2\,, \end{aligned}$$

(2.21)

where $\lambda _{H,A}\equiv \lambda _3+\lambda _4\pm \lambda _5$. It is interesting to note that $\lambda _2$ – the quartic self-coupling of the inert doublet – does not appear in any of these tree-level masses.

As mentioned already, the lightest inert scalar of the IDM – which we assume to be H – is a natural DM candidate, and in our calculations in Sect. 5, we will be considering a particular DM-inspired scenario. One can indeed distinguish [8, 91] two types of scenarios with H as a DM particle: (i) the case where $m_H>m_h$, i.e. all the inert scalars are heavy; or (ii) the case where $m_H\simeq m_h/2$, i.e. the lightest inert scalar is light, while the other two (A and $H^\pm$) can be heavy – this second type of scenario has been discussed in Ref. [62]. On the one hand, for the first case the most natural way to drive the inert scalar masses to high values is to take the mass parameter $\mu _2$ large, but – as we will discuss in detail in the following – this prevents the appearance of large BSM deviations in the Higgs self-couplings. On the other hand, the second type of scenario requires $\mu _2$ to be small to allow $m_H\simeq m_h/2$, and is therefore more interesting from the point of view of obtaining large deviations in the Higgs trilinear and quartic couplings. For this reason, we will in Sect. 5 consider an IDM scenario where $\mu _2=0$ so that $m_H\simeq m_h/2$ while the masses of A and $H^\pm$ can be taken large by increasing the values of the quartic couplings.

2.3 The Higgs-Singlet Model

The third type of model that we consider is an extension of the SM with a real $SU(2)_L$-singlet scalar $\varphi _S$, which we will refer to as “Higgs-Singlet Model” (HSM). Although simple in apparence, the addition of a new singlet scalar can stabilise the Higgs potential [24, 25, 27], and furthermore allows the possibility of a strong first-order EWPT [92]. We expand the SM-like doublet $\varPhi$ and the real singlet $\varphi _S$ as

$$\begin{aligned} \varPhi =\begin{pmatrix} G^+\\ \frac{1}{\sqrt{2}}(v+h+i G)\end{pmatrix}\,,\qquad \text {and}\qquad \varphi _S=v_S+S. \end{aligned}$$

(2.22)

Given the requirement of gauge invariance and assuming that there is no source of CP violation in the HSM scalar sector, the potential of the HSM reads in terms of $\varPhi$ and $\varphi _S$:

$$\begin{aligned} V^{(0)}_\text {HSM}= & {} \mu ^2|\varPhi |^2{+}\frac{1}{2}\mu _S^2\varphi _S^2 +\kappa _1\varphi _S|\varPhi |^2{+}\kappa _2\varphi _S^3+\frac{1}{2}\lambda _H| \varPhi |^4\nonumber \\&+\frac{1}{2}\lambda _{HS}|\varPhi |^2\varphi _S^2+\frac{1}{2}\lambda _S\varphi _S^4\,, \end{aligned}$$

(2.23)

where we have used the freedom to redefine the singlet by a constant shift in order to eliminate the singlet tadpole term that would in principle have been present [92]. The singlet in this model can also become a stable dark matter candidate if we add a $\mathbb {Z}_2$ symmetry under which S changes sign [7]. In this case, the potential reduces to

$$\begin{aligned} V^{(0)}_\text {HSM}= & {} \mu ^2|\varPhi |^2+\frac{1}{2}\mu _S^2S^2+\frac{1}{2} \lambda _H|\varPhi |^4\nonumber \\&+\frac{1}{2}\lambda _{HS}|\varPhi |^2S^2+\frac{1}{2}\lambda _SS^4\,. \end{aligned}$$

(2.24)

Note however that if the sum $\mu _S^2+1/2\lambda _{HS}v^2$ were to be negative, the $\mathbb {Z}_2$ symmetry would be spontaneously broken by a singlet VEV, which also generates again the trilinear couplings $\kappa _1$ and $\kappa _2$ in Eq. (2.23). We will choose to consider such a $\mathbb {Z}_2$-symmetric HSM, ensuring that $\mu _S^2>-1/2\lambda _{HS}v^2$, with the additional motivation⁵ of avoiding mixing between the CP-even component of the Higgs doublet and the singlet.

For $\lambda _{HS}\ge 0$, the HSM tree-level potential is manifestly bounded from below provided that also $\lambda _H$ and $\lambda _S$ are positive. If however $\lambda _{HS}<0$, then one must impose the condition $\lambda _H\lambda _S>1/4 \lambda _{HS}^2$ to avoid the appearance of unstable directions in the potential. Turning next to the counting of parameters in the $\mathbb {Z}_2$-symmetric HSM, one is left with three free parameters, namely

$$\begin{aligned} \mu _S^2,\, \lambda _{HS},\, \lambda _S\,. \end{aligned}$$

(2.25)

Indeed, $\mu ^2$ and $\lambda _H$ can be eliminated respectively with the minimisation condition of the potential and the 125-GeV Higgs mass constraint, while the Higgs VEV v is related to the Fermi constant $G_F$. Moreover, it is also common to trade the quartic coupling $\lambda _{HS}$ for the (tree-level) singlet mass $m_S^2$, using the relation

$$\begin{aligned} m_S^2=\mu _S^2+\frac{1}{2} \lambda _{HS} v^2\,. \end{aligned}$$

(2.26)

Finally, as in the 2HDM and the IDM, the Higgs-field-dependent singlet mass is obtained with the replacement $v\rightarrow v+h$ in the above equation.

Fig. 1

Diagrammatic illustration of the Higgs three-point function $\Gamma _{hhh}$

3 Set-up of the effective potential calculation

3.1 Computation in the $\overline{{\mathrm{MS}}}$ scheme

We investigate in this article the dominant two-loop corrections to the Higgs trilinear and quartic couplings. In principle, the objects that we should compute would then be the three- and four-point functions $\Gamma _{hhh}$ and $\Gamma _{hhhh}$ – the former is shown in Fig. 1. However, these quantities depend on external momenta, making them difficult to compute beyond one loop. Indeed, while closed-form expressions can be obtained for one-loop Passarino-Veltmann functions [109], at two loops one would have to perform numerical integrations (e.g. with SecDec [110]), or derive and solve systems of differential equations between the loop functions – as was done for two-point functions in Ref. [111] (and later implemented in TSIL [112]). This would go well beyond the scope of the present paper in which we limit ourselves to leading two-loop corrections. We therefore choose to neglect the dependence on external momenta and to work in the effective-potential approximation, thereby greatly simplifying the computation. In doing so, we will be missing potential threshold effects, as were found in the complete one-loop calculation in Ref. [12]. However, we can expect the neglected two-loop momentum effects to be subleading, in the light of existing results for scalar mass calculations at two loops – see for instance Refs. [104, 113,114,115,116,117,118], and therefore this setting is sufficient for investigating the possible maximal size of two-loop corrections.

A further approximation that we will make is to neglect contributions from the light scalars, i.e. the 125-GeV Higgs boson and the would-be Goldstone bosons, both at one- and two-loop orders in our calculation. In other words, we will always assume a mass hierarchy of the form

$$\begin{aligned} m_h,\,m_G,\,m_{G^\pm }\ll m_t,\,m_\varPhi , \end{aligned}$$

(3.1)

with $\varPhi$ denoting generically the additional heavy BSM scalars of the 2HDM, IDM, or HSM. We expect that this approximation will not affect our conclusions on the possible size of two-loop contributions from BSM states. Indeed, we know that these contributions grow with increasing BSM-scalar masses, and the large mass limit in which we are interested therefore corresponds precisely to when it is most justified to neglect subleading contributions from light scalars. Moreover, the two-loop diagrams that only involve $h,\,G$, or $G^\pm$ are common with the SM, as we consider here only aligned scenarios, and hence these terms will cancel out from the deviation ratios we will consider in the following. Finally, we should mention that if we choose to include Goldstone contributions, we would encounter infra-red divergences when their running masses become zero or negative – this is the so-called Goldstone Boson Catastrophe [119]. From experience in the case of self-energy calculations [120], we can expect to have to include partial momentum dependence at two loops to solve this technical issue, and we leave this for future work.

We expand the effective potential $V_{{\text {eff}}}$ to successive orders in perturbation theory, up to two loops, as

$$\begin{aligned} V_{{\text {eff}}}\equiv V^{(0)}+\varDelta V_\text {eff}=V^{(0)}+\kappa V^{(1)}+\kappa ^2 V^{(2)}\,, \end{aligned}$$

(3.2)

where $\kappa$ is the loop factor, defined in Eq. (A.1). In terms of $V_{{\text {eff}}}$, we can define effective Higgs trilinear and quartic couplings, and their respective loop expansions, as

$$\begin{aligned} \lambda _{hhh}&\equiv \frac{\partial ^3V_{{\text {eff}}}}{\partial h^3}\bigg |_\text {min}\equiv \lambda _{hhh}^{(0)}+\kappa \delta ^{(1)}\lambda _{hhh}+\kappa ^2 \delta ^{(2)}\lambda _{hhh}\,,\nonumber \\ \lambda _{hhhh}&\equiv \frac{\partial ^4V_{{\text {eff}}}}{\partial h^4}\bigg |_\text {min}\equiv \lambda _{hhhh}^{(0)}+\kappa \delta ^{(1)}\lambda _{hhhh}+\kappa ^2 \delta ^{(2)}\lambda _{hhhh}\,. \end{aligned}$$

(3.3)

Note that these definitions correspond to the following choices of normalisation

$$\begin{aligned} \mathscr {L}\supset -\frac{1}{6}\lambda _{hhh}h^3-\frac{1}{24}\lambda _{hhhh}h^4\,. \end{aligned}$$

(3.4)

Because we are considering scenarios without mixing, we can write the tree-level contributions to $\lambda _{hhh}$ and $\lambda _{hhhh}$ in terms of only the tree-level Higgs mass $m_h$ and Higgs VEV v. It is then convenient to reexpress these in terms of the effective-potential, or curvature, mass of the Higgs, defined as

$$\begin{aligned}{}[M_h^2]_{V_{{\text {eff}}}}&\equiv \mathscr {D}_2 V_{{\text {eff}}}\Big |_\text {min}=m_h^2+\mathscr {D}_2 \varDelta V_{{\text {eff}}}\Big |_\text {min}\,,\qquad \text {where }\nonumber \\ \mathscr {D}_2&\equiv -\frac{1}{v}\frac{\partial }{\partial h}+\frac{\partial ^2}{\partial h^2}\,. \end{aligned}$$

(3.5)

We can then rewrite $\lambda _{hhh}$ as

$$\begin{aligned} \lambda _{hhh}=&\ \frac{3m_h^2}{v}+\frac{\partial ^3\varDelta V_{{\text {eff}}}}{\partial h^3}\bigg |_\text {min}=\frac{3[M_h^2]_{V_{{\text {eff}}}}}{v}+\mathscr {D}_3 \varDelta V_{{\text {eff}}}\Big |_\text {min}\,, \end{aligned}$$

(3.6)

where for the second equality we used Eq. (3.5) and we define

$$\begin{aligned} \mathscr {D}_3\equiv \frac{\partial ^3}{\partial h^3}-\frac{3}{v}\left[ -\frac{1}{v}\frac{\partial }{\partial h}+\frac{\partial ^2}{\partial h^2}\right] \,. \end{aligned}$$

(3.7)

Similarly, we can write

$$\begin{aligned} \lambda _{hhhh}&= \frac{3[M_h^2]_{V_{{\text {eff}}}}}{v^2}+\mathscr {D}_4 \varDelta V_{{\text {eff}}}\Big |_\text {min}\,,\qquad \text {with }\nonumber \\ \mathscr {D}_4&\equiv \frac{\partial ^4}{\partial h^4}-\frac{3}{v^2}\left[ -\frac{1}{v}\frac{\partial }{\partial h}+\frac{\partial ^2}{\partial h^2}\right] \,. \end{aligned}$$

(3.8)

The first derivative term in the definitions of the above differential operators ensures that tadpoles are properly taken into account, by imposing the minimisation condition of the loop corrected potential.

At this point, we should also discuss how renormalisation is performed in this calculation. There are indeed two possible options between which to choose:

(i):

take derivatives of the unrenormalised effective potential, and then perform the renormalisation of the result;

(ii):

renormalise the effective potential first, and take derivatives afterwards.

The two options are of course formally equivalent, but we will prefer here the second one as it conveniently allows us to make use of existing results for two-loop contributions to the effective potential – see e.g. Ref. [121]. These results employ the modified minimal subtraction ($\overline{{\mathrm{MS}}}$) scheme and are expressed in terms of field-dependent tree-level masses. In turn, this implies that the expressions we derive for the Higgs self-couplings using Eq. (3.3) are written in terms of $\overline{{\mathrm{MS}}}$-renormalised parameters.

Before discussing two-loop corrections, we should also review known results for the effective-potential calculation of one-loop corrections to Higgs self-couplings. Using the well-known supertrace formula [122], the dominant one-loop contributions to $V_{{\text {eff}}}$ can be found to be, for the 2HDM, IDM, and HSM

$$\begin{aligned} V^{(1)}= & {} -3m_t^4(h)\left( {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2(h)-\frac{3}{2}\right) \nonumber \\&+\sum _{\varPhi }\frac{n_\varPhi m_\varPhi ^4(h)}{4}\left( {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2(h)-\frac{3}{2}\right) \,, \end{aligned}$$

(3.9)

where the sum on heavy scalars $\varPhi$ includes $\varPhi =H,A,H^\pm$ for the 2HDM, $\varPhi =A,H^\pm$ for the IDM, and $\varPhi =S$ for the HSM. $m_t^2(h)$ and $m_\varPhi ^2(h)$ are the field-dependent masses of the top quark and of the BSM scalars, respectively, and $n_\varPhi =1$ for H and A, and $n_\varPhi =2$ for $H^\pm$. The notation ${{\,\mathrm{\overline{\text {log}}}\,}}x$ is defined in Eq. (A.3). As mentioned above, we have here neglected subleading terms coming from the 125-GeV Higgs and would-be Goldstone bosons. Applying the operators $\mathscr {D}_3$ and $\mathscr {D}_4$, we obtain straightforwardly the leading one-loop corrections to the Higgs trilinear coupling as

$$\begin{aligned} \delta ^{(1)}\lambda _{hhh}={{\,\mathrm{\mathscr {D}_3}\,}}V^{(1)}\Big |_\text {min}= & {} -\frac{48m_t^4}{v^3} +\sum _{\varPhi }\frac{4n_\varPhi m_\varPhi ^4}{v^3}\nonumber \\&\times \bigg (1-\frac{\mathscr {M}^2}{m_\varPhi ^2}\bigg )^3\,, \end{aligned}$$

(3.10)

and for the Higgs quartic coupling,

$$\begin{aligned} \delta ^{(1)}\lambda _{hhhh}=\mathscr {D}_4V^{(1)}\bigg |_\text {min}= & {} -\frac{192m_t^4}{v^4} +\sum _{\varPhi }\frac{8n_\varPhi m_\varPhi ^4}{v^4} \nonumber \\&\times \bigg (1-\frac{\mathscr {M}^2}{m_\varPhi ^2}\bigg )^3 \bigg (2+\frac{\mathscr {M}^2}{m_\varPhi ^2}\bigg )\,.\nonumber \\ \end{aligned}$$

(3.11)

For both equations, we define the shorthand notation $\mathscr {M}$ to denote

$$\begin{aligned} \mathscr {M}=\left\{ \begin{matrix} M\text { for the 2HDM,}\\ \mu _2\text { for the IDM,}\\ \mu _S\text { for the HSM.} \end{matrix} \right. \end{aligned}$$

(3.12)

Fig. 2

Topologies of diagrams contributing to dominant two-loop BSM corrections to the effective potential

Beyond one-loop order, corrections to the effective potential are found not by the supertrace formula, but by computing one-particule-irreducible (1PI) vacuum bubble diagrams [122]. The contributions that we will need in order to investigate the leading BSM effects come from diagrams with only scalars or with scalars and fermions – as shown in Fig. 2. We therefore expand the two-loop potential as

$$\begin{aligned} V^{(2)}=V^{(2)}_{SSS}+V^{(2)}_{SS} + V^{(2)}_{FFS}\,, \end{aligned}$$

(3.13)

where each index S or F indicates a scalar or a Dirac-fermion propagator. Analytic expressions for these, in the $\overline{{\mathrm{MS}}}$ scheme and Landau gauge, can be taken from Ref. [121] (results in a general gauge fixing can be found in Ref. [123]) – note however that as we consider here Dirac instead of Weyl fermions, our definition of $V^{(2)}_{FFS}$ corresponds to the sum of $V_{FFS}$ and $V_{\bar{F}\bar{F}S}$ in Ref. [121]. These expressions only involve the one-loop function A and the two-loop sunrise integral I, and we provide expressions for both in Appendix A. Finally, it should be noted that because we only consider BSM corrections from scalars (neglecting Goldstone bosons) and fermions there is no issue of gauge dependence in our calculations.

3.2 Conversion from $\overline{{\mathrm{MS}}}$ to on-shell renormalisation

Although it makes calculations simple, the $\overline{{\mathrm{MS}}}$ scheme may suffer from a loss of accuracy due to potentially large logarithmic contributions coming from the explicit renormalisation-scale dependence, and hence we choose to convert our calculation to the on-shell (OS) scheme instead. This means that we reexpress our results in terms of physical parameters, namely physical⁶ (or pole) masses and the physical Higgs VEV $v_{\text {phys}}=(\sqrt{2}G_F)^{-1/2}$, and that we moreover must include the effects of finite wave-function renormalisation (WFR). We then obtain OS-renormalised results for the Higgs trilinear and quartic effective couplings, which relate closely to the three- and four-point functions evaluated at vanishing external momenta, as

$$\begin{aligned} \hat{\lambda }_{hhh}&\equiv \left( \frac{Z_h^\text {OS}}{Z_h^{\overline{{\mathrm{MS}}}}} \right) ^{3/2}\lambda _{hhh}=-\Gamma _{hhh}(0,0,0)\,,\nonumber \\ \hat{\lambda }_{hhhh}&\equiv \left( \frac{Z_h^\text {OS}}{Z_h^{\overline{{\mathrm{MS}}}}}\right) ^2 \lambda _{hhhh}=-\Gamma _{hhhh}(0,0,0,0)\,. \end{aligned}$$

(3.14)

In these two equations, $Z_h^\text {OS}$ and $Z_h^{\overline{{\mathrm{MS}}}}$ are respectively the OS- and $\overline{{\mathrm{MS}}}$-scheme Higgs WFR constants, and their ratio can straightforwardly be computed in terms of the corresponding WFR counterterms – $\delta Z_h^\text {OS}$ and $\delta Z_h^{\overline{{\mathrm{MS}}}}$ – as

$$\begin{aligned} \frac{Z_h^\text {OS}}{Z_h^{\overline{{\mathrm{MS}}}}}=1+\delta Z_h^\text {OS}-\delta Z_h^{\overline{{\mathrm{MS}}}} =1+\frac{d}{dp^2}\varPi _{hh}(p^2)\bigg |_{p^2=m_h^2}\,, \end{aligned}$$

(3.15)

where $\varPi _{hh}(p^2)$ is the finite part of the Higgs self-energy, evaluated at external momentum p.

For the masses of the additional scalars and of the top quark, the scheme translation – from $\overline{{\mathrm{MS}}}$ values $m_\varPhi$ ($\varPhi =H,A,H^\pm$ or S) and $m_t$ to physical values $M_\varPhi$ and $M_t$ – also involves the finite part of the corresponding self-energy, i.e.

$$\begin{aligned} M_\varPhi ^2&=m_\varPhi ^2+\varPi _{\varPhi \varPhi }(p^2=M_\varPhi ^2)\quad \text {and}\nonumber \\ M_t^2&=m_t^2+\varPi _{tt}(p^2=M_t^2)\,. \end{aligned}$$

(3.16)

In the case of the 125-GeV Higgs boson, we have already replaced its tree-level mass by its curvature mass, and the latter relates to the physical mass as

$$\begin{aligned} M_h^2=[M_h^2]_{V_{{\text {eff}}}}+\varPi _{hh}(p^2=M_h^2)-\varPi _{hh}(p^2=0)\,. \end{aligned}$$

(3.17)

Finally, the $\overline{{\mathrm{MS}}}$- and OS-renormalised versions of the Higgs VEV satisfy the equation

$$\begin{aligned} v_{\text {phys}}^2=v^2+\kappa \delta ^{(1)}v^2+\kappa ^2\delta ^{(2)}v^2\,. \end{aligned}$$

(3.18)

This general prescription for the scheme conversion is significantly simplified in our case, in particular given that we neglect $m_h,\ m_G,$ and $m_{G^\pm }$ in all loop corrections. First of all, the top quark and BSM scalars only enter the calculation at one loop, and thus we only require one-loop translations for these – the effect of including two-loop corrections to the corresponding self-energies (see Eq. (3.16)) is of three-loop order. Furthermore, as we neglect the 125-GeV-Higgs mass at loop level, we see from Eq. (3.17) that there is then no difference between $M_h^2$ and $[M_h^2]_{V_{{\text {eff}}}}$. Finally, the inclusion of WFR and VEV renormalisation is also simplified: both the multiplication of $\lambda ^{(0)}_{hhh}$ by the finite WFR counterterm, as well as the shift to the Higgs VEV in $\lambda ^{(0)}_{hhh}$ give loop contributions proportional to $M_h^2$ and can therefore be consistently neglected. It is therefore again sufficient to include only one-loop WFR and VEV counterterms here. On the one hand, in the scenarios without mixing that we consider, we find no one-loop correction to the VEV from the BSM scalars, so we only have [19]

$$\begin{aligned} \delta ^{(1)}v^2=-3M_t^2\big (2{{\,\mathrm{\overline{\text {log}}}\,}}M_t^2-1\big )\,. \end{aligned}$$

(3.19)

On the other hand, for the WFR, the additional scalars do give new, model-dependent, contributions – to which we will return in Sect. 5.

4 General $\overline{{\mathrm{MS}}}$ expressions

One may notice from the discussion of models in Sect. 2 that in the scenarios without mixing that we consider, the field-dependent masses always take the form

$$\begin{aligned} m_i^2(h)=\mu _i^2+\frac{1}{2}\hat{\lambda }_i(v+h)^2\,, \end{aligned}$$

(4.1)

where $\mu _i^2$ and $\hat{\lambda }_i$ have respectively mass-dimensions 2 and 0 – $\hat{\lambda }_i$ is either a combination of quartic scalar couplings or a squared Yukawa coupling. This motivates deriving some general expressions for the derivatives of the two-loop integrals contributing to the effective potential, applicable in all scenarios without mixing in the scalar sector. As the potential is renormalised using the $\overline{{\mathrm{MS}}}$ scheme, the results that we derive here are in the same scheme, and a conversion to the OS scheme remains to be done in a model-specific way.

4.1 Eight-shaped diagrams

We consider first the case of the $V^{(2)}_{SS}$ diagrams (see Fig. 2), which involve two scalar propagtors and which we will refer to as eight-shaped diagrams here. They are expressed as [121]

$$\begin{aligned} V^{(2)}_{SS}(m_1^2,m_2^2)=\frac{\lambda _{1122}}{8}A(m_1^2)A(m_2^2)\,, \end{aligned}$$

(4.2)

where $\lambda _{1122}$ is a generic quartic coupling between the two scalars labelled 1 and 2, and $A(m^2)$ is the usual Passarino-Veltmann function [109] (its definition is recalled in Appendix A). Using the differential operators defined in Sect. 3, we obtain in terms of the masses $m_i^2$ and couplings $\hat{\lambda }_i$:

$$\begin{aligned}&\mathscr {D}_2\big [A(m_1^2(h))A(m_2^2(h))\big ]\bigg |_\text {min}\nonumber \\&\quad = \ \hat{\lambda }_1\hat{\lambda }_2v^2+\bigg (\frac{2\hat{\lambda }_1}{m_1^2} +\frac{\hat{\lambda }_2}{m_2^2}\bigg )\hat{\lambda }_2v^2A(m_1^2)\nonumber \\&\qquad +\frac{\hat{\lambda }_1\hat{\lambda }_2}{m_1^2m_2^2}v^2A(m_1^2)A(m_2^2)+(1\leftrightarrow 2)\,,\nonumber \\&\mathscr {D}_3\big [A(m_1^2(h))A(m_2^2(h))\big ]\bigg |_\text {min}\nonumber \\&\quad = \frac{3\hat{\lambda }_1^2\hat{\lambda }_2}{m_1^2}v^3+\bigg (\frac{3 \hat{\lambda }_1}{m_1^2}-\frac{\hat{\lambda }_2}{m_2^2}\bigg ) \frac{\hat{\lambda }_2^2}{m_2^2}v^3A(m_1^2)+(1\leftrightarrow 2)\,,\nonumber \\&\mathscr {D}_4\big [A(m_1^2(h))A(m_2^2(h))\big ]\bigg |_\text {min}\nonumber \\&\quad = \frac{18\hat{\lambda }_1^2\hat{\lambda }_2}{m_1^2}v^2-\frac{4\hat{\lambda }_1^3 \hat{\lambda }_2}{m_1^4}v^4+\frac{3\hat{\lambda }_1^2\hat{\lambda }_2^2}{m_1^2m_2^2}v^4\nonumber \\&\qquad +\bigg [6\hat{\lambda }_2^2v^2\bigg (\frac{3\hat{\lambda }_1}{m_1^2} -\frac{\hat{\lambda }_2}{m_2^2}\bigg )-\frac{2\hat{\lambda }_2^3v^4}{m_2^2} \bigg (\frac{2\hat{\lambda }_1}{m_1^2}-\frac{\hat{\lambda }_2}{m_2^2}\bigg ) \bigg ]\nonumber \\&\qquad \times \,\frac{A(m_1^2)}{m_2^2}+(1\leftrightarrow 2)\,, \end{aligned}$$

(4.3)

where the notation $(1\leftrightarrow 2)$ indicates the permutation of indices 1 and 2.

4.2 Sunrise diagrams

Next, we turn to the sunrise diagrams, corresponding to $V^{(2)}_{SSS}$ and $V^{(2)}_{FFS}$ in Fig. 2. As can be seen for instance in Ref. [121], both types of diagrams are expressed in terms of the sunrise integral I – defined in Eq. (A.9) – as well as products of A functions for $V^{(2)}_{FFS}$. Furthermore, almost all I functions come multiplied by $(v+h)^2$ in models without mixing (because of field-dependent couplings or masses), so we provide here results⁷ for derivatives of the product $(v+h)^2 I(m_1^2(h),m_2^2(h),m_3^2(h))$.

First, for the second derivative we obtain

$$\begin{aligned}&\mathscr {D}_2\big [(v+h)^2I(m_1^2(h),m_2^2(h),m_3^2(h))\big ]\bigg |_{h=0}\nonumber \\&\quad =\ \frac{1}{3}E_1(m_1^2,m_2^2,m_3^2)I(m_1^2,m_2^2,m_3^2)\nonumber \\&\qquad +\frac{1}{3}t_{123}E_1(m_1^2,m_2^2,m_3^2)-E_2(m_1^2,m_2^2,m_3^2)\frac{A(m_1^2)}{m_1^2}\nonumber \\&\qquad -E_3(m_1^2,m_2^2,m_3^2)\frac{A(m_1^2)A(m_2^2)}{m_1^2m_2^2}+(123)\,, \end{aligned}$$

(4.4)

where (123) denotes the sum on cyclical permutations of the indices $\{1,2,3\}$ and

$$\begin{aligned}&E_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \frac{4\hat{\lambda }_1v^2r_{123}}{\varDelta _{123}} -\frac{4J_1(m_1^2,m_2^2,m_3^2)v^4}{\varDelta _{123}^2} +(123)\,,\nonumber \\&E_2(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \frac{4m_1^2v^2(\hat{\lambda }_1r_{123}+(123))}{\varDelta _{123}}+\frac{H_1(m_1^2,m_2^2, m_3^2)v^4}{\varDelta _{123}^2}\,,\nonumber \\&E_3(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \frac{4v^2}{\varDelta _{123}}\big [\hat{\lambda }_1m_2^2r_{231}+\hat{\lambda }_2m_1^2 r_{123}+2\hat{\lambda }_3m_1^2m_2^2\big ]\nonumber \\&\qquad -\frac{2r_{312}v^4}{\varDelta _{123}^2}\big [J_1 (m_1^2,m_2^2,m_3^2)+(123)\big ]\,,\nonumber \\&H_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \chi _{123}\bigg [m_1^2\bigg (\frac{\hat{\lambda }_1^2}{m_1^2}+\frac{\hat{\lambda }_2^2}{m_2^2}+\frac{\hat{\lambda }_3^2}{m_3^2}\bigg )-2\hat{\lambda }_1\big (\hat{\lambda }_2 +\hat{\lambda }_3\big )\bigg ]\nonumber \\&\qquad +4\zeta _{123}\big [\hat{\lambda }_1\hat{\lambda }_2m_3^2 +\hat{\lambda }_1\hat{\lambda }_3m_2^2-\hat{\lambda }_2\hat{\lambda }_3m_1^2\big ]\,, \nonumber \\&J_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \hat{\lambda }_1^2m_2^2m_3^2+\hat{\lambda }_1\hat{\lambda }_2m_3^2r_{312}\,. \end{aligned}$$

(4.5)

Expressions for all the intermediate functions used in these expressions, as well as those in the following, are given in Appendix C.

For the third derivative, we find

$$\begin{aligned}&\mathscr {D}_3\big [(v+h)^2I(m_1^2(h),m_2^2(h),m_3^2(h))\big ]\bigg |_{h=0}\nonumber \\&\quad = -4F_1(m_1^2,m_2^2,m_3^2)I(m_1^2,m_2^2,m_3^2)\nonumber \\&\qquad +6r_{312}F_1(m_1^2,m_2^2,m_3^2) \frac{A(m_1^2)A(m_2^2)}{m_1^2m_2^2}\nonumber \\&\qquad -\bigg [\frac{6H_1(m_1^2,m_2^2,m_3^2)v^3}{\varDelta ^2_{123}} -\frac{H_2(m_1^2,m_2^2,m_3^2)v^5}{\varDelta ^3_{123}}\bigg ]\frac{A(m_1^2)}{m_1^2} \nonumber \\&\qquad -\bigg [\frac{24t_{123}J_1(m_1^2,m_2^2,m_3^2)v^3}{\varDelta _{123}^2} +\frac{J_2(m_1^2,m_2^2,m_3^2)v^5}{\varDelta ^3_{123}}\bigg ]\nonumber \\&\qquad +(123)\,, \end{aligned}$$

(4.6)

where

$$\begin{aligned}&F_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \frac{2J_1(m_1^2,m_2^2,m_3^2)v^3}{\varDelta ^2_{123}} -\frac{L_1(m_1^2,m_2^2,m_3^2)v^5}{\varDelta ^3_{123}}+(123)\,,\nonumber \\&H_2(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv m_1^2\bigg [\varTheta _{123}\frac{\hat{\lambda }_1^3}{m_1^2}+\theta _{123}\frac{\hat{\lambda }_2^3}{m_2^4}+\theta _{132} \frac{\hat{\lambda }_3^3}{m_3^4}\bigg ]+12\rho _{123}\hat{\lambda }_1\hat{\lambda }_2 \hat{\lambda }_3\nonumber \\&\qquad -3m_1^2\bigg [\mu _{123}\hat{\lambda }_2\bigg (\frac{\hat{\lambda }_1^2}{m_1^2} + \frac{\hat{\lambda }_3^2}{m_3^2} \bigg ) +\mu _{132}\hat{\lambda }_3 \bigg (\frac{\hat{\lambda }_1^2}{m_1^2}+ \frac{\hat{\lambda }_2^2}{m_2^2} \bigg )\nonumber \\&\qquad -\nu _{123}\hat{\lambda }_1\bigg (\frac{\hat{\lambda }_2^2}{m_2^2} +\frac{\hat{\lambda }_3^2}{m_3^2}\bigg )\bigg ]\,,\nonumber \\&J_2(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \frac{\hat{\lambda }_1^3}{m_1^2}\tau _{123}r_{123} -16\phi _{123}\hat{\lambda }_1\hat{\lambda }_2\hat{\lambda }_3 \nonumber \\&\qquad -3\big [\varPhi _{123} \hat{\lambda }_1^2\hat{\lambda }_2+\varPhi _{213}\hat{\lambda }_2^2\hat{\lambda }_1\big ]\,,\nonumber \\&L_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \hat{\lambda }_1^3m_2^2m_3^2r_{123} -\hat{\lambda }_1\hat{\lambda }_2m_3^2(\hat{\lambda }_1\omega _{123} +\hat{\lambda }_2\omega _{213})\nonumber \\&\qquad +\frac{1}{3}\varXi _{123}\hat{\lambda }_1 \hat{\lambda }_2\hat{\lambda }_3\,. \end{aligned}$$

(4.7)

For the fourth derivative we have

$$\begin{aligned}&\mathscr {D}_4\big [(v+h)^2I(m_1^2(h),m_2^2(h),m_3^2(h))\big ]\bigg |_{h=0}\nonumber \\&\quad =-\frac{2}{3}G_1(m_1^2,m_2^2,m_3^2)I(m_1^2,m_2^2,m_3^2)\nonumber \\&\qquad +r_{312}G_1 (m_1^2,m_2^2,m_3^2)\frac{A(m_1^2)A(m_2^2)}{m_1^2m_2^2}\nonumber \\&\qquad -\bigg [\frac{36H_1(m_1^2,m_2^2,m_3^2)v^2}{\varDelta _{123}^2} -\frac{14H_2(m_1^2,m_2^2,m_3^2)v^4}{\varDelta _{123}^3}\nonumber \\&\qquad +\frac{2H_3(m_1^2,m_2^2,m_3^2) v^6}{\varDelta _{123}^4}\bigg ]\frac{A(m_1^2)}{m_1^2}\nonumber \\&\qquad -\bigg [\frac{144t_{123}J_1(m_1^2,m_2^2,m_3^2)v^2}{\varDelta _{123}^2} +\frac{14J_2(m_1^2,m_2^2,m_3^2)v^4}{\varDelta _{123}^3}\nonumber \\&\qquad +\frac{J_3(m_1^2,m_2^2,m_3^2) v^6}{\varDelta _{123}^4}\bigg ]+(123)\,, \end{aligned}$$

(4.8)

where

$$\begin{aligned}&G_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \frac{72v^2}{\varDelta _{123}^2}J_1(m_1^2,m_2^2,m_3^2) -\frac{84v^4}{\varDelta _{123}^3}L_1(m_1^2,m_2^2,m_3^2)\nonumber \\&\qquad +\frac{24v^6}{\varDelta _{123}^4}\Big [\hat{\lambda }_1^4m_2^2m_3^2 (\omega _{123}-3m_2^2r_{123}-m_2^2m_3^2)\nonumber \\&\qquad -\hat{\lambda }_1\hat{\lambda }_2m_3^2 (\hat{\lambda }_1^2\xi _{123}+\hat{\lambda }_2^2\xi _{213})\nonumber \\&\qquad +3\hat{\lambda }_1^2\hat{\lambda }_2^2m_3^2(\varXi _{123}-2m_1^2m_2^2m_3^2) {+}\hat{\lambda }_1^2\hat{\lambda }_2\hat{\lambda }_3a_{123}\Big ]{+}(123)\,,\nonumber \\&H_3(m_1^2,m_2^2,m_3^2)\nonumber \\&\equiv m_1^2\big [n_{123}\hat{\lambda }_1^4+p_{123}\hat{\lambda }_2^4 +p_{132}\hat{\lambda }_3^4\big ]\nonumber \\&\qquad +6m_1^4w_{123} \bigg [\frac{\hat{\lambda }_1^2 \hat{\lambda }_2^2}{m_1^2m_2^2}+\frac{\hat{\lambda }_1^2\hat{\lambda }_3^2}{m_1^2m_3^2} +\frac{\hat{\lambda }_2^2\hat{\lambda }_3^2}{m_2^2m_3^2}\bigg ]\nonumber \\&\qquad -2m_1^2\big [2\hat{\lambda }_1^3(\hat{\lambda }_2q_{123}+\hat{\lambda }_3 q_{132})+\hat{\lambda }_1(\hat{\lambda }_2^3u_{123}+\hat{\lambda }_3^3u_{132})\nonumber \\&\qquad +\hat{\lambda }_2\hat{\lambda }_3(\hat{\lambda }_2^2v_{123}+\hat{\lambda }_3^2v_{132}) \big ]\nonumber \\&\qquad +6m_1^2\hat{\lambda }_1\hat{\lambda }_2\hat{\lambda }_3(6\hat{\lambda }_1 A_{123}-\hat{\lambda }_2B_{123}-\hat{\lambda }_3B_{132})\,,\nonumber \\&J_3(m_1^2,m_2^2,m_3^2)\nonumber \\&\quad \equiv \hat{\lambda }_1^4e_{123}+4\hat{\lambda }_1\hat{\lambda }_2 (\hat{\lambda }_1^2f_{123}+\hat{\lambda }_2^2f_{213})-\frac{6\hat{\lambda }_1^2 \hat{\lambda }_2^2}{m_1^2m_2^2}g_{123}\nonumber \\&\qquad -24\hat{\lambda }_1^2\hat{\lambda }_2 \hat{\lambda }_3h_{123}\,. \end{aligned}$$

(4.9)

5 Analytic results for the leading two-loop corrections

In this section, we describe the details of the calculations of leading two-loop corrections to $\hat{\lambda }_{hhh}$ and $\hat{\lambda }_{hhhh}$ in the different models we consider. Analytic expressions for the one-loop effects have been given already in Sect. 3.

5.1 Standard model

We begin with a detailed presentation of the calculation of the dominant SM contributions at two loops – these have been shown already in Refs. [85, 86]. The two dominant contributions to the two-loop SM effective potential can be taken from e.g. Ref. [124], and read

$$\begin{aligned} V^{(2)}_\text {SM}(h)=&\ V^{(2)}_{ttg}(h)+V^{(2)}_{tt\phi }(h)+\cdots \,,\nonumber \\ V^{(2)}_{ttg}(h)=&\ -4g_3^2m_t^2(h)\bigg [4A(m_t^2(h))-8m_t^2(h)-\frac{6A(m_t^2(h))^2}{m_t^2(h)}\bigg ]\,,\nonumber \\ V^{(2)}_{tt\phi }(h)=&\ 3 y_t^2 \bigg [(2 m_t^2(h) - m_h^2(h)) I (m_t^2(h), m_t^2(h), m_h^2(h))\nonumber \\&- \frac{1}{2} m_G^2(h) I (m_t^2(h), m_t^2(h), m_G^2(h)) \nonumber \\&+ (m_t^2(h) - m_G^2(h)) I (m_t^2(h), m_G^2(h), 0)\nonumber \\&+ A (m_t^2(h))^2 - A (m_t^2(h)) A (m_h^2(h))\nonumber \\&- 2 A (m_t^2(h)) A (m_G^2(h))\bigg ]\nonumber \\ \underset{m_h,m_G\ll m_t}{\rightarrow }&\ 3 y_t^2 \bigg [2 m_t^2(h) I (m_t^2(h), m_t^2(h), 0) \nonumber \\&\quad + m_t^2(h) I (m_t^2(h), 0, 0) + A (m_t^2(h))^2\bigg ]\,. \end{aligned}$$

(5.1)

where $g_3$ denotes the $SU(3)_C$ gauge coupling. The dominant contributions to the Higgs trilinear couplings, expressed in terms of $\overline{{\mathrm{MS}}}$ parameters, are then given by

$$\begin{aligned} \delta ^{(2)}\lambda _{hhh}=&\ \mathscr {D}_3\bigg [V^{(2)}_{ttg}(h)+V^{(2)}_{tt\phi }(h)\bigg ]\bigg |_\text {min}\nonumber \\ =&\ \frac{3m_h^2}{v}\bigg [\frac{128 g_3^2 m_t^4 (1 + 6 {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2)}{3 m_h^2 v^2}\nonumber \\&\quad -\frac{8 m_t^4 y_t^2 (-7 + 6 {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2)}{m_h^2 v^2}\bigg ]\,, \end{aligned}$$

(5.2)

which corresponds to Eq. (11) in Ref. [85]. For the quartic coupling we find

$$\begin{aligned} \delta ^{(2)}\lambda _{hhhh}=&\ \mathscr {D}_4\bigg [V^{(2)}_{ttg}(h)+V^{(2)}_{tt\phi }(h)\bigg ]\bigg |_\text {min}\nonumber \\ =&\ \frac{3m_h^2}{v^2}\bigg [\frac{1024 g_3^2 m_t^4 (2 + 3 {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2)}{3m_h^2v^2}\nonumber \\&-\frac{64 y_t^2m_t^4 (-2 + 3 {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2)}{m_h^2v^2}\bigg ]\,. \end{aligned}$$

(5.3)

Fig. 3

Comparison of the dependence on the renormalisation scale Q of our Standard-Model results in the $\overline{{\mathrm{MS}}}$ and OS schemes, at both one- and two-loop orders. Here the quantity $\varDelta \lambda _{hhh}^\text {SM}$ is defined as $\lambda _{hhh}^\text {SM}-(\lambda _{hhh}^{(0)})^\text {SM}$. The numerical inputs used for SM parameters in this paper are as follows: for $g_3$ and $v_{\text {phys}}$ we use values from the PDG [125], respectively $\alpha _S^{\overline{{\mathrm{MS}}}} (Q=M_Z)=g_3^2/4\pi =0.1181$ and $G_F=1/\sqrt{2}v_{\text {phys}}^2=1.1663787\cdot 10^{-5}\text { GeV}^{-2}$, while for the top quark pole mass we take $M_t=173.5\text { GeV}$

These results can be reexpressed straightforwardly in terms of on-shell scheme quantities – $M_t$ and $v_{\text {phys}}$. This conversion only requires (i) the shift to the Higgs VEV given in Eq. (3.19); (ii) the one-loop top quark self-energy in the SM

$$\begin{aligned} \varPi ^{(1)}_{tt}(p^2=M_t^2)\simeq & {} M_t^2\bigg [\frac{4}{3}g_3^2(8-6{{\,\mathrm{\overline{\text {log}}}\,}}M_t^2)\nonumber \\&+\frac{M_t^2}{v_{\text {phys}}^2}(-8+3{{\,\mathrm{\overline{\text {log}}}\,}}M_t^2)\bigg ]\,, \end{aligned}$$

(5.4)

where we have taken the limit $M_h\ll M_t$; and (iii) the derivative with respect to the momentum of the one-loop Higgs self-energy, namely

$$\begin{aligned} \frac{d}{dp^2}\varPi ^{(1)}_{hh}(p^2)\bigg |_{p^2=M_h^2}=\frac{6M_t^2}{v_{\text {phys}}^2}\bigg ({{\,\mathrm{\overline{\text {log}}}\,}}M_t^2+\frac{2}{3}\bigg )\,, \end{aligned}$$

(5.5)

where we have taken the same limit as for $\varPi ^{(1)}_{tt}$. Note also that as the gauge coupling $g_3$ only appears in the corrections to the Higgs self-couplings at two loops, we do not need to specify its renormalisation scheme here. We obtain finally

$$\begin{aligned} \delta ^{(2)}\hat{\lambda }_{hhh} =&\ \frac{72M_t^4}{v_{\text {phys}}^3}\bigg (16g_3^2-\frac{13M_t^2}{v_{\text {phys}}^2}\bigg )\,,\nonumber \\ \delta ^{(2)}\hat{\lambda }_{hhhh}=&\ \frac{384M_t^4}{v_{\text {phys}}^4}\bigg (16g_3^2-\frac{13M_t^2}{v_{\text {phys}}^2}\bigg )\,. \end{aligned}$$

(5.6)

In Fig. 3, we compare the different results that we obtain for the Higgs trilinear coupling, both in the $\overline{{\mathrm{MS}}}$-scheme (dashed curves) and the OS-scheme (solid curves) at one loop (blue lines) and at two loops (red lines). For the $\overline{{\mathrm{MS}}}$ results, we include in this figure only the explicit Q dependence coming from logarithmic terms, and not the running of renormalisation group equations. We can observe, satisfactorily, that this explicit renormalisation scale dependence is reduced when going from one- to two-loop order. Furthermore, if we compare the two-loop $\overline{{\mathrm{MS}}}$ and OS results for $Q=m_t$ – which is the most natural choice of renormalisation scale for the $\overline{{\mathrm{MS}}}$ expression – we find that the two values are extremely close. This provides an important cross-check of our calculation and scheme conversion, while the difference between the results in the two schemes for varying Q provides a rough estimate⁸ of the missing higher-order strong corrections.

5.2 Aligned scenario of a Two-Higgs-Doublet model

The first BSM model that we consider is an aligned scenario of a 2HDM. As discussed in Sect. 2.1, requiring alignment – in other words fixing $\alpha =\beta -\pi /2$ – allows to evade experimental constraints more easily, and on the technical side means we can avoid the complications due to the mixing of the CP-even h and H. In principle, this alignment condition is a tree-level relation and it receives radiative corrections that should be taken into account. However these corrections were studied e.g. in Ref. [104] and were found to be typically very small, so that we will neglect them throughout this work. It should be noted, moreover, that in the presence of four different mass scales – $M_H$, $M_A$, $M_{H^\pm }$, and $\tilde{M}$ – for the BSM scalars, plus the top quark mass $M_t$, the expressions of the radiative corrections at two loops become quite long and cumbersome. We therefore choose to provide complete results in Appendix B.2, and for the main text of this paper we will restrict ourselves to taking the masses of H, A, and $H^\pm$ to be equal. This reduces the number of mass scales and thus allows more compact expressions, without missing any important physical behaviour.

After taking equal the three additional scalar masses, the BSM contributions to the two-loop effective potential of the 2HDM read

$$\begin{aligned} V^{(2)}_{SSS}(h)=&-\frac{4(M^2 - m_\varPhi ^2)^2(v+h)^2}{v^4}I(0,m_\varPhi ^2(h),m_\varPhi ^2(h))\nonumber \\&-\frac{6(M^2 - m_\varPhi ^2)^2\cot ^22\beta (v+h)^2}{v^4}\nonumber \\&\times I(m_\varPhi ^2(h),m_\varPhi ^2(h),m_\varPhi ^2(h))\,,\nonumber \\ V^{(2)}_{SS}(h)=&-\frac{12(M^2 - m_\varPhi ^2)\cot ^22\beta }{v^2}A(m_\varPhi ^2(h))^2\,,\nonumber \\ V^{(2)}_{FFS}(h)=&\ 3 y_t^2c_\beta ^2\Big [A(m_t^2(h))^2 - 3 A(m_\varPhi ^2(h)) A(m_t^2(h))\nonumber \\&+ (m_t^2(h) - m_\varPhi ^2(h)) I(0, m_\varPhi ^2(h),m_t^2(h))\nonumber \\&+ (2m_t^2(h) - m_\varPhi ^2(h)) I(m_\varPhi ^2(h), m_t^2(h), m_t^2(h)) \Big ]\,. \end{aligned}$$

(5.7)

Note that in these expressions we have indicated explicitly which masses should be understood as field-dependent and which should not. The corresponding Feynman diagrams are shown in Fig. 4.

Fig. 4

Dominant two-loop diagrams contributing to the 2HDM effective potential, in the limit of degenerate BSM scalar masses

$\overline{{\mathrm{MS}}}$expressions – Applying then the operators $\mathscr {D}_3$ and $\mathscr {D}_4$, we obtain $\overline{{\mathrm{MS}}}$ expressions for the leading two-loop corrections to the Higgs trilinear and quartic couplings

$$\begin{aligned} \delta ^{(2)}\lambda _{hhh}=&\ \frac{16 m_\varPhi ^4}{v^5}\left( 4+9\cot ^22\beta \right) \left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^4\nonumber \\&\times \left[ -2 M^2 - m_\varPhi ^2 + (M^2 + 2 m_\varPhi ^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\right] \nonumber \\&+\frac{192 m_\varPhi ^6 \cot ^22\beta }{v^5} \left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^4 \left[ 1+2{{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\right] \nonumber \\&+\frac{96 m_\varPhi ^4m_t^2 \cot ^2\beta }{ v^5}\left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^3\nonumber \\&\times \big [-1 + 2 {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2 \big ]+\mathscr {O}\left( \frac{m_\varPhi ^2m_t^4}{v^5}\right) \,, \end{aligned}$$

(5.8)

$$\begin{aligned} \delta ^{(2)}\lambda _{hhhh}=&\ \frac{32m_\varPhi ^2}{v^6}\big (4+9\cot ^22\beta \big ) \left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^4\nonumber \\&\times \big [-5 M^4 - 4 M^2 m_\varPhi ^2 \nonumber \\&+ (2 M^4 + 3 M^2 m_\varPhi ^2 + 4 m_\varPhi ^4) {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\big ]\nonumber \\&+\frac{384m_\varPhi ^6 \cot ^22\beta }{v^6} \left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^4\nonumber \\&\times \big [-M^2 + 4 m_\varPhi ^2 + 2 (M^2 + 2 m_\varPhi ^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\big ]\nonumber \\&+\frac{192m_t^2m_\varPhi ^2\cot ^2\beta }{v^6}\left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^3\nonumber \\&\times \big [-3M^2+2(M^2+2m_\varPhi ^2){{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\big ]\nonumber \\&+\mathscr {O}\left( \frac{m_\varPhi ^2m_t^4}{v^6}\right) \,. \end{aligned}$$

(5.9)

The full expressions for the derivatives of the $V^{(2)}_{FFS}$ diagrams are quite long, so we here give only the leading results – in the third lines of Eqs. (5.8) and (5.9) – and we provide the complete results in Appendix B. The first and second terms of these two equations come respectively from derivatives of sunrise and eight-shaped scalar diagrams in the effective potential – see Fig. 4.

There are several checks that can be performed on these results to verify their validity. First, we have verified that the dependence of the expressions on the renormalisation scale (appearing as $\log Q^2$) is correctly cancelled when including the running of all parameters at one loop – we will return to this when discussing the scheme conversion. Furthermore, the BSM corrections should decouple when taking the mass of the additional scalars to large values, because of the decoupling theorem [126]. More precisely, the expressions in Eqs. (5.8) and (5.9) should tend to zero in the limit $m_\varPhi \rightarrow \infty$. Given that the scalar masses all satisfy a relation of the form $m_\varPhi ^2=M^2+\tilde{\lambda } v^2$, one could in principle achieve the previously-mentioned limit by taking either $M^2$ or $\tilde{\lambda } v^2$ to infinity. However, the latter option would cause a breakdown of perturbativity, and hence the decoupling limit can only be taken properly with the limit $M\rightarrow \infty$, while keeping $\tilde{\lambda }$ fixed. It is then straightforward to see that the above results for the corrections to $\lambda _{hhh}$ and $\lambda _{hhhh}$ decouple properly as all the terms involved are of the form, with $m<n$

$$\begin{aligned} (m_\varPhi ^2)^m \left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^n=\frac{(\tilde{\lambda } v^2)^n}{(M^2+\tilde{\lambda } v^2)^{n-m}}\underset{M\rightarrow \infty }{\xrightarrow {\quad m<n\quad }}0\,. \end{aligned}$$

(5.10)

Conversion to the on-shell scheme – As explained in Sect. 3, we prefer to express our results in terms of pole masses and of the physical Higgs VEV, and hence we convert the expressions in Eqs. (5.8) and (5.9) from the $\overline{{\mathrm{MS}}}$ to the OS scheme. Because we neglect loop corrections proportional to the lightest Higgs mass, it suffices here to translate the parameters that appear in the one-loop corrections – c.f. Eqs. (3.10) and (3.11) – namely the BSM scalar masses $m_\varPhi$, the top-quark mass $m_t$, the Higgs VEV v, and the soft-breaking scale of the 2HDM $\mathbb {Z}_2$ symmetry M. We will return to the discussion of M and its conversion in detail in the following, while for the Higgs VEV, the SM result in Eq. (3.19) is enough as we work in an aligned 2HDM scenario. Furthermore, we must also include finite WFR for the Higgs bosons on the external legs, following Eq. (3.14). The parameter $\tan \beta$ only appears at two loops in our calculation – once again because we work in an aligned scenario – and therefore we will not need a scheme conversion for it here.

The first intermediate results we require for the scheme conversion are the one-loop self-energies of the additional 2HDM scalars, up to leading order in powers of $m_t$. These read

$$\begin{aligned} \varPi _{HH}(p^2)=&-\frac{2(M^2-m_H^2)}{v^2}\cot ^22\beta \nonumber \\&\times \Big [3A(m_H^2)+A(m_A^2) +2A(m_{H^\pm }^2)\Big ]\nonumber \\&-\frac{4(M^2-m_H^2)^2}{v^2}B_0(p^2,0,m_H^2)\nonumber \\&-\frac{18(M^2-m_H^2)^2}{v^2} \cot ^22\beta B_0(p^2,m_H^2,m_H^2)\nonumber \\&-\frac{(m_H^2-m_A^2)^2}{v^2}B_0(p^2,0,m_A^2)\nonumber \\&-\frac{2(m_H^2-m_{H^\pm }^2)^2}{v^2} B_0(p^2,0,m_{H^\pm }^2)\nonumber \\&-\frac{2(M^2-m_H^2)^2}{v^2}\cot ^22\beta \big [B_0(p^2,m_A^2,m_A^2)\nonumber \\&+2B_0(p^2,m_{H^\pm }^2,m_{H^\pm }^2)\big ]\nonumber \\&-\frac{12m_t^2\cot ^2\beta }{v^2}\nonumber \\&\times \bigg [A(m_t^2)-\bigg (2m_t^2- \frac{p^2}{2}\bigg ) B_0(p^2,m_t^2,m_t^2)\bigg ]\,,\nonumber \\ \varPi _{AA}(p^2) =&-\frac{2(M^2-m_H^2)}{v^2}\cot ^22\beta \nonumber \\&\times \Big [A(m_H^2)+3A(m_A^2) +2A(m_{H^\pm }^2)\Big ]\nonumber \\&-\frac{4(M^2-m_A^2)^2}{v^2}B_0(p^2,0,m_A^2)\nonumber \\&-\frac{4(M^2-m_H^2)^2}{v^2}\cot ^22 \beta B_0(p^2,m_A^2,m_H^2)\nonumber \\&-\frac{(m_A^2-m_H^2)^2}{v^2}B_0(p^2,0,m_H^2)\nonumber \\&-\frac{2(m_A^2-m_{H^\pm }^2)^2}{v^2} B_0(p^2,0,m_{H^\pm }^2)\nonumber \\&-\frac{12m_t^2\cot ^2\beta }{v^2}\bigg [A(m_t^2)+ \frac{p^2}{2} B_0(p^2,m_t^2,m_t^2)\bigg ]\,,\nonumber \end{aligned}$$

$$\begin{aligned} \varPi _{H^+H^-}(p^2)=&-\frac{2(M^2-m_H^2)}{v^2}\cot ^22\beta \nonumber \\&\times \Big [A(m_H^2)+A(m_A^2) +4A(m_{H^\pm }^2)\Big ]\nonumber \\&-\frac{4(M^2-m_{H^\pm }^2)^2}{v^2}B_0(p^2,0,m_{H^\pm }^2)\nonumber \\&-\frac{4(M^2-m_H^2)^2}{v^2}\cot ^22\beta B_0(p^2,m_{H^\pm }^2,m_H^2)\nonumber \\&-\frac{(m_{H^\pm }^2-m_H^2)^2}{v^2}B_0(p^2,0,m_H^2)\nonumber \\&-\frac{(m_{H^\pm }^2-m_A^2)^2}{v^2}B_0(p^2,0,m_A^2)\nonumber \\&-\frac{6m_t^2\cot ^2\beta }{v^2}\big [A(m_t^2)+(p^2-m_t^2)\nonumber \\&\quad \times B_0(p^2,0,m_t^2)\big ]\,. \end{aligned}$$

(5.11)

It is important to remember when converting the BSM scalar masses that the parameter points giving either equal tree-level ($\overline{{\mathrm{MS}}}$) masses or equal pole (OS) masses are distinct. Indeed as can be seen from the equations above, the self-energy of H is different from those of A and $H^\pm$, even in the limit of equal masses. Keeping this in mind, we choose nevertheless to show our results in the OS renormalisation scheme also in the limit of equal pole masses.

Turning next to the top quark, although it is required in principle, we do not need the one-loop top-quark self-energy here as we restrict ourselves to terms of order $M_t^2$ only. This is because the dominant BSM corrections to the top-quark self-energy involve the top-quark Yukawa coupling, and are thus proportional to $M_t^2$, and in turn shifting the top-quark mass yields terms of order $M_t^4$. For completeness, we provide the expression of the top-quark self-energy in Eq. (B.3).

At this point, the case of the soft $\mathbb {Z}_2$-symmetry breaking scale M deserves closer attention. This parameter M is not directly related to any physical observable (unlike the scalar and top-quark masses, or the Higgs VEV), and there is thus no clear way to define an on-shell prescription for it. For this reason, one may at first think that there is no use to convert M, and that one may simply continue using its $\overline{{\mathrm{MS}}}$ value. However, the question of the proper decoupling of BSM corrections enters again here, and provides motivation to devise a new “on-shell” prescription for M. Indeed, no matter the renormalisation scheme in which they are expressed, the BSM contributions must vanish in the limit of large BSM scalar masses. In the previous section we found that this is the case for the results written in terms of $\overline{{\mathrm{MS}}}$ parameters, with $m_\varPhi ^2=M^2+\tilde{\lambda } v^2$, when we take the limit $M\rightarrow \infty$. If instead the corrections to the Higgs self-couplings are expressed in terms of pole (i.e. OS-renormalised) masses $M_\varPhi$ and of the soft-breaking scale M still in the $\overline{{\mathrm{MS}}}$ scheme, we must use a one-loop relation between $M_\varPhi$ and M to verify the decoupling behaviour – otherwise, part of the two-loop effects that ensure the decoupling are missed. In practice, this corresponds to having expressions with the additional scalar masses renormalised in the $\overline{{\mathrm{MS}}}$ scheme, but the top-quark mass and Higgs VEV in the OS scheme, and one can straightforwardly find that decoupling is satisfied in this case.

While the need to use a one-loop relation between $M_\varPhi ^2$ and $M^2$ poses no problem, it constitutes a good motivation to define an “on-shell” prescription for M – note that we use here inverted commas for “on-shell” as we are not actually relating M to a physical observable in our prescription. The new quantity that we obtain in this manner, and which we will denote $\tilde{M}$, should be interpreted as the OS-renormalised version of the soft breaking scale of the $\mathbb {Z}_2$ symmetry in the 2HDM. It is, by construction, the parameter that controls the possibility of decoupling of the additional BSM scalars in the 2HDM, when working with all other parameters in the OS scheme. We relate the new $\tilde{M}$ to its $\overline{{\mathrm{MS}}}$ counterpart M with a finite counterterm denoted $\delta ^\text {OS}M^2$ as

$$\begin{aligned} \tilde{M}^2=M^2+\delta ^\text {OS} M^2\,, \end{aligned}$$

(5.12)

and we define this finite counterterm from the requirement that the decoupling behaviour of the BSM corrections to the Higgs trilinear coupling should be apparent when using a relation of the form $M_\varPhi ^2=\tilde{M}^2+\tilde{\lambda } v^2$. With this prescription, we obtain up to one-loop order

$$\begin{aligned} \delta ^\text {OS}M^2&=\frac{3\kappa M^2}{v^2}\bigg [4\cot ^22\beta \left( 1-\frac{M^2}{m_\varPhi ^2}\right) A(m_\varPhi ^2)\nonumber \\&\quad -m_t^2\cot ^2\beta \Big [B_0(m_\varPhi ^2,m_t^2,m_t^2)+B_0(m_\varPhi ^2,0,m_t^2) \Big ]\bigg ]. \end{aligned}$$

(5.13)

Finally, to include the 125-GeV Higgs WFR, we further need the derivative with respect to momentum of the one-loop Higgs self-energy. In the 2HDM, we find for it

$$\begin{aligned} \frac{d}{dp^2}\varPi ^{(1)}_{hh}\Big |_{p^2=0}= & {} \frac{6M_t^2}{v_{\text {phys}}^2}\bigg ({{\,\mathrm{\overline{\text {log}}}\,}}M_t^2+\frac{2}{3}\bigg )\nonumber \\&-\sum _{\varPhi =H,A,H^\pm }\frac{n_\varPhi M_\varPhi ^2}{2v_{\text {phys}}^2}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^2.\nonumber \\ \end{aligned}$$

(5.14)

Combining all these intermediate results, we obtain expressions for the two-loop corrections to the Higgs trilinear and quartic couplings in terms of physical quantities as

$$\begin{aligned} \delta ^{(2)}\hat{\lambda }_{hhh}=&\ \frac{48 M_\varPhi ^6}{v_{\text {phys}}^5}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^4\nonumber \\&\times \left\{ 4+3\cot ^22\beta \left[ 3-\frac{\pi }{\sqrt{3}}\left( \frac{\tilde{M}^2}{M_\varPhi ^2} + 2\right) \right] \right\} \nonumber \\&+\frac{576 M_\varPhi ^6 \cot ^22\beta }{v_{\text {phys}}^5} \left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^4\nonumber \\&+\frac{288 M_\varPhi ^4M_t^2 \cot ^2\beta }{v_{\text {phys}}^5}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^3\nonumber \\&-\frac{48M_\varPhi ^6}{v_{\text {phys}}^5}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^5 +\frac{168M_\varPhi ^4M_t^2}{v_{\text {phys}}^5}\nonumber \\&\times \left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2} \right) ^3+\mathscr {O}\left( \frac{M_\varPhi ^2M_t^4}{v_{\text {phys}}^5}\right) \,, \end{aligned}$$

(5.15)

and,

$$\begin{aligned} \delta ^{(2)}\hat{\lambda }_{hhhh}=&\ \frac{128 M_\varPhi ^6}{v_{\text {phys}}^6}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^4 \bigg [8 + \frac{2 \tilde{M}^2}{M_\varPhi ^2} -\frac{\tilde{M}^4}{M_\varPhi ^4}\bigg ]\nonumber \\&+\frac{576 M_\varPhi ^6\cot ^22\beta }{v_{\text {phys}}^6}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2} \right) ^4\nonumber \\&\times \bigg \{4 + \frac{\tilde{M}^2}{M_\varPhi ^2} - \frac{\tilde{M}^4}{2M_\varPhi ^4} - \frac{\pi }{\sqrt{3}}\bigg [2+\frac{3}{2} \frac{\tilde{M}^2}{M_\varPhi ^2} + \frac{\tilde{M}^4}{M_\varPhi ^4}\bigg ]\bigg \}\nonumber \\&+\frac{384 M_\varPhi ^6 \cot ^22\beta }{v_{\text {phys}}^6}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^4\bigg [\frac{\tilde{M}^2}{M_\varPhi ^2} + 8 \bigg ]\nonumber \\&+\frac{192 M_\varPhi ^4M_t^2 \cot ^2\beta }{v_{\text {phys}}^6}\left( 1-\frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^3 \bigg [\frac{\tilde{M}^2}{M_\varPhi ^2} + 8 \bigg ]\nonumber \\&-\frac{128 M_\varPhi ^6}{v_{\text {phys}}^6} \left( 1- \frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^5 \left( 2+\frac{\tilde{M}^2}{M_\varPhi ^2}\right) \nonumber \\&+ \frac{448 M_\varPhi ^4 M_t^2}{v_{\text {phys}}^6} \left( 1- \frac{\tilde{M}^2}{M_\varPhi ^2}\right) ^3 \left( 2+\frac{\tilde{M}^2}{M_\varPhi ^2} \right) \nonumber \\&+\mathscr {O}\left( \frac{M_\varPhi ^2M_t^4}{v_{\text {phys}}^6}\right) \,. \end{aligned}$$

(5.16)

In each of the two previous equations, the last two terms come from WF and VEV renormalisation. One can observe that all the terms in these expressions have the same form as Eq. (5.10) and therefore decouple for $M_\varPhi ^2=\tilde{M}^2+\tilde{\lambda } v^2$ and $\tilde{M}\rightarrow \infty$, as desired. Importantly, we should point out that, while we define our “on-shell” prescription for $\tilde{M}$ in terms of the calculation of the Higgs trilinear coupling, it also ensures the proper decoupling of BSM corrections to the Higgs quartic coupling, which provides a further validation of our results.

5.3 DM-inspired scenario of Inert-Doublet-model

We now turn to the dark-matter-inspired scenario of the IDM – discussed in Sect. 2.2 – where the CP-even inert state H constitutes a light DM candidate with mass $M_H\simeq M_h/2\ll M_A,\,M_{H^\pm }$, and where $\mu _2=0$. The dominant two-loop corrections to the Higgs self-couplings then come from the pseudoscalar and charged Higgs bosons, which because of their inert nature do not couple to fermions. The relevant diagrams in the two-loop effective potential are shown in Fig. 5, and their expressions read

$$\begin{aligned} V^{(2)}_{SSS}(h)=&-\frac{(v+h)^2}{8} \bigg [2 \lambda _A^2 I(m_A^2, m_A^2, 0) \nonumber \\&+ 4 \lambda _3^2 I(m_{H^\pm }^2, m_{H^\pm }^2, 0) \nonumber \\&+ 2(\lambda _3 - \lambda _A)^2 I(m_A^2, m_{H^\pm }^2,0)\nonumber \\&+ \lambda _A ^2 I(m_A^2, 0,0) +2\lambda _3^2 I(m_{H^\pm }^2,0,0)\bigg ]\,,\nonumber \\ V^{(2)}_{SS}(h)=&\ \frac{1}{8} \lambda _2\Big [3 A(m_A^2)^2 + 4 A(m_A^2) A(m_{H^\pm }^2) + 8 A(m_{H^\pm }^2)^2\Big ]\,, \end{aligned}$$

(5.17)

where all masses are understood to be field-dependent masses.

$\overline{{\mathrm{MS}}}$ expressions – Applying the operators $\mathscr {D}_3$ and $\mathscr {D}_4$, we can present here for the first time complete $\overline{{\mathrm{MS}}}$ expressions for the leading $\mathscr {O}(m_\varPhi ^6/v^5)$ and $\mathscr {O}(\lambda _2 m_\varPhi ^4/v^3)$ corrections – with $m_\varPhi$ being either $m_A$ or $m_{H^\pm }$ – to the Higgs trilinear and quartic couplings:

$$\begin{aligned} \delta ^{(2)}\lambda _{hhh}=&\ \frac{20 m_A^6 (-1 + 2{{\,\mathrm{\overline{\text {log}}}\,}}m_A^2)}{v^5}\nonumber \\&+\frac{40 m_{H^\pm }^6 (-1 + 2{{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2)}{v^5}\nonumber \\&+\frac{8(m_A^2-m_{H^\pm }^2)^2}{v^5}\big [-m_A^2-m_{H^\pm }^2\nonumber \\&+2m_A^2{{\,\mathrm{\overline{\text {log}}}\,}}m_A^2+2m_{H^\pm }^2{{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\big ]\nonumber \\&+\frac{6 \lambda _2 m_A^4 (1 + 2 {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2)}{v^3}\nonumber \\&+\frac{8 \lambda _2 m_A^2 m_{H^\pm }^2 (1 + {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2 + {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2)}{v^3}\nonumber \\&+\frac{16 \lambda _2 m_{H^\pm }^4 (1 + 2 {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2)}{v^3}\,, \end{aligned}$$

(5.18)

and

$$\begin{aligned} \delta ^{(2)}\lambda _{hhhh}=&\ \frac{160 m_A^6 {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2}{v^6}+\frac{320 m_{H^\pm }^6 {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2}{v^6}\nonumber \\&+\frac{64 (m_A^2 - m_{H^\pm }^2)^2}{v^6} \big [m_A^2 {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2 + m_{H^\pm }^2 {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\big ]\nonumber \\&+\frac{48 \lambda _2 m_A^4 (1 +{{\,\mathrm{\overline{\text {log}}}\,}}m_A^2) }{v^4}\nonumber \\&+\frac{32 \lambda _2m_A^2 m_{H^\pm }^2 (2 + {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2 + {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2 )}{v^4}\nonumber \\&+\frac{128 \lambda _2 m_{H^\pm }^4 (1+ {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2)}{v^4}\,. \end{aligned}$$

(5.19)

The results in both Eqs. (5.18) and (5.19) can be understood in terms of their correspondance to the effective-potential diagrams in Fig. 5. The first and second terms come from the derivatives of the leftmost two types of diagrams in Fig. 5 with respectively the pseudoscalar A or the charged $H^\pm$ running in the loops. The third term originates from the third sunrise diagram in Fig. 5, while the last three terms – proportional to $\lambda _2$ – arise from the eight-shaped diagrams.

Fig. 5

Diagrams contributing at leading two-loop order to the IDM effective potential

$\overline{{\mathrm{MS}}}$to OS scheme conversion – For the translation of these expressions to the OS scheme, we require the one-loop self-energies of A and $H^\pm$ in the IDM

$$\begin{aligned} \varPi ^{(1)}_{AA}(p^2)&=\ \frac{3}{2}\lambda _2A(m_A^2)+\lambda _2A(m_{H^\pm }^2)\nonumber \\&\quad -\frac{4m_A^4}{v^2} B_0(p^2,0,m_A^2) -\frac{m_A^4}{v^2}B_0(p^2,0,0)\nonumber \\&\quad -\frac{2(m_A^2-m_{H^\pm }^2)^2}{v^2}B_0(p^2,0,m_{H^\pm }^2)\,,\nonumber \\ \varPi ^{(1)}_{H^+H^-}(p^2)&=\ \frac{1}{2}\lambda _2A(m_A^2)+2\lambda _2A(m_{H^\pm }^2)\nonumber \\&\quad -\frac{4m_{H^\pm }^4}{v^2}B_0(p^2,0,m_{H^\pm }^2)\nonumber \\&\quad -\frac{m_{H^\pm }^4}{v^2}B_0(p^2,0,0)\nonumber \\&\quad -\frac{(m_A^2-m_{H^\pm }^2)^2}{v^2}B_0(p^2,0,m_A^2)\,, \end{aligned}$$

(5.20)

as well as the momentum-derivative of the one-loop self-energy of the 125-GeV Higgs boson

$$\begin{aligned} \frac{d}{dp^2}\varPi ^{(1)}_{hh}(p^2)\bigg |_{p^2=0}= & {} \frac{6M_t^2}{v_{\text {phys}}^2} \bigg ({{\,\mathrm{\overline{\text {log}}}\,}}M_t^2+\frac{2}{3}\bigg )\nonumber \\&-\frac{M_A^2}{3v_{\text {phys}}^2}-\frac{2 M_{H^\pm }^2}{3v_{\text {phys}}^2}\,. \end{aligned}$$

(5.21)

The conversion of the Higgs VEV is the same as in the SM, as given in Eq. (3.19). Moreover, as the coupling $\lambda _2$ only appears at two loops, we do not need any translation for it. Finally, had we not set $\mu _2$ to zero, it would have appeared in the one-loop correction to the Higgs self-couplings – c.f. Eqs. (3.10) and (3.11). However, for the conversion of $M^2$ in the 2HDM we found a shift proportional to $M^2$, and similarly here the tree-level ($\overline{{\mathrm{MS}}}$) relation $\mu _2=0$ will still hold after conversion to the OS scheme.

We then find in terms of on-shell scheme parameters the following expressions for the Higgs trilinear coupling

$$\begin{aligned}&\delta ^{(2)}\hat{\lambda }_{hhh}=\frac{6 \lambda _2}{v_{\text {phys}}^3}\big (3M_A^4+4M_A^2 M_{H^\pm }^2+8M_{H^\pm }^4\big )\nonumber \\&\quad +\frac{60 (M_A^6+2M_{H^\pm }^6)}{v_{\text {phys}}^5} +\frac{24(M_A^2-M_{H^\pm }^2)^2(M_A^2+M_{H^\pm }^2)}{v_{\text {phys}}^5}\nonumber \\&\quad +\frac{24M_t^4(M_A^2+2M_{H^\pm }^2)}{v_{\text {phys}}^5}+\frac{42M_t^2(M_A^4 +2M_{H^\pm }^4)}{v_{\text {phys}}^5}\nonumber \\&\quad -\frac{2(M_A^4+2M_{H^\pm }^4)(M_A^2+2M_{H^\pm }^2)}{v_{\text {phys}}^5}\,, \end{aligned}$$

(5.22)

and for the Higgs quartic coupling

$$\begin{aligned} \delta ^{(2)}\hat{\lambda }_{hhhh}=\frac{16}{3}\frac{\delta ^{(2)}\hat{\lambda }_{hhh}}{v_{\text {phys}}}\,. \end{aligned}$$

(5.23)

It is interesting to note that because of WF and VEV renormalisation – which give the last two lines of Eq. (5.22) – we can find terms involving both the inert-scalar and top-quark masses, even if these do not couple in the IDM. More specifically, these come from the interplay of the fermionic contributions to the Higgs WF and VEV renormalisation with the one-loop scalar contributions to the Higgs self-couplings, as well as of the scalar contributions to Higgs WFR with the one-loop top-quark corrections to the self-couplings.

5.4 A Higgs-Singlet model with $\mathbb {Z}_2$ symmetry

Finally, we turn to the $\mathbb {Z}_2$-symmetric HSM, introduced in Sect. 2.3 and in which we study the effects of the heavy additional real singlet S. This model is quite simple, and only two diagrams contribute to the two-loop effective potential at leading order, as shown in Fig. 6. The potential then is given by

$$\begin{aligned} V^{(2)}(h)= & {} -\frac{1}{4} \lambda _{HS}^2 (v + h)^2 I(m_S^2(h),m_S^2(h),0)\nonumber \\&+ \frac{3}{2} \lambda _S A(m_S^2(h))^2\,. \end{aligned}$$

(5.24)

Fig. 6

Dominant two-loop diagrams contributing to the HSM effective potential

$\overline{{\mathrm{MS}}}$ calculation – Following the same procedure as for the 2HDM and the IDM, we find in terms of $\overline{{\mathrm{MS}}}$ parameters

$$\begin{aligned} \delta ^{(2)}\lambda _{hhh}=&\ \frac{16m_S^4}{v^5} \left( 1 - \frac{\mu _S^2}{m_S^2}\right) ^4 \nonumber \\&\times \big [-m_S^2 - 2 \mu _S^2 + (2 m_S^2 + \mu _S^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_S^2\big ]\nonumber \\&+\frac{24m_S^4}{v^3} \left( 1 - \frac{\mu _S^2}{m_S^2}\right) ^3 \lambda _S \big [1 + 2 {{\,\mathrm{\overline{\text {log}}}\,}}m_S^2\big ]\,, \end{aligned}$$

(5.25)

for the corrections to the Higgs trilinear coupling, and

$$\begin{aligned} \delta ^{(2)}\lambda _{hhhh}=&\ \frac{32 m_S^2}{v^6} \left( 1-\frac{\mu _S^2}{m_S^2}\right) ^4 \big [-\mu _S^2 (4 m_S^2 + 5 \mu _S^2)\nonumber \\&+ (4 m_S^4 + 3 m_S^2 \mu _S^2 + 2 \mu _S^4) {{\,\mathrm{\overline{\text {log}}}\,}}m_S^2\big ]\nonumber \\&+\frac{48m_S^2}{v^4} \left( 1 - \frac{\mu _S^2}{m_S^2}\right) ^3 \lambda _S \big [4 m_S^2 - \mu _S^2\nonumber \\&+ 2 (2 m_S^2 + \mu _S^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_S^2\big ]\,, \end{aligned}$$

(5.26)

for those to the Higgs quartic coupling. In both equations, the first and second terms correspond respectively to the left and right effective-potential diagrams in Fig. 6. As a (partial) cross-check of our calculation, we have confirmed that if we take the fourth derivative $\partial ^4/\partial h^4$ of $V^{(2)}$ (instead of applying $\mathscr {D}_4$), we reproduce the same result as in Eq. (29) of Ref. [35] for the two-loop corrections to the Higgs quartic coupling (after taking the limit of $m_h\rightarrow 0$ in said equation) – note that the results in Ref. [35] were obtained in a diagrammatic computation (at zero external momentum), using results generated by the Mathematica package SARAH [73,74,75,76,77,78].

Conversion to the OS scheme – To convert the above expressions to the on-shell scheme, we need here first the one-loop self-energy of S, which reads

$$\begin{aligned} \varPi ^{(1)}_{SS}(p^2)=6\lambda _SA(m_S^2)-\frac{4(m_S^2-\mu _S^2)^2}{v^2} B_0(p^2,0,m_S^2)\,, \end{aligned}$$

(5.27)

after neglecting the light scalar masses.

Moreover, as in the 2HDM, the mass parameter $\mu _S$ appears in the one-loop corrections to the Higgs couplings. Therefore, we also need a finite “OS” counterterm – that we denote $\delta ^\text {OS}\mu _S^2$ – for it, which we define as

$$\begin{aligned} \tilde{\mu }_S^2=\mu _S^2+\delta ^\text {OS}\mu _S^2\,, \end{aligned}$$

(5.28)

where $\tilde{\mu }_S$ and $\mu _S$ are the OS- and $\overline{{\mathrm{MS}}}$-renormalised versions of the mass parameter. As in the 2HDM, we determine $\delta ^\text {OS}\mu _S^2$ by requiring that it ensures the proper decoupling of the two-loop corrections to the Higgs trilinear when using a relation of the form $M_S^2=\tilde{\mu }_S^2+\frac{1}{2}\lambda _S v^2$ and taking the limit $\tilde{\mu }_S\rightarrow \infty$ while keeping $\lambda _S$ fixed. We find eventually

$$\begin{aligned} \delta ^\text {OS}\mu _S^2=6\kappa \lambda _S\tilde{\mu }_S^2({{\,\mathrm{\overline{\text {log}}}\,}}M_S^2-1)\,. \end{aligned}$$

(5.29)

We use the corrections to the Higgs trilinear coupling to determine the finite counterterm $\delta ^\text {OS}\mu _S^2$, but verifying that it also fulfills the above requirement for the quartic coupling provides an important cross-check of our result.

Finally, we require the derivative of the one-loop 125-GeV Higgs-boson self-energy

$$\begin{aligned} \frac{d}{dp^2}\varPi ^{(1)}_{hh}(p^2)\bigg |_{p^2=0}= & {} \frac{6M_t^2}{v_{\text {phys}}^2}\bigg ({{\,\mathrm{\overline{\text {log}}}\,}}M_t^2+\frac{2}{3}\bigg )\nonumber \\&-\frac{M_S^2}{3v_{\text {phys}}^2}\left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^2\,. \end{aligned}$$

(5.30)

Using all the above results, we obtain finally the following OS-renormalised expressions for the two-loop corrections to the Higgs self-couplings

$$\begin{aligned} \delta ^{(2)}\hat{\lambda }_{hhh}=&\ \frac{48 M_S^6}{v_{\text {phys}}^5}\left( 1 - \frac{\tilde{\mu }_S^2}{M_S^2}\right) ^4+\frac{72\lambda _S M_S^4}{v_{\text {phys}}^3}\left( 1 - \frac{\tilde{\mu }_S^2}{M_S^2}\right) ^3 \nonumber \\&+\frac{24 M_S^2 M_t^4}{v_{\text {phys}}^5}\left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^2 \nonumber \\&+\frac{42 M_S^4 M_t^2}{v_{\text {phys}}^5} \left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^3-\frac{2M_S^6}{v_{\text {phys}}^5} \left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^5\,,\nonumber \\ \delta ^{(2)}\hat{\lambda }_{hhhh}=&\ \frac{32 M_S^6}{v_{\text {phys}}^6}\left( 1 - \frac{\tilde{\mu }_S^2}{M_S^2}\right) ^4\bigg [8 + \frac{2 \tilde{\mu }_S^2}{M_S^2} - \frac{\tilde{\mu }_S^4}{M_S^4}\bigg ]\nonumber \\&+\frac{48 \lambda _S M_S^4}{v_{\text {phys}}^4}\left( 1 - \frac{\tilde{\mu }_S^2}{M_S^2}\right) ^3\bigg [8 + \frac{\tilde{\mu }_S^2}{M_S^2}\bigg ]\nonumber \\&+\frac{128 M_S^2M_t^4}{v_{\text {phys}}^6}\left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^2 \nonumber \\&+\frac{112 M_S^4 M_t^2}{v_{\text {phys}}^6} \left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^3\left( 2 + \frac{\tilde{\mu }_S^2}{M_S^2}\right) \nonumber \\&-\frac{16 M_S^6}{3v_{\text {phys}}^6}\left( 1-\frac{\tilde{\mu }_S^2}{M_S^2}\right) ^5\left( 2 + \frac{\tilde{\mu }_S^2}{M_S^2}\right) \,. \end{aligned}$$

(5.31)

As we had observed already in the IDM, even if the additional scalar S does not couple directly to the top quark, the finite Higgs WF and VEV renormalisations introduce terms that involve both $M_S$ and $M_t$.

Fig. 7

Illustrations of the decoupling of the deviations of ${\hat{\lambda }}_{hhh}$ calculated in the 2HDM with respect to its SM prediction. (Left side): BSM deviations $\delta R^{1\ell }$ and $\delta R^{2\ell }$ – defined in Eq. (6.2) – respectively at one loop (solid blue curve) and at two loops (dot-dashed red curve) as a function of $\tilde{M}$. Results are shown for $\tan \beta =1.5$ and for several values of the difference $M_\varPhi ^2-\tilde{M}^2={\tilde{\lambda }}v^2$, namely $(200\text { GeV})^2$, $(300\text { GeV})^2$, and $(400\text { GeV})^2$, where $M_\varPhi$ is the degenerate mass of the heavy 2HDM scalars. (Right side): Behaviour of the two-loop BSM contributions to ${\hat{\lambda }}_{hhh}$ in the 2HDM, shown here with $\delta R^{2\ell }-\delta R^{1\ell }$, as a function of $\tilde{M}$ and for several values of $\tan \beta$, up to the higher value allowed under the criterion of tree-level perturbative unitarity. For this figure, the masses of the additional BSM scalar are once again taken to be degenerate, and $M_\varPhi ^2-\tilde{M}^2=(400\text { GeV})^2$

6 Numerical examples

We now turn to the discussion of the numerical behaviour of the BSM corrections computed in Sect. 5. Before looking at concrete examples, a comment should be made about the theoretical and experimental constraints that we include in our analysis. On the theory side, in addition of the potential being bounded from below (as discussed in Sect. 2), we require that unitarity should not be violated. For this, we choose to take as our criterion⁹ that tree-level perturbative unitarity [80] should hold, and we employ for the 2HDM and the IDM the results of Refs. [81, 82] and for the HSM those of Ref. [35, 127]. On the experimental side, we here use the public program HiggsBounds-5.3.2beta [128,129,130,131] to take into account constraints from searches at LEP, the Tevatron, and the LHC on the allowed parameter spaces of the BSM scenarios we investigate. To obtain the input files for HiggsBounds in the different models that we study, specific spectrum generators based on SPheno [71, 72] are created using SARAH.

6.1 Decoupling limit

A natural first point to study numerically is the decoupling behaviour of the two-loop BSM corrections. In Sect. 5, we had discussed the decoupling of BSM effects in terms of analytical expressions, finding that the corrections¹⁰ to Higgs self-couplings in the $\overline{{\mathrm{MS}}}$ scheme have the form

$$\begin{aligned} (m_\varPhi ^2)^m\left( 1-\frac{\mathscr {M}^2}{m_\varPhi ^2}\right) ^n,\text { with }m<n\,, \end{aligned}$$

(6.1)

where $\mathscr {M}$ is defined in Eq. (3.12), thereby ensuring that the radiative corrections indeed vanish when decoupling additional states in the extended Higgs sectors. We have also devised schemes for the BSM mass parameters $\mathscr {M}$ – i.e. M in the 2HDM and $\mu _S$ in the HSM – so that OS-scheme expressions also have a similar form in terms of the OS-renormalised parameters $M_\varPhi$ and either $\tilde{M}$ or $\tilde{\mu }_S$. As the expressions for the 2HDM and the HSM are very similar, we will for concreteness consider the case of the 2HDM in the following.

The left part of Fig. 7 illustrates the decoupling of the one- and two-loop corrections to the Higgs trilinear coupling ${\hat{\lambda }}_{hhh}$ in the 2HDM, when expressed in terms of OS parameters. More precisely, the plot shows, as a function of $\tilde{M}$, the BSM deviations $\delta R$, defined – in the alignment limit – at one and two loops respectively as

$$\begin{aligned} \delta R^{1\ell }&=\frac{({\hat{\lambda }}_{hhh}^\text {2HDM})^{(1)}}{({\hat{\lambda }}_{hhh}^\text {SM})^{(1)}}-1=\frac{\kappa (\delta ^{(1)}{\hat{\lambda }}_{hhh}^\text {2HDM}-\delta ^{(1)}{\hat{\lambda }}_{hhh}^\text {SM})}{{\hat{\lambda }}_{hhh}^\text {tree}+\kappa \delta ^{(1)}{\hat{\lambda }}_{hhh}^\text {SM}}\,,\nonumber \\ \delta R^{2\ell }&=\frac{({\hat{\lambda }}_{hhh}^\text {2HDM})^{(2)}}{({\hat{\lambda }}_{hhh}^\text {SM})^{(2)}}-1\nonumber \\&=\frac{\kappa (\delta ^{(1)}{\hat{\lambda }}_{hhh}^\text {2HDM}-\delta ^{(1)}{\hat{\lambda }}_{hhh}^\text {SM})+\kappa ^2(\delta ^{(2)}{\hat{\lambda }}_{hhh}^\text {2HDM}-\delta ^{(2)}{\hat{\lambda }}_{hhh}^\text {SM})}{{\hat{\lambda }}_{hhh}^\text {tree}+\kappa \delta ^{(1)}{\hat{\lambda }}_{hhh}^\text {SM}+\kappa ^2 \delta ^{(2)}{\hat{\lambda }}_{hhh}^\text {SM}} \,, \end{aligned}$$

(6.2)

for an example point with degenerate BSM scalar masses and $\tan \beta =1.5$. For each of the values of the difference $M_\varPhi ^2-\tilde{M}^2={\tilde{\lambda }}v^2$ that we consider – namely $(200\text { GeV})^2$, $(300\text { GeV})^2$, and $(400\text { GeV})^2$ – we can indeed observe that the BSM effects decouple rapidly for increasing $\tilde{M}$, both at the one-loop level (blue solid curves) and at the two-loop level (red dot-dashed curves).

Another interesting point that we can study is the new $\tan \beta$ dependance at two loops and its impact on the decoupling of BSM corrections. Indeed, as can be seen from Eqs. (5.15) and (5.16), the two-loop contributions to the Higgs self-couplings involve $\cot ^2\beta$ and $\cot ^22\beta$, even in the alignment limit, because these appear in tree-level scalar couplings. While $\cot ^2\beta$ obviously vanishes very quickly with increasing $\tan \beta$, $\cot ^22\beta$ grows very fast – like $\tan ^2\beta$ for large $\tan \beta$. This implies possible enhancements of the two-loop corrections for large values of $\tan \beta$, only limited by the upper bound that the constraint of perturbative unitarity puts on $\tan \beta$. We show on the right side of Fig. 7 the magnitude of the two-loop deviations – obtained as $\delta R^{2\ell }-\delta R^{1\ell }$ – as a function of $\tilde{M}$. We present results for four values of $\tan \beta$, $\tan \beta =1.75$ being close to the maximal value allowed under the criterion of tree-level perturbative unitarity [81] for $\tilde{M}=0$ and $M_\varPhi =400\text { GeV}$. As expected, while $\delta R^{1\ell }$ does not depend on $\tan \beta$, we can observe significant enhancements of $\delta R^{2\ell }$ when increasing $\tan \beta$, especially for small values of $\tilde{M}$. When $\tilde{M}$ grows, the effect of the terms in $\delta ^{(2)}{\hat{\lambda }}_{hhh}$ proportional to $\cot ^22\beta$ diminishes because these have higher powers of the suppression factor than some of the other, non-$\tan \beta$-enhanced terms. All in all, even at the limit of the region of parameter space allowed by unitarity, the decoupling of the BSM corrections from additional scalars occurs rapidly – not significantly slower than for smaller $\tan \beta$.

Fig. 8

Illustration of the non-decoupling behaviour of the BSM corrections to the Higgs trilinear coupling ${\hat{\lambda }}_{hhh}$, at one and two loop(s). (Upper left): Comparison of the results for $\delta R$ in the 2HDM (solid), the IDM (dot-dashed), and the HSM (dashed) at both one-loop (blue) and two-loop (red) orders, as a function of the degenerate mass $M_\varPhi$ of all BSM scalars. We also choose: $\tan \beta =1.1$ for the 2HDM; $\lambda _2=0.5$ for the IDM; and $\lambda _S=0.5$ for the HSM. (Upper right): Non-decoupling behaviour of the BSM effects in the 2HDM as a function of $M_\varPhi$, at one loop (blue) and two loops (red). At two-loop order, the two curves correspond to $\tan \beta =1$ (dashed) and $\tan \beta =1.4$ (dot-dashed), the latter being the maximal value of $\tan \beta$ allowed under the criterion of tree-level unitarity [80, 81] for $M_\varPhi =500\text { GeV}$. (Lower left): Similar plot as upper right one, for the case of the IDM. For the two-loop curves, $\lambda _2=0$ (dashed) and $\lambda _2=6$ (dot-dashed) correspond to the extreme values of $\lambda _2$ allowed respectively from stability of the potential and tree-level unitarity. (Lower right): Similar plot as upper right one, for the case of the HSM – as there is only one scalar S, we denote its mass $M_S$ instead of $M_\varPhi$ here. Again, $\lambda _S=3.7$ is the largest possible value to fulfill the criterion of tree-level unitarity [35, 80] until $M_S=500\text { GeV}$

6.2 Non-decoupling effects

After having verified the proper decoupling behaviour of our results, we can now turn to the more important question of the non-decoupling effects and of the maximal possible size they can reach. Due to the presence of the reduction factors, as shown in Eq. (6.1), in the corrections to the Higgs self-couplings, the largest BSM effects will be found for smaller values of $\mathscr {M}$ – i.e. when the BSM scalars obtain their masses mostly from the Higgs VEV and quartic couplings, and can therefore not be decoupled.

Figure 8 illustrates our results for the BSM deviations $\delta R$ of the Higgs trilinear coupling ${\hat{\lambda }}_{hhh}$ in this limit. In the upper left plot, we compare the magnitude of the deviations obtained in the three different models considered in this paper, at one loop (blue curves) and at two loops (red curves), for scenarios where the additional scalar are degenerate in mass (and with the remaining parameters fixed as indicated in the caption of Fig. 8). As could be expected from the analytical expressions found in Sect. 5, the one- and two-loop corrections have very similar behaviours in the three models – the two-loop effects giving additional positive contributions to the BSM deviation. However, one can immediately notice the numerical discrepancies between theories, arising mainly from the different number of BSM degrees of freedom: one for the HSM, three for the IDM, and four for the 2HDM (we recall that the charged Higgs in the 2HDM and IDM is associated with two degrees of freedom). In all cases, we can observe that the two-loop corrections grow faster than the one-loop one, due to the $M_\varPhi ^6$ dependence of part of the two-loop terms – see Eqs. (5.15), (5.22), and (5.31) – but importantly they remain well below the size of the one-loop effects for the entire mass range considered, for which we have verified that perturbative unitarity is not violated.

In the other three subplots of Fig. 8 (upper right and lower ones), we show the behaviour of the BSM corrections as a function of the BSM scalar mass separetely for the three models. In each case, the one-loop results are presented in blue, while the possible range of two-loop values is given between the red curves. Each model includes an additional parameter – $\tan \beta$ for the 2HDM, $\lambda _2$ for the IDM, and $\lambda _S$ for the HSM – that enters the expressions for the Higgs self-couplings only at two loops. We choose here to vary these parameters from the smallest possible value they can take, for which the eight-shaped diagrams in the effective potential vanish, respectively $\tan \beta =1$, $\lambda _2=0$, and $\lambda _S=0$, up to the maximal value for which the condition of tree-level unitarity is verified all the way to $M_\varPhi =500\text { GeV}$ – this gives the maximal values $\tan \beta =1.4$, $\lambda _2=6.0$, and $\lambda _S=3.7$. Turning first to the 2HDM, we observe that varying $\tan \beta$ in the allowed range does not produce any significant enhancement of the two-loop corrections; in other words the eight-shaped diagrams in $V_{{\text {eff}}}$ only have a limited impact in the 2HDM. Indeed requiring that tree-level unitarity is preserved up to $M_\varPhi =500\text { GeV}$ and for $\tilde{M}=0$ strongly constrains the maximal possible value of $\tan \beta$ – we will return to this point in the next section. In contrast to the 2HDM, in the IDM and the HSM the parameters $\lambda _2$ and $\lambda _S$ – the quartic couplings of the inert doublet and of the inert singlet respectively – are less severely constrained by unitarity. For this reason, they can cause large two-loop contributions, as can be seen from the orange-shaded areas in the lower subplots of Fig. 8.

Regarding the numerical impact of the two-loop corrections, if for example we consider the same parameter choices as for Fig. 8, we can note first that the two-loop corrections are minute for small masses, say $M_\varPhi \lesssim v$; however this is also where the approximations made in our calculation are least justified. Considering next the points for $M_\varPhi =400\text { GeV}$ – i.e. well within the region where perturbative unitarity conditions are fulfilled – we find that the two-loop corrections can at most become as large as 25%, 52%, and 62% of the one-loop corrections (respectively for the 2HDM, the IDM, and the HSM). If we turn off the eight-shaped diagrams and consider only the effect of the sunrise diagrams in $V_{{\text {eff}}}$, these shifts are reduced to 22%, 24%, and 21% respectively. Even if we push our computation to the limit of the parameter region allowed by perturbative unitary – here at $M_\varPhi =500\text { GeV}$ for the maximal values of $\tan \beta =1.4$, $\lambda _2=6$, and $\lambda _S=3.7$ – the two-loop corrections remain smaller than the one-loop ones, with respectively 38%, 64%, and 73% of the one-loop results. In conclusion, on the one hand, in the non-decoupling limit the two-loop corrections are not always negligible before the one-loop corrections and in the case of the IDM and HSM new enhancements from the quartic couplings $\lambda _2$ and $\lambda _S$ can appear. On the other hand, however, as long as perturbative unitarity is not violated, the two-loop corrections to $\hat{\lambda }_{hhh}$ are smaller than their one-loop counterparts, and the validity of the perturbation expansion is not in doubt.

As a final remark, we can comment also about the experimental limits on the parameter regions considered in Fig. 8. First, for the 2HDM, comprehensive studies of experimental constraints can be found e.g. in Refs. [99, 132]. In addition to searches of charged and neutral scalars at the LHC (which are included in HiggsBounds), results from flavour Physics – in particular decays like $b\rightarrow s\gamma$ – can also limit severely the allowed values of $M_{H^\pm }$ and $\tan \beta$, and for this reason we choose to work in a 2HDM of type I (as mentioned already in Sect. 2) where these constraints are weakest. Indeed in type I, the lower limit on $M_{H^\pm }$ decreases from approximately $440\text { GeV}$ for $\tan \beta =1$ to below $200\text { GeV}$ for $\tan \beta =1.5$ [99]. Additionally, with HiggsBounds-5.3.2beta, we have verified that for $\tilde{M}=0$ and $\tan \beta \in [1,1.4]$, BSM scalar masses above 355 GeV are still allowed. For the IDM, as we consider a scenario where H is a DM candidate of mass $M_H\simeq M_h/2$, collider and DM searches do not constrain the mass range $M_\varPhi =M_A=M_{H^\pm }$ between 200 and 500 GeV [133]. Finally, for the HSM, we have considered a scenario with an additional $\mathbb {Z}_2$ symmetry under which S is charged, and is thus inert. This considerably weakens the existing constraints (that can become strong if h and S are able to mix), and the range of masses $M_S$ between 200 and 500 GeV is also not constrained here [134].

6.3 Maximal possible size of the BSM corrections

Another question that we can consider is by how much the two-loop result for the Higgs trilinear in an extended sector $\hat{\lambda }_{hhh}^\text {BSM}$ can deviate from its SM prediction $\hat{\lambda }_{hhh}^\text {SM}$, given the constraint of perturbative unitarity, in particular if we consider a broader range of BSM masses and parameters. For concreteness, we concentrate here on the case of the aligned scenario of the 2HDM with mass-degenerate additional scalars. Figure 9 shows the maximal possible size of the BSM deviation at two loops $\delta R^{2\ell }$ in the plane of $M_\varPhi$ and $\tan \beta$, considering now significantly larger ranges than in the previous section. Once these two parameters are fixed, the value of $\tilde{M}$ alone determines how large the corrections can be: the criterion of tree-level perturbative unitarity essentially yields upper bounds on the Lagrangian scalar quartic coupling, which via the relations of the form $M_\varPhi ^2=\tilde{M}^2+{\tilde{\lambda }}v^2$ translate into a lower bound on $\tilde{M}$. For $M_\varPhi$ sufficiently small, as considered in Fig. 8 for instance, $\tilde{M}$ can be taken to zero without violating the unitarity conditions. However, if we now increase $M_\varPhi$ we reach a point when the lower bound on $\tilde{M}$ is non-vanishing, and then the reduction factors shown in Eq. (6.1) enter into play and diminish the size of the BSM effects. Additionally, the value of $M_\varPhi$ until which $\tilde{M}$ can remain zero diminishes with larger $\tan \beta$ as the unitarity relations then become increasingly stringent. In turn, this limits the maximal size the BSM corrections can reach for larger $\tan \beta$.

This is indeed what we can observe if we consider a horizontal (i.e. constant $\tan \beta$) section of Fig. 9: first when $M_\varPhi$ increases, $\delta R^{2\ell }$ grows rapidly and deviations of more than 400% from the SM prediction are possible for small values of $\tan \beta$. Then the point where $\tilde{M}$ cannot remain vanishing is reached and the BSM corrections decrease in size for even larger $M_\varPhi$. Therefore, large deviations in the Higgs trilinear coupling are only possible in a relatively limited area of parameter space, for low values of $\tan \beta$ and intermediate BSM scalar masses. Given the levels of precision envisioned for the measurement of the Higgs trilinear coupling at future colliders (as discussed in the introduction), the blue-shaded region in Fig. 9 corresponds (roughly) to the part of parameter space that could be probed at the HL-LHC, while the green-shaded one illustrates the potential reach of high-energy lepton colliders, such as the 1-TeV ILC, or the 3-TeV CLIC. This discussion also shows that investigating the Higgs trilinear coupling can allow probing the low $\tan \beta$ and intermediate $M_\varPhi$ region of the 2HDM parameter space, in complementarity to $H^{\pm }\rightarrow t b$ searches at high-energy colliders.

Fig. 9

Contour plot of the maximal possible BSM deviation of ${\hat{\lambda }}_{hhh}$ evaluated at two loops in the 2HDM from the SM result when requiring tree-level perturbative unitarity to hold, in the plane of $M_\varPhi$ and $\tan \beta$. Here again, $M_\varPhi$ is the common mass of the heavy BSM scalars in the 2HDM, and the figure is made in the alignment limit (i.e. $s_{\beta -\alpha }=1$)

Fig. 10

Comparison of the BSM corrections $\varDelta \lambda _{hhh}^\text {IDM}\equiv \lambda _{hhh}^\text {IDM}-\lambda _{hhh}^\text {SM}$ calculated at one- and two-loop orders (blue and red curves respectively) in the $\overline{{\mathrm{MS}}}$ scheme (dashed curves) and in the OS scheme (solid curves). To obtain the $\overline{{\mathrm{MS}}}$ results, we convert the OS-renormalised input parameters ($M_A,\,M_{H^\pm },\,M_t,\,v_{\text {phys}}$) to $\overline{{\mathrm{MS}}}$, using the expressions presented in Secs. 3 and 5. (Left side): $\varDelta \lambda _{hhh}^\text {IDM}$ as a function of the pole mass $M_\varPhi$, for parameters points where $M_A=M_{H^\pm }=M_\varPhi$ and for $\lambda _2=0.5$. For the $\overline{{\mathrm{MS}}}$ results, we fix $Q=M_\varPhi$ here. (Right side): $\varDelta \lambda _{hhh}^\text {IDM}$ as a function of the renormalisation scale Q, for the parameter point defined by $M_A=M_{H^\pm }=M_\varPhi =400\text { GeV}$ and $\lambda _2=0.5$

6.4 An estimate of the theoretical uncertainty

Having discussed the possible size of the dominant two-loop BSM corrections to $\hat{\lambda }_{hhh}$ and $\hat{\lambda }_{hhhh}$, it is natural to conclude by a discussion of the theoretical uncertainty associated with our results. First of all, we need to clarify the type of effects that we want to estimate here, and for concreteness, we consider the case of the IDM, studied in Sect. 5.3. The types of corrections that we have included in our calculations are the $\mathscr {O}(M_\varPhi ^4/v_{\text {phys}}^3)$ effects at one loop, and the $\mathscr {O}(M_\varPhi ^6/v_{\text {phys}}^5)$ and $\mathscr {O}(\lambda _2M_\varPhi ^4/v_{\text {phys}}^3)$ ones at two loops – where by $M_\varPhi$ we mean here either $M_A$ or $M_{H^\pm }$, or some combination of both. Therefore, we choose to gauge the size of the leading three-loop corrections – of the form $\mathscr {O}(M_\varPhi ^8/v_{\text {phys}}^7)$, $\mathscr {O}(\lambda _2M_\varPhi ^6/v_{\text {phys}}^5)$, and $\mathscr {O}(\lambda _2^2M_\varPhi ^4/v_{\text {phys}}^3)$ – by comparing the results obtained in the $\overline{{\mathrm{MS}}}$ and OS schemes for the BSM corrections to the Higgs trilinear coupling.

For this comparison, we take the OS-renormalised parameters $M_A,\,M_{H^\pm },\,M_t,\,v_{\text {phys}}$ as our inputs: we use them directly for the OS scheme results, and convert them into $\overline{{\mathrm{MS}}}$-renormalised parameters – using Eqs. (3.19), (5.4), and (5.20) – for use in our $\overline{{\mathrm{MS}}}$ expressions. Figure 10 illustrates the results we obtain for the BSM corrections $\varDelta \lambda _{hhh}^\text {IDM}$ calculated either in the OS scheme, i.e. $\hat{\lambda }_{hhh}^\text {IDM}-\hat{\lambda }_{hhh}^\text {SM}$, or in the $\overline{{\mathrm{MS}}}$ scheme, i.e. $\lambda _{hhh}^\text {IDM}-\lambda _{hhh}^\text {SM}$, at both one and two loop(s). The left-hand plot shows the comparison of these different results as a function of the degenerate mass $M_\varPhi$ of the heavy inert scalars A and $H^\pm$. While the difference between the one-loop curves (in blue), which measures the typical size of two-loop corrections, is very large – up to 75% for $M_\varPhi =500\text { GeV}$, the two-loop results (red curves) are much closer, and differ by at most 13% (again for $M_\varPhi =500\text { GeV}$). If we consider for example the situation for $M_\varPhi =400\text { GeV}$, we find $m_A(Q=M_\varPhi )=m_{H^\pm }(Q=M_\varPhi )=434\text { GeV}$, $m_t(Q=M_\varPhi )=157\text { GeV}$, and $v(Q=M_\varPhi )=241\text { GeV}$. In turn, with these values, we obtain for $\varDelta \lambda _{hhh}^\text {IDM}$ at one loop 192 GeV ($\overline{{\mathrm{MS}}}$) and 130 GeV (OS) respectively, and at two loops 180 GeV ($\overline{{\mathrm{MS}}}$) and 168 GeV (OS). This indicates a significant decrease in the theoretical uncertainty on the prediction of $\hat{\lambda }_{hhh}$. It should be noted that the apparant small size of the two-loop corrections in the $\overline{{\mathrm{MS}}}$ scheme comes in part from a cancellation between the terms involving $\lambda _2$ and those independent of it, which come with opposite signs (for $Q=m_\varPhi$).

The right-hand plot of Fig. 10 shows also the renormalisation scale dependence of the $\overline{{\mathrm{MS}}}$ results for $\varDelta \lambda _{hhh}^\text {IDM}$ at one- and two-loop orders, compared to the OS values. A first positive point to notice is the reduced renormalisation scale dependence of the two-loop result (where we do not truncate the scale dependence arising at one loop) compared to the one-loop one. Furthermore, for most of the range of scales shown in the figure, the agreement between schemes is significantly better at two loops than at one loop – note , however, that for large values of Q the $\overline{{\mathrm{MS}}}$ calculation breaks down because of unphysically large logarithmic terms of the form $\log m_\varPhi ^2/Q^2$. Varying the renormalisation scale in the $\overline{{\mathrm{MS}}}$ calculation by factors of 1/2 and 2 around the natural value $Q=400\text { GeV}$ yields changes of up to $4.6\%$, which gauge the size of three-loop (subleading) logarithmic terms of the form $\mathscr {O}(M_\varPhi ^8/v_{\text {phys}}^7)$, $\mathscr {O}(\lambda _2M_\varPhi ^6/v_{\text {phys}}^5)$, and $\mathscr {O}(\lambda _2^2M_\varPhi ^4/v_{\text {phys}}^3)$. In conclusion, we keep as our estimate of the theoretical uncertainty on our two-loop result the value of approximately 5% (of the total result).

7 Discussion

Beyond their inherent calculational aspects, our new results for the radiative corrections to Higgs self-couplings, and in particular the Higgs trilinear coupling, have important consequences for the physics of the considered BSM scenarios. First and foremost, we should emphasize that these are important examples of scenarios where some new physics does exist close to the electroweak scale, but is hidden from observation either by alignment or by some global $\mathbb {Z}_2$ symmetry. It is then only via the precise study of Higgs boson properties – and possible non-decoupling effects – that these scenarios can be distinguished from situations where new states are heavy and thus decoupled. Our results demonstrate that these non-decoupling effects – found first at one-loop order – are by no means calculational artefacts caused by a breakdown of perturbation theory, but true physical effects derived in a properly converging perturbative expansion. Furthermore, we have shown that the new corrections at two loops, while always well below the size of one-loop effects and typically mild, give positive enhancements of the Higgs self-coupling. It should be noted however that, given the prospects for the measurement of $\lambda _{hhh}$ at future colliders these two-loop contributions could potentially be distinguishable experimentally in the future. For instance, supposing a deviation of $\mathscr {O}(100\%)$ in the Higgs trilinear coupling (from its SM prediction) were to be found at the HL-LHC, the accuracy obtained at lepton colliders – down to $10\%$, c.f. the discussion in the introduction – would require theoretical predictions at two loops to properly interpret the measurements in terms of the parameter space of BSM models. One may expect that this observation would hold also for other Higgs couplings (e.g. to gauge bosons or fermions) because, even if the non-decoupling effects these couplings exhibit are smaller, the predicted accuracy of their future measurements is significantly better than for $\lambda _{hhh}$ – in turn, this should also motivate us to study higher-order BSM corrections to these other couplings of the Higgs boson.

Returning to the case of the Higgs trilinear coupling, a natural avenue for further work would be to continue the calculation of new two-loop corrections, with both effective-potential and diagrammatic methods. In particular, this would imply investigating the impact of non-zero external momenta at two loops – with extensive work necessary in this direction – as well as the size of subleading two-loop effects. However, for the latter, we may expect that when we include Goldstone bosons in the $\overline{{\mathrm{MS}}}$ calculations, we will encounter IR divergences in the limit $m_{G,G^\pm }\rightarrow 0$ even if we include the dependence on external momenta (this is the so-called “Goldstone Boson Catastrophe” [119]), and it will be unavoidable to employ an OS renormalisation scheme for the Goldstone boson masses as was done in Ref. [120] (for other solutions see also [118, 135,136,137,138,139]).

The precise calculation of the Higgs trilinear coupling is also of great importance to relate the properties of theories with extended Higgs sectors to BSM phenomena in these models, such as for example the possibility of a strong first-order electroweak phase transition. In particular, our findings in this work motivate considering the calculation of two-loop contributions to $\lambda _{hhh}$ at finite temperature in models with extended sectors, in order to study their effect on the strength of the EWPT. This was considered already for the case of the IDM in Ref. [85], and a mild weakening of the EWPT was found. Moreover, the complementarity and synergy between the measurement of $\lambda _{hhh}$ at future colliders and the future measurement of gravitational waves from a strong first-order EWPT at LISA and DECIGO will explore the shape of the Higgs potential, and further clarify the physics behind EWSB – see for instance Refs. [140,141,142,143].

Turning finally to the IDM and the HSM (with an exact $\mathbb {Z}_2$ symmetry), these are two examples of models with inert sectors that can host a candidate of DM particle while evading current experimental searches. In both models, the inert scalars have a quartic self-coupling – respectively $\lambda _2$ for the IDM and $\lambda _S$ for the HSM – that is difficult to probe experimentally but plays an important role for the physics of the inert sector. Interestingly, these couplings appear in the radiative corrections to Higgs self-couplings, starting from the two-loop order. Moreover, they will also appear in other couplings of the Higgs boson, such as the $h\gamma \gamma$ coupling, which is already measured to a high accuracy and can be accessed to percent level at future lepton colliders such as the ILC (see e.g. Ref. [46]). Calculating higher-order corrections to Higgs couplings therefore offers a new way to probe the hidden dynamics of theories with inert sectors.

8 Conclusion

In this paper, we have investigated two-loop corrections to the Higgs self-couplings, building on our work in Ref. [86] and giving details about our methods and calculations.

We have presented new general results, in terms of $\overline{{\mathrm{MS}}}$-renormalised parameters, for the derivatives of integrals appearing in the effective potential, which can be applied to further models (in the absence of scalar mixing) and also served for important cross-checks of our model-specific computations. Indeed, we have also calculated the dominant two-loop corrections to the Higgs trilinear and quartic couplings in three particular BSM theories, namely a Two-Higgs-Doublet Model in the alignment limit, the Inert-Doublet Model, and the Higgs Singlet Model (with an exact $\mathbb {Z}_2$ symmetry). We have provided expressions for these corrections both in the $\overline{{\mathrm{MS}}}$ scheme and in the on-shell scheme. In particular, we have explained our modified “on-shell” prescription (originally presented in Ref. [86]) for the renormalisation of the soft-breaking scale $\tilde{M}$ of the $\mathbb {Z}_2$ symmetry in the 2HDM, ensuring the explicit decoupling of the BSM corrections, expressed in terms of OS-renormalised parameters, when taking the limit $\tilde{M}\rightarrow \infty$ and using a relation of the form $M_\varPhi ^2=\tilde{M}^2+{\tilde{\lambda }}v^2$. We have also extended this result to the case of $\mu _S$ in the HSM, and interestingly we found (for both the 2HDM and the HSM) that while our prescription is defined in terms of the corrections to the Higgs trilinear coupling, it directly works for the corrections to the Higgs quartic coupling.

Furthermore, we have performed an extended numerical analysis of our results for the Higgs trilinear coupling, confirming that our prescription to renormalise the additional mass scale ($\tilde{M}$ or $\tilde{\mu }_S$) works properly, and examining the behaviour of the two-loop corrections. We have shown that the leading two-loop corrections, when expressed in terms of OS-renormalised parameters, give a positive enhancement of the one-loop results. We have moreover investigated the maximal size that these corrections can reach, as well as the possibility of large new effects at two loops. Indeed, in all three studied models, a new parameter appears in the radiative corrections at two loops and this raises the question of whether new large corrections are possible. We find that, with the requirement of tree-level unitarity, there is no large effect in the 2HDM, while in the IDM and HSM the corrections originating from the eight-shaped diagrams in the effective potential can become the leading two-loop contributions. However, in all three cases, the two-loop corrections do remain smaller than their one-loop counterparts, at least as long as unitarity is preserved. All in all, while the relative size of the corrections at two loops with respect to one loop depends greatly on the choices of parameters, we can take as a typical estimate of the relative size of the two-loop contributions to be $\mathscr {O}(\sim 20\%)$ of the one-loop ones – noting the caveat for the IDM and the HSM that large values of the scalar quartic couplings $\lambda _2$ or $\lambda _S$ respectively can result in numerically important additional contributions. For the 2HDM, we have also found that the largest possible deviations of the Higgs trilinear coupling with respect to its SM prediction occur for low $\tan \beta$ and intermediate masses – i.e. $M_\varPhi$ between 500 and 900 GeV. Finally, we have carried out a preliminary estimate of the theoretical uncertainty associated with our two-loop computation of ${\hat{\lambda }}_{hhh}$, taking this time the example of the IDM, and we found a conservative estimate of about 5% (of the total result).

Appendices

Notations and loop functions

This appendix summarises the notations and definitions of loop functions employed throughout this paper.

First of all, the loop factor is defined as

$$\begin{aligned} \kappa =\frac{1}{16\pi ^2}\,. \end{aligned}$$

(A.1)

We will make use of a number of loop integrals, which we define in $d=4-2\epsilon$ dimensions, with the shorthand notations

$$\begin{aligned} \hat{\alpha }\equiv 16\pi ^2 \frac{\mu ^{2\epsilon }}{i(2\pi )^d}\text { and }\int _{q}=\hat{\alpha } \int d^dq\,, \end{aligned}$$

(A.2)

with $\mu$ denoting the regularisation scale. We also use the definition

$$\begin{aligned} {{\,\mathrm{\overline{\text {log}}}\,}}x\equiv \log \frac{x}{Q^2}\,, \end{aligned}$$

(A.3)

where Q is the renormalisation scale, defined as $Q^2=4\pi e^{-\gamma _E}\mu ^2$. At one-loop order, we encounter the two integrals

$$\begin{aligned} \mathbf {A}(x)&\equiv -\int _k \frac{1}{k^2-x}\,,\nonumber \\ \mathbf {B}(p^2,x,y)&\equiv \int _k\frac{1}{(k^2-x)((p-k)^2-y)}\,, \end{aligned}$$

(A.4)

and from their finite (and $\epsilon$-independent) part, we obtain the usual Passarino-Veltmann functions [109] as

$$\begin{aligned} A(x)&\equiv \lim _{\epsilon \rightarrow 0}\bigg [\mathbf {A}(x)+\frac{x}{\epsilon }\bigg ]=x({{\,\mathrm{\overline{\text {log}}}\,}}x-1)\,,\nonumber \\ B(p^2,x,y)&\equiv \lim _{\epsilon \rightarrow 0}\bigg [\mathbf {B}(p^2,x,y)-\frac{1}{\epsilon }\bigg ]\nonumber \\&=-{{\,\mathrm{\overline{\text {log}}}\,}}p^2-f_B(x_+)-f_B(x_-), \end{aligned}$$

(A.5)

where

$$\begin{aligned} f_B(x)&=\log (1-x)-x\log \left( 1-\frac{1}{x}\right) -1\,,\nonumber \\ x_\pm&=\frac{p^2+x+y\pm \sqrt{(p^2+x+y)^2-4p^2x}}{2p^2}\,. \end{aligned}$$

(A.6)

The following limits are also useful

$$\begin{aligned} B(0,x,y)&=-\frac{A(x)-A(y)}{x-y}\,,\nonumber \\ B(x,0,0)&=2-{{\,\mathrm{\overline{\text {log}}}\,}}x\,,\nonumber \\ B(0,x,x)&=-{{\,\mathrm{\overline{\text {log}}}\,}}x\,,\nonumber \\ B(x,0,x)&=2-{{\,\mathrm{\overline{\text {log}}}\,}}x\,,\nonumber \\ B(x,x,x)&=2-\frac{\pi }{\sqrt{3}}-{{\,\mathrm{\overline{\text {log}}}\,}}x\,,\nonumber \\ B(x,0,y)&=2-{{\,\mathrm{\overline{\text {log}}}\,}}y+\left( \frac{y}{x}-1\right) \log \left( 1-\frac{x}{y}\right) \,. \end{aligned}$$

(A.7)

For loop diagrams involving fermions, we also use the function $B_1$ given by

$$\begin{aligned} B_1(p^2,x,y)=\frac{1}{2p^2}\big [A(x)-A(y)+(p^2+x-y)B_0(p^2,x,y)\big ]\,. \end{aligned}$$

(A.8)

Because we neglect all external momenta, we only need one master integral at two loops, namely the sunrise integral (see e.g. [124])

$$\begin{aligned} \mathbf {I}(x,y,z)\equiv - \int _{k_1}\int _{k_2} \frac{1}{(k_1^2-x)(k_2^2-y)((k_1+k_2)^2-z)}\,. \end{aligned}$$

(A.9)

Its finite, $\epsilon$-independent, part is obtained as [111]

$$\begin{aligned} I(x,y,z)\equiv&\ \lim _{\epsilon \rightarrow 0}\bigg [\mathbf {I}(x,y,z)-\frac{1}{\epsilon }[\mathbf {A}(x) +\mathbf {A}(y)+\mathbf {A}(z)]\nonumber \\&-\frac{1}{2}\left( \frac{1}{\epsilon ^2} -\frac{1}{\epsilon }\right) [x+y+z]\bigg ]\,. \end{aligned}$$

(A.10)

Several expressions for I – only differing by dilogarithm identities – can be found in numerous references in the literature, see e.g. Refs. [111, 121, 144]. Furthermore, a number of useful limits can be found in Refs. [111, 120], among which

$$\begin{aligned} I(0,0,x)\equiv&\ -\frac{1}{2} x {{\,\mathrm{\overline{\text {log}}}\,}}^2 x + 2 x {{\,\mathrm{\overline{\text {log}}}\,}}x - \frac{5}{2} x - \frac{\pi ^2}{6} x\,,\nonumber \\ I(0,x,x)\equiv&\ -x {{\,\mathrm{\overline{\text {log}}}\,}}^2 x + 4 x {{\,\mathrm{\overline{\text {log}}}\,}}x - 5 x\,,\nonumber \\ I(x,x,x)\equiv&\ \frac{3}{2} x (- {{\,\mathrm{\overline{\text {log}}}\,}}^2 x + 4 {{\,\mathrm{\overline{\text {log}}}\,}}x -5 + c_{xxx})\,, \end{aligned}$$

(A.11)

where $c_{xxx}$ is a numerical constant defined as

$$\begin{aligned} c_{xxx}\equiv - \frac{i}{\sqrt{3}} \left[ \frac{\pi ^2}{9} - 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{i\sqrt{3}}{2}\right) \right] \approx 2.3439\,, \end{aligned}$$

(A.12)

$\mathrm {Li}_2$ being the dilogarithm function.

Complete expressions for the 2HDM

We give in this appendix some expressions for the 2HDM that we have deemed too long for the discussion in the main text.

1.1 Derivatives of $V^{(2)}_{FFS}$ for degenerate BSM scalar masses

The complete expressions for the derivatives (with the operators $\mathscr {D}_3$ and $\mathscr {D}_4$) of the 2HDM effective potential diagrams involving both the top quark and BSM scalars are:

$$\begin{aligned}&\mathscr {D}_3V^{(2)}_{FFS}=\,\frac{48m_t^2\cot ^2\beta }{v^5} \Bigg \{-\frac{2(m_\varPhi ^2-m_t^2)(m_\varPhi ^2-M^2)^3}{m_\varPhi ^2(m_\varPhi ^2-4m_t^2)}\nonumber \\&\quad +\frac{(28 M^6 - 72 M^4 m_\varPhi ^2 + 72 M^2 m_\varPhi ^4 - 31 m_\varPhi ^6) m_t^4 + 28 m_\varPhi ^4 m_t^6}{m_\varPhi ^4(m_\varPhi ^2-4m_t^2)}\nonumber \\&\quad +\bigg [12M^4- 6 M^2 (2 m_\varPhi ^2 - 3 m_t^2)+ 4 m_\varPhi ^2 (m_\varPhi ^2 - 3 m_t^2)\nonumber \\&\quad +\frac{2 M^6 (-2 m_\varPhi ^8 + 15 m_\varPhi ^6 m_t^2 - 24 m_\varPhi ^4 m_t^4 - 10 m_\varPhi ^2 m_t^6 + 12 m_t^8)) }{m_\varPhi ^4 (m_\varPhi ^2 - m_t^2) (m_\varPhi ^2 - 4 m_t^2)^2}\bigg ]{{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\nonumber \\&\quad +6m_t^4\bigg [-1+\frac{ M^6 (m_\varPhi ^4 - 10 m_\varPhi ^2 m_t^2 + 12 m_t^4) }{m_\varPhi ^4 (m_\varPhi ^2 - m_t^2) (m_\varPhi ^2 - 4 m_t^2)^2}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2 \nonumber \\&\quad -\frac{24M^6m_t^8}{m_\varPhi ^5(m_\varPhi ^2-4m_t^2)^{5/2}}\bigg [\frac{\pi ^2}{3} -\log ^2\frac{m_t^2}{m_\varPhi ^2}\nonumber \\&\quad +2\log ^2\left( \frac{1}{2}-\frac{1}{2} \sqrt{1-\frac{4m_t^2}{m_\varPhi ^2}}\right) \nonumber \\&\quad -4\mathrm {Li}_2\left( \frac{1}{2}-\frac{1}{2}\sqrt{1-\frac{4m_t^2}{m_\varPhi ^2}} \right) \bigg ]\Bigg \}\,, \end{aligned}$$

(B.1)

and,

$$\begin{aligned} \mathscr {D}_4V^{(2)}_{FFS}&=\frac{192m_t^2\cot ^2\beta }{v^6}\Bigg \{-4 m_t^2(3 m_\varPhi ^2 - m_t^2) -3 M^2 (m_\varPhi ^2 - 9 m_t^2)+\frac{ 9 M^4}{m_\varPhi ^2} (m_\varPhi ^2 - m_t^2)\nonumber \\&-\frac{ M^6 (9 m_\varPhi ^4 - 21 m_\varPhi ^2 m_t^2 - 62 m_t^4)}{m_\varPhi ^4(m_\varPhi ^2 - 4 m_t^2)} -\frac{ M^8 (-3 m_\varPhi ^8 + 18 m_\varPhi ^6 m_t^2 + 12 m_\varPhi ^4 m_t^4 - 160 m_\varPhi ^2 m_t^6 + 124 m_t^8)}{m_\varPhi ^6 (m_\varPhi ^2 - m_t^2)(m_\varPhi ^2 - 4 m_t^2)^2 }\nonumber \\&+\bigg [4 (m_\varPhi ^4 - 3 m_\varPhi ^2 m_t^2) + 5 M^2 (-2 m_\varPhi ^2 + 3 m_t^2) +6 M^4 + \frac{M^6 (m_\varPhi ^4 - 2 m_\varPhi ^2 m_t^2 - 2 m_t^4) (2 m_\varPhi ^4 - 11 m_\varPhi ^2 m_t^2 + 6 m_t^4)}{m_\varPhi ^4 (m_\varPhi ^2 - m_t^2) (m_\varPhi ^2 - 4 m_t^2)^2} \nonumber \\&+ \frac{M^8 (-2 m_\varPhi ^{12} + 22 m_\varPhi ^{10} m_t^2 - 71 m_\varPhi ^8 m_t^4 + 44 m_\varPhi ^6 m_t^6 + 194 m_\varPhi ^4 m_t^8 - 280 m_\varPhi ^2 m_t^{10} + 120 m_t^{12})}{m_\varPhi ^6 (m_\varPhi ^2 - 4 m_t^2)^3 (m_\varPhi ^2 - m_t^2)^2}\bigg ]{{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\nonumber \\&+3m_t^4\bigg [-2 - \frac{M^6 (m_\varPhi ^4 - 10 m_\varPhi ^2 m_t^2 + 12 m_t^4}{m_\varPhi ^4 (m_\varPhi ^2 - 4 m_t^2)^2 (m_\varPhi ^2 - m_t^2)} \nonumber \\&+ \frac{M^8 (3 m_\varPhi ^8 - 46 m_\varPhi ^6 m_t^2 + 162 m_\varPhi ^4 m_t^4 - 216 m_\varPhi ^2 m_t^6 + 88 m_t^8)}{m_\varPhi ^6 (m_\varPhi ^2 - 4 m_t^2)^3 (m_\varPhi ^2 - m_t^2)^2}\bigg ]{{\,\mathrm{\overline{\text {log}}}\,}}m_t^2\nonumber \\&+\frac{12 M^6 m_t^8 (m_\varPhi ^4 - 4 m_\varPhi ^2 m_t^2 - 5 M^2 (m_\varPhi ^2 - 2 m_t^2))}{m_\varPhi ^7 (m_\varPhi ^2 - 4 m_t^2)^{7/2}} \times \bigg [\frac{\pi ^2}{3}-\log ^2\frac{m_t^2}{m_\varPhi ^2}+2\log ^2\left( \frac{1}{2}-\frac{1}{2}\sqrt{1-\frac{4m_t^2}{m_\varPhi ^2}}\right) \nonumber \\&-4\mathrm {Li}_2\left( \frac{1}{2}-\frac{1}{2}\sqrt{1-\frac{4m_t^2}{m_\varPhi ^2}}\right) \bigg ]\Bigg \}\,. \end{aligned}$$

(B.2)

When working beyond $\mathscr {O}(M_t^4)$, the scheme translation also requires new BSM contributions to the one-loop top-quark self-energy. Indeed, with respect to its SM expression – given in equation (5.4) – it receives the following new one-loop terms

$$\begin{aligned}&\varPi _{tt}^{(1)}(p^2)\Big |^\text {2HDM}\nonumber \\&\quad =-\frac{2m_t^4\cot ^2\beta }{v^2} \Big \{p^2\big [B_1(p^2,m_t^2,m_H^2)+B_1(p^2,m_t^2,m_A^2)\nonumber \\&\qquad +B_1(p^2,0,m_{H^\pm }^2)\big ]\nonumber \\&\qquad +m_t^2\big [B_0(p^2,m_t^2,m_H^2)-B_0(p^2,m_t^2,m_A^2)\big ]\Big \}\,. \end{aligned}$$

(B.3)

1.2 Results for arbitrary BSM scalar masses

We give in this section results for the 2HDM, when not assuming the masses of the BSM scalars to be degenerate as in the main text.

First, the dominant BSM contributions to the 2HDM effective potential at two loops read

$$\begin{aligned} V^{(2)}_{SSS}(h)=&-\sum _{\varPhi =H,A,H^\pm }\frac{n_\varPhi (M^2-m_\varPhi ^2)^2(v+h)^2}{v^4}\nonumber \\&\times I(0,m_\varPhi ^2(h),m_\varPhi ^2(h))\nonumber \\&-\sum _{\varPhi =A,H^\pm }\frac{n_\varPhi (M^2-m_H^2)^2\cot ^22\beta (v+h)^2}{v^4}\nonumber \\&\times I(m_H^2(h),m_\varPhi ^2(h),m_\varPhi ^2(h))\nonumber \\&-\frac{3(M^2-m_H^2)^2\cot ^22\beta (v+h)^2}{v^4}\nonumber \\&\times I(m_H^2(h),m_H^2(h),m_H^2(h))\nonumber \\&-\frac{(v+h)^2}{2v^4}\bigg [(m_H^2-m_A^2)^2I(0,m_H^2(h),m_A^2(h))\nonumber \\&+2(m_H^2-m_{H^\pm }^2)^2I(0,m_H^2(h),m_{H^\pm }^2(h))\nonumber \\&+2(m_A^2-m_{H^\pm }^2)^2I(0,m_A^2(h),m_{H^\pm }^2(h))\bigg ] \end{aligned}$$

(B.4)

$$\begin{aligned} V^{(2)}_{SS}(h)=&-\frac{(M^2 - m_H^2) \cot ^22\beta }{v^2} \bigg [\frac{3}{2} A(m_H^2(h))^2 \nonumber \\&+ \frac{3}{2} A(m_A^2(h))^2 + 4 A(m_{H^\pm }^2(h))^2\nonumber \\&+ A(m_A^2(h)) A(m_H^2(h)) \nonumber \\&+ 2 A(m_H^2(h)) A(m_{H^\pm }^2(h))\nonumber \\&+ 2 A(m_A^2(h)) A(m_{H^\pm }^2(h))\bigg ] \end{aligned}$$

(B.5)

$$\begin{aligned} V^{(2)}_{FFS}(h)=&-\frac{3}{2}y_t^2c_\beta ^2\bigg [2A(m_t^2(h))A(m_H^2(h))\nonumber \\&-A(m_t^2(h))^2-(4m_t^2(h)-m_H^2(h))\nonumber \\&\times I(m_H^2(h),m_t^2(h),m_t^2(h))\bigg ]\nonumber \\&-\frac{3}{2}y_t^2c_\beta ^2\bigg [2A(m_t^2(h))A(m_A^2(h))\nonumber \\&-A(m_t^2(h))^2 +m_A^2(h)I(m_A^2(h),m_t^2(h),m_t^2(h))\bigg ]\nonumber \\&-3y_t^2c_\beta ^2\bigg [A(m_t^2(h))A(m_{H^\pm }^2(h))\nonumber \\&-(m_t^2(h)-m_{H^\pm }^2(h)) I(0,m_{H^\pm }^2(h),m_t^2(h))\bigg ] \end{aligned}$$

(B.6)

Then, applying the operators $\mathscr {D}_3$ and $\mathscr {D}_4$ to these effective potential contributions, we obtain the leading two-loop corrections to $\lambda _{hhh}$ and $\lambda _{hhhh}$, in terms of $\overline{{\mathrm{MS}}}$-renormalised parameters, as

$$\begin{aligned} \delta ^{(2)}\lambda _{hhh}=&\sum _{\varPhi =H,A,H^\pm }\frac{16n_\varPhi m_\varPhi ^4}{v^5} \left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^4 \big [-2 M^2 - m_\varPhi ^2 + (M^2 + 2 m_\varPhi ^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\big ] \nonumber \\&+\sum _{\varPhi =A,H^\pm }\frac{8n_\varPhi m_H^4\cot ^22\beta }{v^5}\left( 1-\frac{M^2}{m_H^2}\right) ^2 \nonumber \\&\times \Bigg \{-m_H^2 - 2 m_\varPhi ^2+ \frac{3 M^4 (2 m_H^4 - 7 m_H^2 m_\varPhi ^2 - 4 m_\varPhi ^4)}{m_H^2 m_\varPhi ^2 (m_H^2 - 4 m_\varPhi ^2)}-\frac{2 M^6 (m_H^6 - 5 m_H^4 m_\varPhi ^2 - 2 m_H^2 m_\varPhi ^4 - 3 m_\varPhi ^6) }{m_H^4 m_\varPhi ^4 (m_H^2 - 4 m_\varPhi ^2)} \nonumber \\&\quad +\bigg [-3 M^2 + 2 m_H^2 + \frac{M^6 (2 m_H^8 - 20 m_H^6 m_\varPhi ^2 + 31 m_H^4 m_\varPhi ^4 - 16 m_H^2 m_\varPhi ^6 + 12 m_\varPhi ^8)}{m_H^4 m_\varPhi ^4 (m_H^2 - 4 m_\varPhi ^2)^2}\bigg ]{{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+2 \bigg [-3 M^2 + 2 m_\varPhi ^2 + \frac{M^6 (2 m_H^2 + m_\varPhi ^2) (m_H^4 + 2 m_\varPhi ^4)}{m_H^4 m_\varPhi ^2 (m_H^2 - 4 m_\varPhi ^2)^2}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\nonumber \\&+\frac{6 M^6 (m_H^6 - 3 m_H^2 m_\varPhi ^4 + 2 m_\varPhi ^6) }{m_H^5 (m_H^2 - 4 m_\varPhi ^2)^{5/2}}\bigg [\frac{\pi ^2}{3} - \log ^2\frac{m_\varPhi ^2}{m_H^2} + 2 \log ^2\left( \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{4 m_\varPhi ^2}{m_H^2}}\right) \nonumber \\&- 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{4 m_\varPhi ^2}{m_H^2}}\right) \bigg ]\Bigg \} \nonumber \\&+\frac{72 m_H^4\cot ^22\beta }{v^5} \left( 1-\frac{M^2}{m_H^2}\right) ^4 \big [-2 M^2 - m_H^2 + (M^2 + 2 m_H^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\big ] +\frac{4(m_H^2 - m_A^2)^2}{v^5}\nonumber \\&\times \Bigg \{-\left( \frac{1}{m_H^2} + \frac{1}{m_A^2}\right) \left[ -3 M^4 + m_H^2 m_A^2 + M^6 \left( \frac{1}{m_H^2} + \frac{1}{m_A^2}\right) \right] + \bigg [\frac{M^6}{m_A^4} - 3 M^2 + 2 m_H^2\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+\bigg [\frac{M^6}{m_H^4} - 3 M^2 + 2 m_A^2\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2\Bigg \} +\frac{8(m_H^2 - m_{H^\pm }^2)^2}{v^5}\Bigg \{-\left( \frac{1}{m_H^2} + \frac{1}{m_{H^\pm }^2}\right) \nonumber \\&\times \left[ -3 M^4 + m_H^2 m_{H^\pm }^2 + M^6 \left( \frac{1}{m_H^2} + \frac{1}{m_{H^\pm }^2}\right) \right] + \bigg [\frac{M^6}{m_{H^\pm }^4} - 3 M^2 + 2 m_H^2\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+\bigg [\frac{M^6}{m_H^4} - 3 M^2 + 2 m_{H^\pm }^2\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\Bigg \} \nonumber \\&\quad +\frac{8(m_A^2 - m_{H^\pm }^2)^2}{v^5}\Bigg \{-\left( \frac{1}{m_A^2} + \frac{1}{m_{H^\pm }^2}\right) \left[ -3 M^4 + m_A^2 m_{H^\pm }^2 + M^6 \left( \frac{1}{m_A^2} + \frac{1}{m_{H^\pm }^2}\right) \right] \nonumber \\&+ \bigg [\frac{M^6}{m_{H^\pm }^4} - 3 M^2 + 2 m_A^2\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2+\bigg [\frac{M^6}{m_A^4} - 3 M^2 + 2 m_{H^\pm }^2\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\Bigg \}\\&+\frac{8m_H^2\cot ^22\beta }{v^5}\left( 1-\frac{M^2}{m_H^2}\right) \nonumber \\&\times \Bigg \{3m_H^4\left( 1-\frac{M^2}{m_H^2}\right) ^3\big [1+2{{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\big ] +3m_A^4\left( 1-\frac{M^2}{m_A^2}\right) ^3\big [1+2{{\,\mathrm{\overline{\text {log}}}\,}}m_A^2\big ] \end{aligned}$$

$$\begin{aligned}&+8m_{H^\pm }^4\left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^3\big [1+2{{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\big ] +m_H^2 m_A^2 \bigg [\left( 1 - \frac{M^2}{m_H^2}\right) ^3 + \left( 1 - \frac{M^2}{m_A^2}\right) ^3\bigg ]\nonumber \\&+\left( 1 - \frac{M^2}{m_H^2}\right) ^2 \big [2 m_H^2 m_A^2 + M^2 (m_A^2 - 3 m_H^2)\big ] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2 + \left( 1 - \frac{M^2}{m_A^2}\right) ^2 \big [2 m_H^2 m_A^2 + M^2 (m_H^2 - 3 m_A^2)\big ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+2m_H^2 m_{H^\pm }^2 \bigg [\left( 1 - \frac{M^2}{m_H^2}\right) ^3 + \left( 1 - \frac{M^2}{m_{H^\pm }^2}\right) ^3\bigg ] +2\left( 1 - \frac{M^2}{m_H^2}\right) ^2 \big [2 m_H^2 m_{H^\pm }^2 + M^2 (m_{H^\pm }^2 - 3 m_H^2)\big ] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\nonumber \\&+ 2\left( 1 - \frac{M^2}{m_{H^\pm }^2}\right) ^2 \big [2 m_H^2 m_{H^\pm }^2 + M^2 (m_H^2 - 3 m_{H^\pm }^2)\big ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2 +2m_A^2 m_{H^\pm }^2 \bigg [\left( 1 - \frac{M^2}{m_A^2}\right) ^3 + \left( 1 - \frac{M^2}{m_{H^\pm }^2}\right) ^3\bigg ]\nonumber \\&+2\left( 1 - \frac{M^2}{m_H^2}\right) ^2 \big [2 m_A^2 m_{H^\pm }^2 + M^2 (m_{H^\pm }^2 - 3 m_A^2)\big ] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2 + 2\left( 1 - \frac{M^2}{m_{H^\pm }^2}\right) ^2 \big [2 m_A^2 m_{H^\pm }^2 + M^2 (m_A^2 - 3 m_{H^\pm }^2)\big ] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2\bigg \}\nonumber \\&+\frac{24 m_t^2 \cot ^2\beta }{v^5}\Bigg \{ -m_H^4\left( 1-\frac{M^2}{m_H^2}\right) ^3 - 2 m_t^2 (m_H^2-4 M^2) \left( 1-\frac{M^2}{m_H^2}\right) ^2 + 10 m_t^4 \nonumber \\&+\bigg [ \frac{2m_H^6 }{m_H^2 - 4 m_t^2}\left( 1-\frac{M^2}{m_H^2}\right) ^3 + \frac{2m_H^2m_t^2}{m_H^2 - 4 m_t^2} (M^2 - 10 m_H^2) \left( 1-\frac{M^2}{m_H^2}\right) ^2 \nonumber \\&+ \frac{12m_t^4}{m_H^4(m_H^2 - 4 m_t^2)} (M^6 - 6 M^2 m_H^4 + 4 m_H^6) \bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2 + \frac{12 m_t^4}{m_H^4(m_H^2 - 4 m_t^2)}\big [M^6 - m_H^6 + 4 m_H^4 m_t^2 \big ]{{\,\mathrm{\overline{\text {log}}}\,}}m_t^2\nonumber \\&+\frac{12 M^6 m_t^6}{m_H^5 (m_H^2 - 4 m_t^2)^{3/2}} \bigg [\frac{\pi ^2}{3}-\log ^2\frac{m_t^2}{m_H^2} + 2 \log ^2\left( \frac{1}{2} - \frac{1}{2}\sqrt{1 - \frac{4 m_t^2}{m_H^2}}\right) - 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{1}{2}\sqrt{1 - \frac{4 m_t^2}{m_H^2}}\right) \bigg ]\Bigg \}\nonumber \\&+\frac{24 m_t^2 \cot ^2\beta }{v^5}\Bigg \{ -m_A^4\left( 1-\frac{M^2}{m_A^2}\right) ^3 - 6m_t^2m_A^2\left( 1-\frac{M^2}{m_A^2}\right) ^2 +2m_t^4\bigg [1 - \frac{4 M^6}{m_A^4(m_A^2 - 4 m_t^2)}\bigg ]\nonumber \\&+\bigg [ 2m_A^2 (m_A^2 - 2 m_t^2) - 6 M^2 (m_A^2 - m_t^2) + 6 M^4 - \frac{2M^6 (m_A^4 - 7 m_A^2 m_t^2 + 6 m_t^4)}{m_A^2 (m_A^2 - 4 m_t^2)^2}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2 \nonumber \\&+4m_t^4\bigg [1 -\frac{ M^6 (m_A^2 + 8 m_t^2)}{m_A^4 (m_A^2 - 4 m_t^2)^2}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2 -\frac{12 M^6 m_t^6}{m_A^3 (m_A^2 - 4 m_t^2)^{5/2}} \bigg [\frac{\pi ^2}{3}-\log ^2\frac{m_t^2}{m_A^2} \nonumber \\&+ 2\log ^2\left( \frac{1}{2} - \frac{1}{2}\sqrt{1 - \frac{4 m_t^2}{m_A^2}}\right) - 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{1}{2}\sqrt{1 - \frac{4 m_t^2}{m_A^2}}\right) \bigg ]\Bigg \} +\frac{48 m_t^2 \cot ^2\beta }{v^5}\Bigg \{- m_{H^\pm }^4 \left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^3\nonumber \\&- 2 m_t^2m_{H^\pm }^2\left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^3 + m_t^4 +\bigg [ \frac{2m_{H^\pm }^6 }{m_{H^\pm }^2 - m_t^2}\left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^3\nonumber \\&- \frac{6m_{H^\pm }^4 m_t^2}{ m_{H^\pm }^2 - m_t^2} \left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^2 + \frac{2( 2 m_{H^\pm }^2-3 M^2) m_t^4}{m_{H^\pm }^2 - m_t^2}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2 +\frac{2 m_t^4}{(m_{H^\pm }^2 - m_t^2)} \bigg [\frac{M^6}{m_{H^\pm }^4} - m_{H^\pm }^2 + m_t^2 \bigg ]{{\,\mathrm{\overline{\text {log}}}\,}}m_t^2\Bigg \}\,, \end{aligned}$$

(B.7)

and,

$$\begin{aligned} \delta ^{(2)}\lambda _{hhhh}=&\sum _{\varPhi =H,A,H^\pm } \frac{32n_\varPhi m_\varPhi ^2}{v^6}\left( 1-\frac{M^2}{m_\varPhi ^2}\right) ^4 \big [-5 M^4 - 4 M^2 m_\varPhi ^2 + (2 M^4 + 3 M^2 m_\varPhi ^2 + 4 m_\varPhi ^4) {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\big ]\nonumber \\&+\sum _{\varPhi =A,H^\pm }\frac{16n_\varPhi m_H^4\cot ^22\beta }{v^6}\left( 1-\frac{M^2}{m_H^2}\right) ^2 \times \Bigg \{-12 M^2 + 3 M^4 \left( \frac{1}{m_H^2} + \frac{2}{m_\varPhi ^2}\right) \nonumber \\&+ \frac{2 M^6 (5 m_H^6 - 21 m_H^4 m_\varPhi ^2 - 11 m_\varPhi ^6)}{m_H^4 m_\varPhi ^4 (m_H^2 - 4 m_\varPhi ^2)} - \frac{M^8 (4 m_H^{10} - 28 m_H^8 m_\varPhi ^2 + 60 m_H^6 m_\varPhi ^4 + 53 m_H^4 m_\varPhi ^6 - 14 m_H^2 m_\varPhi ^8 + 60 m_\varPhi ^{10})}{m_H^6 m_\varPhi ^6 (m_H^2 - 4 m_\varPhi ^2)^2}\nonumber \\&-\frac{6 M^6 (m_H^2 - m_\varPhi ^2)^2}{m_H^7 (m_H^2 - 4 m_\varPhi ^2)^{7/2}}\big [m_H^6 - 2 m_H^4 m_\varPhi ^2 - 8 m_H^2 m_\varPhi ^4 + 4 M^2 (2 m_H^4 + 2 m_H^2 m_\varPhi ^2 + 5 m_\varPhi ^4)\big ]\nonumber \\&\times \bigg [\frac{\pi ^2}{3} - \log ^2\frac{m_\varPhi ^2}{m_H^2} + 2 \log ^2\left( \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{4 m_\varPhi ^2}{m_H^2}}\right) - 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{4 m_\varPhi ^2}{m_H^2}}\right) \bigg ]\nonumber \\&+\bigg [-5 M^2 + 4 m_H^2 + \frac{ M^6 (-2 m_H^8 + 20 m_H^6 m_\varPhi ^2 - 31 m_H^4 m_\varPhi ^4 + 16 m_H^2 m_\varPhi ^6 - 12 m_\varPhi ^8)}{m_H^4 m_\varPhi ^4 (m_H^2 - 4 m_\varPhi ^2)^2} \nonumber \\&+ \frac{2 M^8 (2 m_H^{12} - 26 m_H^{10} m_\varPhi ^2 + 124 m_H^8 m_\varPhi ^4 - 131 m_H^6 m_\varPhi ^6 - 22 m_H^4 m_\varPhi ^8 + 86 m_H^2 m_\varPhi ^{10} - 60 m_\varPhi ^{12})}{m_H^6 m_\varPhi ^6 (m_H^2 - 4 m_\varPhi ^2)^3}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+2\bigg [-5 M^2 + 4 m_\varPhi ^2 - \frac{M^6 (2 m_H^2 + m_\varPhi ^2) (m_H^4 + 2 m_\varPhi ^4)}{m_H^4 m_\varPhi ^2 (m_H^2 - 4 m_\varPhi ^2)^2} \nonumber \\&+ \frac{2 M^8 (m_H^{10} - 14 m_H^8 m_\varPhi ^2 + 2 m_H^6 m_\varPhi ^4 + 5 m_H^4 m_\varPhi ^6 - 19 m_H^2 m_\varPhi ^8 - 2 m_\varPhi ^{10})}{m_H^6 m_\varPhi ^4 (m_H^2 - 4 m_\varPhi ^2)^3}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_\varPhi ^2\Bigg \}\nonumber \\&+\frac{144 m_H^2\cot ^22\beta }{v^6}\left( 1-\frac{M^2}{m_H^2}\right) ^4 \big [-5 M^4 - 4 M^2 m_H^2 + (2 M^4 + 3 M^2 m_H^2 + 4 m_H^4) {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\big ]\nonumber \\&+\frac{8(m_H^2 - m_A^2)^2}{v^6}\Bigg \{ -8 M^2 + 3 M^4 \left( \frac{1}{m_H^2} + \frac{1}{m_A^2}\right) + M^6 \left( \frac{5}{m_H^4} + \frac{2}{m_H^2 m_A^2} + \frac{5}{m_A^4}\right) \nonumber \\&- M^8 \left( \frac{1}{m_H^2} + \frac{1}{m_A^2}\right) \left( \frac{2}{m_H^4} + \frac{1}{m_H^2 m_A^2} + \frac{2}{m_A^4}\right) +\left[ \frac{2 M^8}{m_A^6} - \frac{M^6}{m_A^4} - 5 M^2 + 4 m_H^2\right] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+\left[ \frac{2M^8}{m_H^6} -\frac{M^6}{m_H^4} - 5 M^2 + 4 m_A^2\right] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2\Bigg \} +\frac{16(m_H^2 - m_{H^\pm }^2)^2}{v^6}\Bigg \{ -8 M^2 + 3 M^4 \left( \frac{1}{m_H^2} + \frac{1}{m_{H^\pm }^2}\right) \nonumber \\&+ M^6 \left( \frac{5}{m_H^4} + \frac{2}{m_H^2 m_{H^\pm }^2} + \frac{5}{m_{H^\pm }^4}\right) - M^8 \left( \frac{1}{m_H^2} + \frac{1}{m_{H^\pm }^2}\right) \left( \frac{2}{m_H^4} + \frac{1}{m_H^2 m_{H^\pm }^2} + \frac{2}{m_{H^\pm }^4}\right) \nonumber \\&+\left[ \frac{2 M^8}{m_{H^\pm }^6} - \frac{M^6}{m_{H^\pm }^4} - 5 M^2 + 4 m_H^2\right] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&+\left[ \frac{2M^8}{m_H^6} -\frac{M^6}{m_H^4} - 5 M^2 + 4 m_{H^\pm }^2\right] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\Bigg \} +\frac{16(m_A^2 - m_{H^\pm }^2)^2}{v^6}\Bigg \{ -8 M^2 + 3 M^4 \left( \frac{1}{m_A^2} + \frac{1}{m_{H^\pm }^2}\right) \nonumber \\&+ M^6 \left( \frac{5}{m_A^4} + \frac{2}{m_A^2 m_{H^\pm }^2} + \frac{5}{m_{H^\pm }^4}\right) \end{aligned}$$

$$\begin{aligned}&- M^8 \left( \frac{1}{m_A^2} + \frac{1}{m_{H^\pm }^2}\right) \left( \frac{2}{m_A^4} + \frac{1}{m_A^2 m_{H^\pm }^2} + \frac{2}{m_{H^\pm }^4}\right) +\left[ \frac{2 M^8}{m_{H^\pm }^6} - \frac{M^6}{m_{H^\pm }^4} - 5 M^2 + 4 m_A^2\right] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2\nonumber \\&+\left[ \frac{2M^8}{m_A^6} -\frac{M^6}{m_A^4} - 5 M^2 + 4 m_{H^\pm }^2\right] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\Bigg \}+\frac{16m_H^2\cot ^22\beta }{v^6}\left( 1-\frac{M^2}{m_H^2}\right) \Bigg \{3m_H^2\left( 1-\frac{M^2}{m_H^2}\right) ^3\big [-M^2 + 4 m_H^2\nonumber \\&+ 2 (M^2 + 2 m_H^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\big ]+3m_A^2\left( 1-\frac{M^2}{m_A^2}\right) ^3\big [-M^2 + 4 m_A^2 + 2 (M^2 + 2 m_A^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2\big ]+8m_{H^\pm }^2\left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^3\nonumber \\&\times \big [-M^2 + 4 m_{H^\pm }^2 + 2 (M^2 + 2 m_{H^\pm }^2) {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2\big ]+\left( 1-\frac{M^2}{m_H^2}\right) \left( 1-\frac{M^2}{m_A^2}\right) \bigg [16 m_H^2 m_A^2 - 7 M^2 (m_H^2 + m_A^2) - 4M^4\bigg (\frac{m_H^2}{m_A^2} - \frac{3}{2} + \frac{m_A^2}{m_H^2}\bigg )\bigg ]\nonumber \\&+ \left( 1-\frac{M^2}{m_A^2}\right) ^2 \bigg [4 m_H^2 m_A^2 + M^2 (3 m_H^2-5 m_A^2) + 2M^4 \bigg (\frac{m_H^2}{m_A^2}-2\bigg )\bigg ][{{\,\mathrm{\overline{\text {log}}}\,}}m_H^2-1]+ \left( 1-\frac{M^2}{m_H^2}\right) ^2 \bigg [4 m_H^2 m_A^2 + M^2 (3 m_A^2-5 m_H^2)\nonumber \\&+ 2M^4 \bigg (\frac{m_A^2}{m_H^2}-2\bigg )\bigg ][{{\,\mathrm{\overline{\text {log}}}\,}}m_A^2-1]+2\left( 1-\frac{M^2}{m_H^2}\right) \left( 1-\frac{M^2}{m_{H^\pm }^2}\right) \bigg [16 m_H^2 m_{H^\pm }^2 - 7 M^2 (m_H^2 + m_{H^\pm }^2) - 4M^4\bigg (\frac{m_H^2}{m_{H^\pm }^2} - \frac{3}{2} + \frac{m_{H^\pm }^2}{m_H^2}\bigg )\bigg ]\nonumber \\&+ 2\left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^2 \bigg [4 m_H^2 m_{H^\pm }^2 + M^2 (3 m_H^2-5 m_{H^\pm }^2) + 2M^4 \bigg (\frac{m_H^2}{m_{H^\pm }^2}-2\bigg )\bigg ][{{\,\mathrm{\overline{\text {log}}}\,}}m_H^2-1]+ 2\left( 1-\frac{M^2}{m_H^2}\right) ^2 \nonumber \\&\times \bigg [4 m_H^2 m_{H^\pm }^2 + M^2 (3 m_{H^\pm }^2-5 m_H^2) + 2M^4 \bigg (\frac{m_{H^\pm }^2}{m_H^2}-2\bigg )\bigg ][{{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2-1]+2\left( 1-\frac{M^2}{m_A^2}\right) \left( 1-\frac{M^2}{m_{H^\pm }^2}\right) \bigg [16 m_A^2 m_{H^\pm }^2\nonumber \\&-7 M^2 (m_A^2 + m_{H^\pm }^2) - 4M^4\bigg (\frac{m_A^2}{m_{H^\pm }^2} - \frac{3}{2} + \frac{m_{H^\pm }^2}{m_A^2}\bigg )\bigg ]+ 2\left( 1-\frac{M^2}{m_{H^\pm }^2}\right) ^2 \bigg [4 m_A^2 m_{H^\pm }^2 + M^2 (3 m_A^2-5 m_{H^\pm }^2) + 2M^4\nonumber \\&\times \bigg (\frac{m_A^2}{m_{H^\pm }^2}-2\bigg )\bigg ][{{\,\mathrm{\overline{\text {log}}}\,}}m_A^2-1]+ 2\left( 1-\frac{M^2}{m_A^2}\right) ^2 \bigg [4 m_A^2 m_{H^\pm }^2 + M^2 (3 m_{H^\pm }^2-5 m_A^2) + 2M^4 \bigg (\frac{m_{H^\pm }^2}{m_A^2}-2\bigg )\bigg ][{{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2-1]\Bigg \}\nonumber \\&+\frac{12 m_t^2 \cot ^2\beta }{v^6}\Bigg \{32 m_t^2 (m_t^2-2 m_H^2 ) + 4 M^2 (44 m_t^2-3 m_H^2 ) + 36M^4 \left( 1 - \frac{2 m_t^2}{m_H^2}\right) - 4 M^6 \left( \frac{9}{m_H^2} + \frac{32 m_t^2}{m_H^4}\right) \nonumber \\&+ M^8 \left( \frac{9}{m_H^4} + \frac{76 m_t^2}{m_H^6} + \frac{3}{m_H^2 (m_H^2 - 4 m_t^2)}\right) +8 \bigg [2 m_H^2(m_H^2 - 6 m_t^2) - 5 M^2 (m_H^2 - 3 m_t^2) + 3 M^4\nonumber \\&+ \frac{M^6 (m_H^2 - 3 m_t^2) (m_H^2 + 2 m_t^2)}{m_H^4(m_H^2 - 4 m_t^2)} - \frac{M^8 }{(m_H^2 - 4 m_t^2)^2} \left( 1 - \frac{2 m_t^2}{m_H^2} - \frac{14 m_t^4}{m_H^4} + \frac{60 m_t^6}{m_H^6}\right) \bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_H^2\nonumber \\&-48 m_t^4\bigg [2 + \frac{M^6}{m_H^4 (m_H^2 - 4 m_t^2)} - \frac{3 M^8(m_H^2 - 2 m_t^2)}{m_H^6 (m_H^2 - 4 m_t^2)^2} \bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2-\frac{48 M^6 m_t^6 }{m_H^3 (m_H^2 - 4 m_t^2)^{5/2} }\bigg [1 - \frac{4 m_t^2}{m_H^2} - \frac{M^2}{m_H^4} (4 m_H^2 - 10 m_t^2)\bigg ]\nonumber \\&\times \bigg [\frac{\pi ^2}{3}-\log ^2\frac{m_t^2}{m_H^2} + 2 \log ^2\left( \frac{1}{2} - \frac{1}{2}\sqrt{1 - \frac{4 m_t^2}{m_H^2}}\right) - 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{4 m_t^2}{m_H^2}}\right) \bigg ]\Bigg \}+\frac{12 m_t^2 \cot ^2\beta }{v^6}\Bigg \{4 \bigg [8 m_t^2 (m_t^2-2 m_A^2)\nonumber \\&+ M^2 (28 m_t^2-3 m_A^2) + 3 M^4 \left( 3 - \frac{2 m_t^2}{m_A^2}\right) - \frac{M^6 }{m_A^2 - 4 m_t^2}\left( 9 - \frac{28 m_t^2}{ m_A^2} -\frac{ 40 m_t^4}{m_A^4}\right) + \frac{M^8 }{(m_A^2 - 4 m_t^2)^2}\left( 3 -\frac{ 22 m_t^2}{m_A^2} + \frac{4 m_t^4}{m_A^4} +\frac{ 64 m_t^6}{m_A^6}\right) \bigg ]\nonumber \\&+8 \bigg [2 m_A^2 (m_A^2 - 2 m_t^2) + 5 M^2 (-m_A^2 + m_t^2) + 3 M^4 + \frac{M^6m_A^2 }{(m_A^2 - 4 m_t^2)^2}\left( 1 - \frac{7 m_t^2}{m_A^2} + \frac{6 m_t^4}{m_A^4}\right) \nonumber \\&- \frac{M^8 m_A^2 }{(m_A^2 - 4 m_t^2)^3}\left( 1- \frac{10 m_t^2}{m_A^2} + \frac{18 m_t^4}{m_A^4} - \frac{36 m_t^6}{m_A^6}\right) \bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_A^2+16 m_t^4 \bigg [2 +\frac{ M^6 (m_A^2 + 8 m_t^2)}{m_A^4 (m_A^2 - 4 m_t^2)^2} \end{aligned}$$

$$\begin{aligned}&-\frac{M^8 (3 m_A^4 + 34 m_A^2 m_t^2 - 64 m_t^4)}{m_A^6 (m_A^2 - 4 m_t^2)^3)} \bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2+\frac{48 M^6 m_t^6 }{m_A (m_A^2 - 4 m_t^2)^{7/2} }\bigg [1 - \frac{4 m_t^2}{m_A^2} - \frac{M^2}{m_A^4} (4 m_A^2 - 6 m_t^2)\bigg ]\nonumber \\&\times \bigg [\frac{\pi ^2}{3}-\log ^2\frac{m_t^2}{m_A^2} + 2 \log ^2\left( \frac{1}{2} - \frac{1}{2}\sqrt{1 - \frac{4 m_t^2}{m_A^2}}\right) - 4 \mathrm {Li}_2\left( \frac{1}{2} - \frac{1}{2} \sqrt{1 - \frac{4 m_t^2}{m_A^2}}\right) \bigg ]\Bigg \}\nonumber \\&-\frac{96 m_t^2 \cot ^2\beta }{v^6}\Bigg \{8 m_{H^\pm }^2 m_t^2 + 3 (m_{H^\pm }^2 - 6 m_t^2) M^2 - 3 M^4\left( 3 - \frac{2 m_t^2}{m_{H^\pm }^2}\right) +\frac{ M^6}{m_{H^\pm }^2}\left( 9 + \frac{10 m_t^2}{m_{H^\pm }^2}\right) \nonumber \\&- \frac{ M^8}{m_{H^\pm }^2(m_{H^\pm }^2 - m_t^2)}\left( 3 + \frac{3 m_t^2}{m_{H^\pm }^2} - \frac{4 m_t^4}{m_{H^\pm }^4}\right) -2 \bigg [ 2 m_{H^\pm }^2 (m_{H^\pm }^2 - 2 m_t^2) - 5 M^2 (m_{H^\pm }^2 - m_t^2) + 3 M^4 + \frac{M^6}{m_{H^\pm }^2 - m_t^2}\nonumber \\&-\frac{M^8}{(m_{H^\pm }^2 - m_t^2)^2)}\bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_{H^\pm }^2+2 m_t^4 \bigg [2+ \frac{M^6}{m_{H^\pm }^4 (m_{H^\pm }^2 - m_t^2)} + \frac{M^8 ( 2 m_t^2-3 m_{H^\pm }^2 )}{m_{H^\pm }^6 (m_{H^\pm }^2 - m_t^2)^2} \bigg ] {{\,\mathrm{\overline{\text {log}}}\,}}m_t^2 \Bigg \}\,. \end{aligned}$$

(B.8)

Although somewhat tedious, the conversion of these results into the OS scheme poses no conceptual problem, and can be performed using the expressions in equations (3.19), (5.11), (B.3), as well as

$$\begin{aligned} \kappa ^{-1}\frac{\delta ^\text {OS}M^2}{M^2}=&\ \frac{2m_H^2\cot ^22\beta }{v^2}\left( 1-\frac{M^2}{m_H^2}\right) \nonumber \\&\times \bigg [\frac{A(m_H^2)}{m_H^2}+\frac{3A(m_A^2)}{m_A^2}+\frac{2A(m_{H^\pm }^2)}{m_{H^\pm }^2}\bigg ]\nonumber \\&-\frac{3m_t^2\cot ^2\beta }{2v^2}\bigg [B_0(m_H^2,m_t^2,m_t^2)\nonumber \\&+B_0(m_A^2,m_t^2,m_t^2)+2B_0(m_{H^\pm }^2,0,m_t^2)\bigg ]\,. \end{aligned}$$

(B.9)

Intermediate functions for the general expressions

1.1 Non-degenerate mass case

While expressions for derivatives of the effective potential simplify greatly when some masses are degenerate, in the general case, they are more involved and we have defined the following functions to reduce the length of our results.

$$\begin{aligned} \varDelta _{xyz}\equiv&\, (m_x - m_y - m_z) (m_x + m_y - m_z)\nonumber \\&\quad (m_x - m_y + m_z) (m_x + m_y + m_z)\,, \nonumber \\ \omega _{xyz}\equiv&\, m_x^4 - 2 m_y^4 + m_z^4 + m_x^2 m_y^2 - 2 m_x^2 m_z^2 + m_y^2 m_z^2\,,\nonumber \\ \varXi _{xyz}\equiv&\, m_x^6 + m_y^6 + m_z^6 - m_x^4 m_y^2 - m_x^2 m_y^4 - m_x^4 m_z^2\nonumber \\&- m_x^2 m_z^4 - m_y^4 m_z^2 - m_y^2 m_z^4 + 6 m_x^2 m_y^2 m_z^2\,,\nonumber \\ \chi _{xyz}\equiv&\, m_x^6 - 3 m_x^4 m_y^2 + 3 m_x^2 m_y^4 - m_y^6 - 3 m_x^4 m_z^2\nonumber \\&- 2 m_x^2 m_y^2 m_z^2 + m_y^4 m_z^2 + 3 m_x^2 m_z^4 + m_y^2 m_z^4 - m_z^6\,,\nonumber \\ \end{aligned}$$

$$\begin{aligned} \varTheta _{xyz}\equiv&\, m_x^8 - 4 m_x^6 m_y^2 + 6 m_x^4 m_y^4 - 4 m_x^2 m_y^6 + m_y^8\nonumber \\&- 4 m_x^6 m_z^2 + 4 m_x^4 m_y^2 m_z^2 - 8 m_x^2 m_y^4 m_z^2 + 8 m_y^6 m_z^2\nonumber \\&+ 6 m_x^4 m_z^4 - 8 m_x^2 m_y^2 m_z^4 - 18 m_y^4 m_z^4\nonumber \\&- 4 m_x^2 m_z^6 + 8 m_y^2 m_z^6 + m_z^8\,,\nonumber \\ \theta _{xyz}\equiv&\, (m_x^2 - 2 m_y^2) (m_x^2 - m_y^2)^4 + (-5 m_x^8 + 10 m_x^6 m_y^2\nonumber \\&+ 8 m_x^4 m_y^4 - 14 m_x^2 m_y^6 + m_y^8) m_z^2\nonumber \\&+ 2 (5 m_x^6 + 3 m_x^4 m_y^2 + 9 m_x^2 m_y^4 + 5 m_y^6) m_z^4 \nonumber \\&- 2 (5 m_x^4 + 9 m_x^2 m_y^2 + 8 m_y^4) m_z^6 + (5 m_x^2 + 8 m_y^2) m_z^8 - m_z^{10}\,,\nonumber \\ \mu _{xyz}\equiv&\, m_x^8 - 4 m_x^6 m_y^2 + 6 m_x^4 m_y^4 - 4 m_x^2 m_y^6 + m_y^8 - 2 m_x^6 m_z^2\nonumber \\ \end{aligned}$$

$$\begin{aligned}&- 6 m_x^4 m_y^2 m_z^2 + 6 m_x^2 m_y^4 m_z^2 + 2 m_y^6 m_z^2\nonumber \\&- 4 m_x^4 m_z^4 - 12 m_x^2 m_y^2 m_z^4 - 12 m_y^4 m_z^4 + 10 m_x^2 m_z^6\nonumber \\&+ 14 m_y^2 m_z^6 - 5 m_z^8\,,\nonumber \\ \nu _{xyz}\equiv&\, ((m_x^2 - m_y^2)^3 + (-3 m_x^4 - 10 m_x^2 m_y^2 + m_y^4) m_z^2\nonumber \\&+ (3 m_x^2 + m_y^2) m_z^4 - m_z^6)\zeta _{xyz}/m_x^2\nonumber \\ =&\,\frac{(\chi _{xyz}-8m_x^2m_y^2m_z^2)\zeta _{xyz}}{m_x^2}\,,\nonumber \\ \rho _{xyz}\equiv&\, m_x^8 - 4 m_x^6 m_y^2 + 4 m_x^4 m_y^4 - m_y^8 - 4 m_x^6 m_z^2 \nonumber \\&- 4 m_x^4 m_y^2 m_z^2 + 4 m_y^6 m_z^2 + 4 m_x^4 m_z^4 - 6 m_y^4 m_z^4 \nonumber \\&+ 4 m_y^2 m_z^6 - m_z^8\,,\nonumber \\ \tau _{xyz}\equiv&\, m_x^8 - 4 m_x^6 m_y^2 + 6 m_x^4 m_y^4 - 4 m_x^2 m_y^6 + m_y^8\nonumber \\&- 4 m_x^6 m_z^2 - 8 m_x^4 m_y^2 m_z^2 - 8 m_x^2 m_y^4 m_z^2\nonumber \\&- 4 m_y^6 m_z^2 + 6 m_x^4 m_z^4 - 8 m_x^2 m_y^2 m_z^4 + 6 m_y^4 m_z^4\nonumber \\&- 4 m_x^2 m_z^6 - 4 m_y^2 m_z^6 + m_z^8\,,\nonumber \\ \varPhi _{xyz}\equiv&\, \big \{(m_x^2 - m_y^2)^5 - (m_x^2 - m_y^2) (7 m_x^6 + 11 m_x^4 m_y^2\nonumber \\&+ m_x^2 m_y^4 + 5 m_y^6) m_z^2 \nonumber \\&\quad +2 (3 m_x^6 + m_x^4 m_y^2 + 5 m_x^2 m_y^4 - 5 m_y^6) m_z^4\nonumber \\&+ 2 (3 m_x^4 - 2 m_x^2 m_y^2 + 5 m_y^4) m_z^6\nonumber \\&- (7 m_x^2 + 5 m_y^2) m_z^8 + m_z^{10}\big \}/m_x^2\,,\nonumber \\ \end{aligned}$$

$$\begin{aligned} \phi _{xyz}\equiv&\, m_x^2 m_y^2(m_x^2 - m_y^2)^2+m_x^2 m_z^2(m_x^2 - m_z^2)^2\nonumber \\&+m_y^2 m_z^2(m_y^2 - m_z^2)^2\,,\nonumber \\ \zeta _{xyz}\equiv&\, m_x^2(m_y^2+m_z^2)-(m_y^2-m_z^2)^2\,,\nonumber \\ r_{xyz}\equiv&\, m_x^2-m_y^2-m_z^2\,,\nonumber \\ t_{xyz}\equiv&\, m_x^2+m_y^2+m_z^2\,,\nonumber \\ \xi _{xyz}\equiv&\, m_x^6+m_x^4m_y^2-5m_x^2m_y^4+3m_y^6-3m_x^4m_z^2\nonumber \\&+4m_x^2m_y^2m_z^2+3m_y^4 m_z^2 +3m_x^2m_z^4\nonumber \\&-5 m_y^2m_z^4-m_z^6\,,\nonumber \\ a_{xyz}\equiv&\, m_x^8 - m_x^6 m_y^2 - 3 m_x^4 m_y^4 + 5 m_x^2 m_y^6 - 2 m_y^8\nonumber \\&- m_x^6 m_z^2+ 16 m_x^4 m_y^2 m_z^2 - 11 m_x^2 m_y^4 m_z^2 \nonumber \\&\qquad - 4 m_y^6 m_z^2 - 3 m_x^4 m_z^4 - 11 m_x^2 m_y^2 m_z^4\nonumber \\&+ 12 m_y^4 m_z^4 + 5 m_x^2 m_z^6 - 4 m_y^2 m_z^6 - 2 m_z^8\,,\nonumber \\ e_{xyz}\equiv&\, -3 m_x^{10} + 19 m_x^8 (m_y^2 + m_z^2) + (m_y^2 - m_z^2)^6\nonumber \\&\quad (m_y^2 + m_z^2)/m_x^4 - (m_y^2 - m_z^2)^4\nonumber \\&\quad (9 m_y^4 + 26 m_y^2 m_z^2 + 9 m_z^4)/m_x^2\nonumber \\&- m_x^6 (51 m_y^4 + 38 m_y^2 m_z^2 + 51 m_z^4) \nonumber \\&+ m_x^4 (m_y^2 + m_z^2) (75 m_y^4 + 26 m_y^2 m_z^2 + 75 m_z^4) \nonumber \\&+ (m_y^2 + m_z^2) (33 m_y^8 + 44 m_y^6 m_z^2\nonumber \\&+ 86 m_y^4 m_z^4 + 44 m_y^2 m_z^6 + 33 m_z^8)\nonumber \\&- m_x^2 (65 m_y^8 + 164 m_y^6 m_z^2 + 134 m_y^4 m_z^4\nonumber \\&+ 164 m_y^2 m_z^6 + 65 m_z^8)\,,\nonumber \\ \end{aligned}$$

$$\begin{aligned} f_{xyz}\equiv&\, 3 m_x^{10} - (m_y^2 - m_z^2)^7/m_x^4 - 19 m_x^8 (m_y^2 + m_z^2)\nonumber \\&+ (m_y^2 - m_z^2)^5 (9 m_y^2 + 13 m_z^2)/m_x^2\nonumber \\&+ m_x^6 (51 m_y^4 + 13 m_z^4) \nonumber \\&+ m_x^4 (-75 m_y^6 + 29 m_y^4 m_z^2 - 17 m_y^2 m_z^4 + 55 m_z^6)\nonumber \\&+ m_x^2 (65 m_y^8 - 22 m_y^4 m_z^4 + 8 m_y^2 m_z^6 - 99 m_z^8) \nonumber \\&+ (-33 m_y^{10} + 15 m_y^8 m_z^2 - 74 m_y^6 m_z^4\nonumber \\&+ 54 m_y^4 m_z^6 - 21 m_y^2 m_z^8 + 59 m_z^{10})\,,\nonumber \\ g_{xyz}\equiv&\, (m_x^2 - m_y^2)^6 (m_x^2 + m_y^2) - (m_x^2 - m_y^2)^2 (5 m_x^8\nonumber \\&+ 8 m_x^6 m_y^2 + 70 m_x^4 m_y^4 + 8 m_x^2 m_y^6 + 5 m_y^8) m_z^2\nonumber \\&+ (m_x^2 + m_y^2) (9 m_x^8 - 68 m_x^6 m_y^2 + 70 m_x^4 m_y^4\nonumber \\&- 68 m_x^2 m_y^6 + 9 m_y^8) m_z^4 \nonumber \\&+ (-5 m_x^8 + 124 m_x^6 m_y^2 + 2 m_x^4 m_y^4\nonumber \\&+ 124 m_x^2 m_y^6 - 5 m_y^8) m_z^6 \nonumber \\&- (m_x^2 + m_y^2) (5 m_x^4 + 54 m_x^2 m_y^2 + 5 m_y^4) m_z^8\nonumber \\&+ (9 m_x^4 + 2 m_x^2 m_y^2 + 9 m_y^4) m_z^{10}\nonumber \\&- 5 (m_x^2 + m_y^2) m_z^{12} + m_z^{14}\,,\nonumber \\ h_{xyz}\equiv&\, 2 m_x^{10} - (m_y^2 - m_z^2)^6/m_x^2 - 19 m_x^8 (m_y^2 + m_z^2)\nonumber \\&+ 4 m_x^2 (m_y^2 + m_z^2)^2 (m_y^4 + m_y^2 m_z^2 + m_z^4) \nonumber \\&+ (m_y^2 - m_z^2)^2 (m_y^2 + m_z^2) (5 m_y^4 + 14 m_y^2 m_z^2 + 5 m_z^4)\nonumber \\&- 2 m_x^4 (m_y^2 + m_z^2) (17 m_y^4 - 6 m_y^2 m_z^2 + 17 m_z^4)\nonumber \\&+ m_x^6 (43 m_y^4 + 14 m_y^2 m_z^2 + 43 m_z^4)\,,\nonumber \\ n_{xyz}\equiv&\, \big \{(m_x^2 - m_y^2)^5 - (m_x^2 - m_y^2)^2 (5 m_x^4 - 8 m_x^2 m_y^2\nonumber \\&+ 27 m_y^4) m_z^2 + 2 (5 m_x^6 - 24 m_x^4 m_y^2 - 7 m_x^2 m_y^4 + 14 m_y^6) m_z^4\nonumber \\&+ 2 (-5 m_x^4 + 31 m_x^2 m_y^2 + 14 m_y^4) m_z^6\nonumber \\ \end{aligned}$$

$$\begin{aligned}&+ (5 m_x^2 - 27 m_y^2) m_z^8 - m_z^{10}\big \}/m_x^2\,,\nonumber \\ p_{xyz}\equiv&\, \big \{(m_x^2 - m_y^2)^5 (m_x^4 - 3 m_x^2 m_y^2 + 3 m_y^4) \nonumber \\&- (m_x^2 - m_y^2)^2 (7 m_x^8 - 17 m_x^6 m_y^2 - 8 m_x^4 m_y^4\nonumber \\&+ 41 m_x^2 m_y^6 + m_y^8) m_z^2\nonumber \\&+ (21 m_x^{10} - 35 m_x^8 m_y^2 - 30 m_x^6 m_y^4 - 56 m_x^4 m_y^6\nonumber \\&+ 37 m_x^2 m_y^8 + 39 m_y^{10}) m_z^4\nonumber \\&- (35 m_x^8 + 10 m_x^6 m_y^2 - 8 m_x^4 m_y^4 + 52 m_x^2 m_y^6\nonumber \\&+ 87 m_y^8) m_z^6 + (35 m_x^6 + 50 m_x^4 m_y^2 + 66 m_x^2 m_y^4 + 83 m_y^6) m_z^8\nonumber \\&- (21 m_x^4 + 37 m_x^2 m_y^2 + 39 m_y^4) m_z^{10} \nonumber \\&+ (7 m_x^2 + 9 m_y^2) m_z^{12} - m_z^{14}\big \}/m_y^6\,,\nonumber \\ q_{xyz}\equiv&\, \big \{(m_x^2 - m_y^2)^5 - 8 m_y^2 (m_x^6 - 3 m_x^2 m_y^4 + 2 m_y^6) m_z^2\nonumber \\&- 4 (4 m_x^6 + 4 m_x^4 m_y^2 + 9 m_x^2 m_y^4 - 11 m_y^6) m_z^4\nonumber \\&+ 2 (19 m_x^4 + 20 m_x^2 m_y^2 - 13 m_y^4) m_z^6\nonumber \\&- 11 (3 m_x^2 + m_y^2) m_z^8 + 10 m_z^{10}\big \}/m_x^2\,,\nonumber \\ u_{xyz}\equiv&\, \big \{-m_x^{10} m_z^2 - (m_y^2 - m_z^2)^5 (2 m_y^4\nonumber \\&+ 5 m_y^2 m_z^2 - m_z^4)/m_x^2 + 2 m_x^8 (m_y^4 + 6 m_y^2 m_z^2 + 3 m_z^4)\nonumber \\&+ 2 (m_y^2 - m_z^2)^3 (5 m_y^6 - 19 m_y^4 m_z^2 + 5 m_y^2 m_z^4 - 3 m_z^6)\nonumber \\&- m_x^6 (10 m_y^6 + 25 m_y^4 m_z^2 + 38 m_y^2 m_z^4 + 15 m_z^6)\nonumber \\&+ 4 m_x^4 (5 m_y^8 - 10 m_y^6 m_z^2 - 14 m_y^4 m_z^4 + 8 m_y^2 m_z^6 + 5 m_z^8)\nonumber \\&+ m_x^2 (-20 m_y^{10} + 117 m_y^8 m_z^2 - 120 m_y^6 m_z^4\nonumber \\&+ 26 m_y^4 m_z^6 + 12 m_y^2 m_z^8 - 15 m_z^{10})\big \} /m_y^4\,,\nonumber \\ v_{xyz}\equiv&\, \big \{m_x^{12} - 2 m_x^{10} (5 m_y^2 + 3 m_z^2) - (m_y^2 - m_z^2)^4\nonumber \\&\quad (17 m_y^4 + 8 m_y^2 m_z^2 - m_z^4) \nonumber \\&+ m_x^8 (25 m_y^4 + 28 m_y^2 m_z^2 + 15 m_z^4)\nonumber \\&+ 2 m_x^2 (m_y^2 - m_z^2)^2 (25 m_y^6 - 23 m_y^4 m_z^2\nonumber \\&+ 13 m_y^2 m_z^4 - 3 m_z^6) - 4 m_x^6 (2 m_y^6 - 2 m_y^4 m_z^2\nonumber \\&+ 3 m_y^2 m_z^4 + 5 m_z^6)\nonumber \\ \end{aligned}$$

$$\begin{aligned}&+ m_x^4 (-41 m_y^8 + 56 m_y^6 m_z^2 + 50 m_y^4 m_z^4\nonumber \\&- 32 m_y^2 m_z^6 + 15 m_z^8)\big \}/m_y^4\,,\nonumber \\ w_{xyz}\equiv&\, m_x^8 (m_y^2 + m_z^2) - 5 m_x^6 (m_y^4 + m_z^4) - (m_y^2 - m_z^2)^4\nonumber \\&\quad (m_y^4 + 10 m_y^2 m_z^2 + m_z^4)/m_x^2\nonumber \\&+ 2 m_x^4 (5 m_y^6 - 12 m_y^4 m_z^2 - 12 m_y^2 m_z^4 + 5 m_z^6)\nonumber \\&+ (m_y^2 - m_z^2)^2 (5 m_y^6 + m_y^4 m_z^2 + m_y^2 m_z^4 + 5 m_z^6)\nonumber \\&- 2 m_x^2 (5 m_y^8 - 19 m_y^6 m_z^2 + 16 m_y^4 m_z^4 - 19 m_y^2 m_z^6 + 5 m_z^8)\,,\nonumber \\ A_{xyz}\equiv&\, m_x^8 - 4 m_x^6 (m_y^2 + m_z^2) + 2 (m_y^2 - m_z^2)^4\nonumber \\&(m_y^2 + m_z^2)/m_x^2 + 4 m_x^4 (m_y^4 + m_z^4)\nonumber \\&- (m_y^2 - m_z^2)^2 (5 m_y^4 - 6 m_y^2 m_z^2 + 5 m_z^4)\nonumber \\&+ 2 m_x^2 (m_y^6 - 3 m_y^4 m_z^2 - 3 m_y^2 m_z^4 + m_z^6)\,,\nonumber \\ B_{xyz}\equiv&\, \big \{m_x^{10} + 4 m_x^8 (m_y^2 - m_z^2) + 4 m_z^2 (-7 m_y^2 + m_z^2)\nonumber \\&(-m_y^2 + m_z^2)^3 + (m_y^2 - m_z^2)^5 (5 m_y^2 + m_z^2)/m_x^2\nonumber \\&+ m_x^6 (-33 m_y^4 - 20 m_y^2 m_z^2 + 5 m_z^4)\nonumber \\&+ 4 m_x^4 (15 m_y^6 - 10 m_y^4 m_z^2 - 3 m_y^2 m_z^4)\nonumber \\&+ m_x^2 (-37 m_y^8 + 60 m_y^6 m_z^2 - 86 m_y^4 m_z^4\nonumber \\&+ 68 m_y^2 m_z^6 - 5 m_z^8)\big \}/m_y^2\,. \end{aligned}$$

(C.1)

1.2 Limits for degenarate masses

We find when all mass arguments are equal

$$\begin{aligned} \varDelta _{mmm}&=-3m^4,&\omega _{mmm}&=0,&\varXi _{mmm}&=3m^6,\nonumber \\ \chi _{mmm}&=-m^6,&\varTheta _{mmm}&=-15m^8,&\theta _{mmm}&=12m^{10},\nonumber \\ \mu _{mmm}&=-9m^8,&\nu _{mmm}&=-18m^8,&\rho _{mmm}&=-3m^8,\nonumber \\ \tau _{mmm}&=-27m^8,&\varPhi _{mmm}&=9m^8,&\phi _{mmm}&=0,\nonumber \\ \zeta _{mmm}&=2m^4,&r_{mmm}&=-m^2,&t_{mmm}&=3m^2\nonumber \\ \xi _{mmm}&=m^6,&a_{mmm}&=-3m^8,&e_{mmm}&=135m^{10},\nonumber \\ f_{mmm}&=-27m^{10},&g_{mmm}&=27m^{14},&h_{mmm}&=0,\nonumber \\ n_{mmm}&=33m^8,&p_{mmm}&=-48m^8,&q_{mmm}&=-6m^8,\nonumber \\ u_{mmm}&=-93m^8,&v_{mmm}&=69m^8,&w_{mmm}&=-12m^{10},\nonumber \\ A_{mmm}&=-7m^8,&B_{mmm}&=-39m^8. \end{aligned}$$

(C.2)

When there are two different masses, we have

$$\begin{aligned} \varDelta _{m_1m_1m_2}&=m_2^2(m_2^2-4m_1^4)=\varDelta _{m_1m_2m_1}=\varDelta _{m_2m_1m_1},\nonumber \\ \omega _{m_1m_1m_2}&=m_2^2(m_2^2-m_1^2)=\omega _{m_2m_1m_1},\nonumber \\ \omega _{m_1m_2m_1}&=2m_2^2(m_1^2-m_2^2),\nonumber \\ \varXi _{m_1m_1m_2}&=m_2^2(4m_1^4-2m_1^2m_2^2+m_2^4)\nonumber \\&=\varXi _{m_1m_2m_1}=\varXi _{m_2m_1m_1},\nonumber \\ \chi _{m_1m_1m_2}&=-m_2^2(m_2^2-2m_1^2)^2=\chi _{m_1m_2m_1},\nonumber \\ \chi _{m_2m_1m_1}&=4 m_1^4 m_2^2 - 6 m_1^2 m_2^4 + m_2^6\nonumber \\ \varTheta _{m_1m_1m_2}&=-20 m_1^4 m_2^4 + 4 m_1^2 m_2^6 + m_2^8=\varTheta _{m_1m_2m_1},\nonumber \\ \varTheta _{m_2m_1m_1}&=-24 m_1^6 m_2^2 + 16 m_1^4 m_2^4 - 8 m_1^2 m_2^6 + m_2^8,\nonumber \\ \theta _{m_1m_1m_2}&=-m_2^4 (-44 m_1^6 + 44 m_1^4 m_2^2 - 13 m_1^2 m_2^4 + m_2^6),\nonumber \\ \theta _{m_1m_2m_1}&=-2 m_2^4 (-12 m_1^6 + 10 m_1^4 m_2^2 - 5 m_1^2 m_2^4 + m_2^6),\nonumber \\ \theta _{m_2m_1m_1}&=-12 m_1^6 m_2^4 + 34 m_1^4 m_2^6 - 11 m_1^2 m_2^8 + m_2^{10},\nonumber \\ \mu _{m_1m_1m_2}&=-28 m_1^4 m_2^4 + 24 m_1^2 m_2^6 - 5 m_2^8,\nonumber \\ \mu _{m_1m_2m_1}&=-8 m_1^6 m_2^2 - 2 m_1^2 m_2^6 + m_2^8,\nonumber \\ \mu _{m_2m_1m_1}&=-4 m_1^4 m_2^4 - 6 m_1^2 m_2^6 + m_2^8,\nonumber \\ \nu _{m_1m_1m_2}&=-36 m_1^4 m_2^4 + 24 m_1^2 m_2^6 - 7 m_2^8\nonumber \\&+\frac{m_2^{10}}{m_1^2}=\nu _{m_1m_2m_1},\nonumber \\ \nu _{m_2m_1m_1}&=2m_1^2 m_2^2 (-4 m_1^4 - 6 m_1^2 m_2^2 + m_2^4),\nonumber \\ \rho _{m_1m_1m_2}&=-m_2^2 (4 m_1^6 + 2 m_1^4 m_2^2 - 4 m_1^2 m_2^4 + m_2^6)\nonumber \\ \end{aligned}$$

$$\begin{aligned}&=\rho _{m_1m_2m_1}, \nonumber \\ \rho _{m_2m_1m_1}&=4m_1^4 m_2^4 - 8 m_1^2 m_2^6 + m_2^8,\nonumber \\ \tau _{m_1m_1m_2}&=-24 m_1^6 m_2^2 + 4 m_1^4 m_2^4 - 8 m_1^2 m_2^6 + m_2^8\nonumber \\&=\tau _{m_1m_2m_1}=\tau _{m_2m_1m_1},\nonumber \\ \varPhi _{m_1m_1m_2}&=8 m_1^4 m_2^4 + 12 m_1^2 m_2^6 - 12 m_2^8 + \frac{m_2^{10}}{m_1^2},\nonumber \\ \varPhi _{m_1m_2m_1}&=-16 m_1^6 m_2^2 + 40 m_1^4 m_2^4 - 24 m_1^2 m_2^6\nonumber \\&\quad + 10 m_2^8 - \frac{m_2^{10}}{m_1^2},\nonumber \\ \varPhi _{m_2m_1m_1}&=8 m_1^6 m_2^2 + 12 m_1^4 m_2^4 - 12 m_1^2 m_2^6 + m_2^8,\nonumber \\ \phi _{m_1m_1m_2}&=2 m_1^2 m_2^2 (m_1^2 - m_2^2)^2=\phi _{m_1m_2m_1} =\phi _{m_2m_1m_1},\nonumber \\ \zeta _{m_1m_1m_2}&=3 m_1^2 m_2^2 - m_2^4=\zeta _{m_1m_2m_1},\nonumber \\ \zeta _{m_2m_1m_1}&=2m_1^2m_2^2\nonumber \\ r_{m_1m_1m_2}&=-m_2^2=r_{m_1m_2m_1},\qquad r_{m_2m_1m_1}=-2 m_1^2 + m_2^2,\nonumber \\ t_{m_1m_1m_2}&=2m_1^2+m_2^2=t_{m_1m_2m_1}=t_{m_2m_1m_1},\nonumber \\ \xi _{m_1m_1m_2}&=4 m_1^4 m_2^2 - 2 m_1^2 m_2^4 - m_2^6,\nonumber \\ \xi _{m_1m_2m_1}&=- 2 m_1^2 m_2^4 + 3 m_2^6,\nonumber \\ \xi _{m_2m_1m_1}&=2 m_1^4 m_2^2 - 2 m_1^2 m_2^4 + m_2^6,\nonumber \\ a_{m_1m_1m_2}&=-2 m_1^4 m_2^4 + m_1^2 m_2^6 - 2 m_2^8=a_{m_1m_2m_1},\nonumber \\ a_{m_2m_1m_1}&=-12 m_1^6 m_2^2 + 10 m_1^4 m_2^4 - 2 m_1^2 m_2^6 + m_2^8,\nonumber \\ e_{m_1m_1m_2}&=m_2^4 \bigg (96 m_1^6 - 48 m_1^4 m_2^2 + 48 m_1^2 m_2^4 + 52 m_2^6\nonumber \\&\quad - \frac{14 m_2^8}{m_1^2} + \frac{m_2^{10}}{m_1^4}\bigg )=e_{m_1m_2m_1},\nonumber \\ \end{aligned}$$

$$\begin{aligned} e_{m_2m_1m_1}&=480 m_{1}^{10} - 592 m_{1}^8 m_2^2 + 352 m_{1}^6 m_2^4\nonumber \\&\quad - 140 m_{1}^4 m_2^6 + 38 m_{1}^2 m_2^8 - 3 m_2^{10},\nonumber \\ f_{m_1m_1m_2}&=m_2^4 \bigg (-96 m_1^6 + 192 m_1^4 m_2^2 - 240 m_1^2 m_2^4\nonumber \\&\quad + 136 m_2^6 - \frac{20 m_2^8}{m_1^2} + \frac{m_2^{10}}{m_1^4}\bigg ),\nonumber \\ f_{m_1m_2m_1}&=m_2^4 \bigg (48 m_1^6 - 144 m_1^4 m_2^2 + 140 m_1^2 m_2^4\nonumber \\&\quad - 86 m_2^6 + \frac{16 m_2^8}{m_1^2} - \frac{m_2^{10}}{m_1^4}\bigg ),\nonumber \\ f_{m_2m_1m_1}&=-48 m_1^8 m_2^2 - 8 m_1^6 m_2^4 + 64 m_1^4 m_2^6\nonumber \\&\quad - 38 m_1^2 m_2^8 + 3 m_2^{10},\nonumber \\ g_{m_1m_1m_2}&=m_2^4 \bigg (-96 m_1^{10} + 240 m_1^8 m_2^2 - 128 m_1^6 m_2^4\nonumber \\&\quad + 20 m_1^4 m_2^6 - 10 m_1^2 m_2^8 + m_2^{10}\bigg )\nonumber \\&=g_{m_1m_2m_1}=g_{m_2m_1m_1},\nonumber \\ h_{m_1m_1m_2}&=m_2^4 \bigg (8 m_1^6 - 16 m_1^4 m_2^2 - 2 m_1^2 m_2^4\nonumber \\&\quad + 11 m_2^6 - \frac{m_2^8}{m_1^2}\bigg )=h_{m_1m_2m_1},\nonumber \\ h_{m_2m_1m_1}&=2 m_2^2 (24 m_1^8 - 56 m_1^6 m_2^2 + 50 m_1^4 m_2^4\nonumber \\&- 19 m_1^2 m_2^6 + m_2^8),\nonumber \\ n_{m_1m_1m_2}&=-24 m_1^4 m_2^4 + 80 m_1^2 m_2^6 - 22 m_2^8\nonumber \\&\quad - \frac{m_2^{10}}{m_1^2}=n_{m_1m_2m_1},\nonumber \\ n_{m_2m_1m_1}&=(10 m_1^4 - 8 m_1^2 m_2^2 + m_2^4) \nonumber \\&\quad (12 m_1^4 - 2 m_1^2 m_2^2 + m_2^4),\nonumber \\ \end{aligned}$$

$$\begin{aligned} p_{m_1m_1m_2}&=-m_2^4\big \{ 24 m_1^{10} + 176 m_1^8 m_2^2 - 234 m_1^6 m_2^4\nonumber \\&+ 97 m_1^4 m_2^6 - 16 m_1^2 m_2^8 + m_2^{10}\big \}/m_1^6,\nonumber \\ p_{m_1m_2m_1}&=-120 m_1^8 + 104 m_1^6 m_2^2 - 46 m_1^4 m_2^4\nonumber \\&+ 17 m_1^2 m_2^6 - 3 m_2^8,\nonumber \\ p_{m_2m_1m_1}&=m_2^4 \bigg (24 m_1^4 + 20 m_1^2 m_2^2 - 158 m_2^4\nonumber \\&\quad + \frac{80 m_2^6}{m_1^2} - \frac{15 m_2^8}{m_1^4} + \frac{m_2^{10}}{m_1^6}\bigg ),\nonumber \\ q_{m_1m_1m_2}&=-24 m_1^4 m_2^4 + 52 m_1^2 m_2^6 - 44 m_2^8 + \frac{10 m_2^{10}}{m_1^2},\nonumber \\ q_{m_1m_2m_1}&=-52 m_1^4 m_2^4 + 58 m_1^2 m_2^6 - 11 m_2^8 - \frac{m_2^{10}}{m_1^2},\nonumber \\ q_{m_2m_1m_1}&=12 m_1^6 m_2^2 - 14 m_1^4 m_2^4 - 5 m_1^2 m_2^6 + m_2^8,\nonumber \\ u_{m_1m_1m_2}&=-m_2^4 \big \{48 m_1^{10} + 152 m_1^8 m_2^2 - 168 m_1^6 m_2^4\nonumber \\&\quad + 76 m_1^4 m_2^6 - 16 m_1^2 m_2^8 + m_2^{10}\big \}/m_1^6,\nonumber \\ u_{m_1m_2m_1}&=-280 m_1^6 m_2^2 + 256 m_1^4 m_2^4 - 82 m_1^2 m_2^6\nonumber \\&\quad + 15 m_2^8 - \frac{2 m_2^{10}}{m_1^2},\nonumber \\ u_{m_2m_1m_1}&=-24 m_1^4 m_2^4 - 88 m_1^2 m_2^6 + 20 m_2^8 - \frac{m_2^{10}}{m_1^2},\nonumber \\ \end{aligned}$$

$$\begin{aligned} v_{m_1m_1m_2}&=m_2^4 (38 m_1^4 - 16 m_1^2 m_2^2 + m_2^4)\nonumber \\&\quad (4 m_1^4 - 2 m_1^2 m_2^2 + m_2^4)/m_1^4,\nonumber \\ v_{m_1m_2m_1}&=232 m_1^6 m_2^2 - 256 m_1^4 m_2^4 + 110 m_1^2 m_2^6 - 17 m_2^8,\nonumber \\ v_{m_2m_1m_1}&=m_2^4 \bigg (48 m_1^4 - 32 m_1^2 m_2^2 + 68 m_2^4 - \frac{16 m_2^6}{m_1^2} + \frac{m_2^8}{m_1^4}\bigg ),\nonumber \\ w_{m_1m_1m_2}&=-m_2^4 \big \{24 m_1^8 - 14 m_1^4 m_2^4 + m_1^2 m_2^6 + m_2^8\big \}/\nonumber \\&m_1^2=w_{m_1m_2m_1},\nonumber \\ w_{m_2m_1m_1}&=2 m_1^2 m_2^2 (12 m_1^6 - 14 m_1^4 m_2^2 - 5 m_1^2 m_2^4 + m_2^6),\nonumber \\ A_{m_1m_1m_2}&=-20 m_1^4 m_2^4 + 22 m_1^2 m_2^6 - 11 m_2^8 + \frac{2 m_2^{10}}{m_1^2}=A_{m_1m_2m_1},\nonumber \\ A_{m_2m_1m_1}&=-8 m_1^6 m_2^2 + 8 m_1^4 m_2^4 - 8 m_1^2 m_2^6 + m_2^8,\nonumber \\ B_{m_1m_1m_2}&=m_2^4 \bigg (-136 m_1^4 + 124 m_1^2 m_2^2 - 30 m_2^4 + \frac{4 m_2^6}{m_1^2} - \frac{m_2^8}{m_1^4}\bigg ),\nonumber \\ B_{m_1m_2m_1}&=-48 m_1^6 m_2^2 - 8 m_1^4 m_2^4 + 36 m_1^2 m_2^6 - 24 m_2^8 + \frac{5 m_2^{10}}{m_1^2},\nonumber \\ B_{m_2m_1m_1}&=m_2^4 \{8 m_1^6 - 48 m_1^4 m_2^2 + m_2^6\}/m_1^2. \end{aligned}$$

(C.3)