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

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{\text{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{\text{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.


I. 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 Z 2 symmetry [6-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 λ hhh and its quartic coupling λ hhhh . 1 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 selfcouplings 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, λ hhh must deviate from its SM value by at least 20 − 30%. λ 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.
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 λ hhh , some limits are already available: using single Higgs production data from LHC Run 2, the ATLAS collaboration set a limit on the ratio as −3.2 < κ λ < 11.9 at 95% confidence level (CL) [36], while with double Higgs production, the best intervals obtained at 95% CL are respectively −5.0 < κ λ < 12.1 from ATLAS [37] (see also Ref. [38]) and −11 < κ λ < 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 κ λ 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 highluminosity upgrade of the LHC (the HL-LHC) could reach 0.5 < κ λ < 1.6 (at 68% CL) with an integrated luminosity of 3 ab −1 [42] (see also Ref. [43]). A possible high-energy version of the LHC (the HE-LHC), with centre-of-mass energy of √ s = 27 TeV, was found in Ref. [44] to be able to attain 0.54 < λ hhh /λ SM hhh < 1.46 (at 68% CL) using 15 ab −1 of data. 2 Turning next to the case of possible future lepton colliders, the initial stage of the International Linear Collider (ILC) running at √ s = 250 GeV cannot access λ 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 ab −1 of data [41]. With the data from further ILC extensions to 500 GeV (4 ab −1 ) or even 1 TeV (8 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 ab −1 of data at 380 GeV, 2.5 ab −1 at 1.5 TeV, and 5 ab −1 at 3 TeV, obtain a final accuracy of 0.93 < κ λ < 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 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 ab −1 of data, the FCC-hh would only be able to constrain the ratio λ exp. hhhh /λ SM hhhh to be approximately between −4 and 16 [52] (at 95% confidence level). On the theoretical side, the first one-loop calculations of λ 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 λ 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 2 Note that a different analysis was performed in Ref. [45], including fewer sources of background, and found a possible accuracy of 15% on the Higgs trilinear coupling, still at 68% CL. Both analyses, in Refs. [44] and [45], derive limits using double-Higgs production and the bbγγ channel, however expected improvements for the other available channels -bbbb, bbτ τ and bbV V -should better the accuracy reached for λ hhh down to 10 − 20% [42,44]. radiative corrections involving additional BSM scalars can, in the non-decoupling regime, cause a significant enhancement of λ 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 λ 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 λ hhh to a similar level of accuracy. 3 In Refs. [83,84], the leading O(α s α t ) corrections to λ 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 λ 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 λ 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 λ hhh , we also give new formulae for the two-loop corrections to λ 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 Sec. II, before describing the set-up of our calculations in Sec. III. In Sec. IV, we give general results for derivatives of the effective potential, expressed in the MS scheme. Then we present our analytical results for the SM and three BSM scenarios, in both MS and on-shell schemes, in Sec.

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

A. 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 Z 2 symmetry under which the scalar doublets transform respectively as Φ 1 → Φ 1 and Φ 2 → −Φ 2 , and the different families of fermions have charges ±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. 4 The Z 2 symmetry can be broken softly -i.e. without reintroducing dangerous FCNCs -via an off-diagonal mass term m 2 3 Φ † 2 Φ 1 + h.c. . We follow the conventions of Ref. [12] and write the tree-level scalar potential as 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] We expand each of the scalar doublets as [12] Φ Here v 1 and v 2 denote the VEVs of the neutral components of the scalar doublets, and satisfy the We will further assume that the values of the parameters in the potential ensure that both v 1 = 0 and v 2 = 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 2 1 and m 2 2 -using the minimisation conditions of the potential (i.e. the tadpole equations), which read at tree level The angle β in the above two equations is defined from the ratio of VEVs v 2 /v 1 ≡ tan β, and we make use of the following shorthand notations λ 345 ≡ λ 3 +λ 4 +λ 5 , c x ≡ cos x, and s x ≡ sin x. Once the tadpole equations have been applied, six free parameters remain in the 2HDM Higgs sector, namely m 2 3 , λ i (i = 2, · · · , 5), tan β , (II. 6) where we note that one of the quartic couplings -here we choose λ 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 Z 2 -breaking scale M defined as M 2 ≡ 2m 2 3 /s 2β . Moreover, while the angle β 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 β , with to the component fields of the two doublets by the angle β, we obtain new states as In this new basis -often refered to as the Higgs basis -only φ 1 carries a VEV, i.e. φ 1 = v and φ 2 = 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) (II.9) φ 1 and φ 2 are, however, not mass eigenstates in this basis, and their mass matrix reads (II.10) An additional rotation of angle (α − β) is necessary to obtain tree-level CP-even mass eigenstates, which will be denoted h and H The latter of the two relations shows that α is defined as the CP-even Higgs mixing angle. In turn the tree-level mass eigenvalues for h and H can be found as 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 β−α = 1 or c β−α = 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 and therefore we find that all the masses of the additional Higgs bosons Φ = H, A, H ± take the form m 2 Φ = M 2 +λv 2 , (II.14) whereλ denotes some simple function of Lagrangian quartic couplings (and tan β) -as given in equations (II.9) and (II.13). Moreover, in the alignment limit we can obtain the h-field dependent masses of the additional scalars with the simple replacement v → 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 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 ± , and the CP-even mixing angle α (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 We should emphasise here that as these are tree-level relations, they can only be used if the masses m Φ are tree-level 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) 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.

B. 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 Z 2 symmetry acts only on one of the two Higgs doublet -say Φ 2 to fix the notation -and remains unbroken after EWSB. This condition forbids the presence of a mass term (Φ † 2 Φ 1 + h.c.), as well as the appearance of a non-zero VEV for the neutral component of Φ 2 . In this case, the Z 2 -odd doublet Φ 2 cannot mix with the SM-like doublet Φ 1 , nor can it couple to the fermion sector.
We follow the conventions of Ref. [107] and we expand the two scalar doublets as where the notations for the component fields are common with the 2HDM. Because they do not couple to fermions, the components of Φ 2 -i.e. H, A, H ± -are referred to as inert scalars. The lightest of the two neutral of these 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 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 where all parameters are real. As only one of the two Higgs doublets acquires a VEV, we have a single tadpole equation which we use to eliminate the mass parameter µ 1 . We are then left with five free parameters in the Higgs sector, namely while v is related to the Fermi constant and λ 1 is constrained by m h = 125 GeV. We note that for the 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 equation (II.2).
We can obtain the tree-level, field-dependent, masses of the inert scalars as where λ H,A ≡ λ 3 + λ 4 ± λ 5 . It is interesting to note that λ 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 Sec. V, we will be considering a particular DM-inspired scenario. One can indeed distinguish [8,91]  -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 µ 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 µ 2 to be small to allow m H 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 Sec. V consider an IDM scenario where µ 2 = 0 so that m H m h /2 while the masses of A and H ± can be taken large by increasing the values of the quartic couplings.

C. 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 ϕ 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 Φ and the real singlet ϕ S as 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 Φ and ϕ S : 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 Z 2 symmetry under which S changes sign [7]. In this case, the potential reduces to Note however that if the sum µ 2 S + 1/2λ HS v 2 were to be negative, the Z 2 symmetry would be spontaneously broken by a singlet VEV, which also generates again the trilinear couplings κ 1 and κ 2 in equation (II.23). We will choose to consider such a Z 2 -symmetric HSM, ensuring that µ 2 S > −1/2λ HS v 2 , with the additional motivation 5 of avoiding mixing between the CP-even component of the Higgs doublet and the singlet.
For λ HS ≥ 0, the HSM tree-level potential is manifestly bounded from below provided that also λ H and λ S are positive. If however λ HS < 0, then one must impose the condition λ H λ S > 1/4λ 2 HS to avoid the appearance of unstable directions in the potential. Turning next to the counting of parameters in the Z 2 -symmetric HSM, one is left with three free parameters, namely Indeed, µ 2 and λ 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 λ HS for the (tree-level) singlet mass m 2 S , using the relation Finally, as in the 2HDM and the IDM, the Higgs-field-dependent singlet mass is obtained with the FIG. 1. Diagrammatic illustration of the Higgs three-point function Γ hhh .

A. Computation in the 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 Γ hhh and Γ 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 with Φ 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 ± 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 eff to successive orders in perturbation theory, up to two loops, as where κ is the loop factor, defined in eq. (A.1). In terms of V eff , we can define effective Higgs trilinear and quartic couplings, and their respective loop expansions, as Note that these definitions correspond to the following choices of normalisation Because we are considering scenarios without mixing, we can write the tree-level contributions to λ hhh and λ 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 We can then rewrite λ hhh as where for the second equality we used equation (III.5) and we define Similarly, we can write 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 (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. (III.3) are written in terms of MSrenormalised parameters.
Before discussing two-loop corrections, we should also review known results for the effectivepotential calculation of one-loop corrections to Higgs self-couplings. Using the well-known supertrace formula [122], the dominant one-loop contributions to V eff can be found to be, for the 2HDM, IDM, and HSM where the sum on heavy scalars Φ includes Φ = H, A, H ± for the 2HDM, Φ = A, H ± for the IDM, and Φ = S for the HSM. the BSM scalars, respectively, and n Φ = 1 for H and A, and n Φ = 2 for H ± . The notation log x is defined in equation (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 D 3 and D 4 , we obtain straightforwardly the leading one-loop corrections to the Higgs trilinear coupling as and for the Higgs quartic coupling, For both equations, we define the shorthand notation M to denote 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 where each index S or F indicates a scalar or a Dirac-fermion propagator. Analytic expressions for these, in the 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) F F S corresponds to the sum of V F F S and VFF 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.

B. Conversion from MS to on-shell renormalisation
Although it makes calculations simple, the 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 6 (or pole) masses and the physical Higgs VEV v phys = ( √ 2G 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, aŝ In these two equations, Z OS h and Z MS h are respectively the OS-and MS-scheme Higgs WFR constants, and their ratio can straightforwardly be computed in terms of the corresponding WFR where Π hh (p 2 ) is the finite part of the Higgs self-energy, evaluated at external momentum p.
(III. 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 Finally, the MS-and OS-renormalised versions of the Higgs VEV satisfy the equation This general prescription for the scheme conversion is significantly simplified in our case, in par-  h 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] δ (1) On the other hand, for the WFR, the additional scalars do give new, model-dependent, contributions -to which we will return in Sec. V.

IV. GENERAL MS EXPRESSIONS
One may notice from the discussion of models in Sec. II that in the scenarios without mixing that we consider, the field-dependent masses always take the form where µ 2 i andλ i have respectively mass-dimensions 2 and 0 -λ 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 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.

A. 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] where λ 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 Sec. III, we obtain in terms of the masses m 2 i and couplingsλ i : where the notation (1 ↔ 2) indicates the permutation of indices 1 and 2.

B. Sunrise diagrams
Next, we turn to the sunrise diagrams, corresponding to V 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 equation (A.9) -as well as products of A functions for V (2) F F S . Furthermore, almost all I functions come multiplied by (v + h) 2 in models without mixing (because of fielddependent couplings or masses), so we provide here results 7 for derivatives of the product (v + h) 2 I(m 2 1 (h), m 2 2 (h), m 2 3 (h)). 7 Expressions for derivatives of I(m 2 1 (h), m 2 2 (h), m 2 3 (h)) alone can then be obtained straightforwardly.
First, for the second derivative we obtain where (123) denotes the sum on cyclical permutations of the indices {1, 2, 3} and 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 For the fourth derivative we have

A. 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 where g 3 denotes the SU (3) C gauge coupling. The dominant contributions to the Higgs trilinear couplings, expressed in terms of MS parameters, are then given by which corresponds to equation (11) in Ref. [85]. For the quartic coupling we find where we have taken the same limit as for Π (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 In Fig. 3 The first BSM model that we consider is an aligned scenario of a 2HDM. As discussed in Sec. II A, requiring alignment -in other words fixing α = β − π/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 ± , andM -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 ± 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 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 The full expressions for the derivatives of the V F F S diagrams are quite long, so we here give only the leading results -in the third lines of eqs. (V.8) and (V.9) -and we provide the complete results in Appendix B. The first and second lines 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. (V.8) and (V.9) should tend to zero in the limit m Φ → ∞. Given that the scalar masses all satisfy a relation of the form m 2 Φ = M 2 +λv 2 , one could in principle achieve the previously-mentioned limit by taking either M 2 orλ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 → ∞, while keepingλ fixed. It is then straightforward to see that the above results for the corrections to λ hhh and λ hhhh decouple properly as all the terms involved are of the form, with m < n Furthermore, we must also include finite WFR for the Higgs bosons on the external legs, following equation (III.14). The parameter tan β 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 It is important to remember when converting the BSM scalar masses that the parameter points For this reason, one may at first think that there is no use to convert M , and that one may simply continue using its 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 MS parameters, with m 2 Φ = M 2 +λv 2 , when we take the limit M → ∞. If instead the corrections to the Higgs self-couplings are expressed in terms of pole (i.e. OS-renormalised) masses M Φ and of the soft-breaking scale M still in the MS scheme, we must use a one-loop relation between M Φ 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 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 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 denoteM , should be interpreted as the OS-renormalised version of the soft breaking scale of the 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 newM to its MS counterpart M with a finite counterterm denoted δ OS M 2 as 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 2 Φ =M 2 +λv 2 . With this prescription, we obtain up to one-loop order (V.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 (V.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 phys , (V. 15) and, phys .
(V. 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 equation (V.10) and therefore decouple for M 2 Φ =M 2 +λv 2 andM → ∞, as desired. Importantly, we should point out that, while we define our "on-shell" prescription forM 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.
where all masses are understood to be field-dependent masses. being either m A or m H ± -to the Higgs trilinear and quartic couplings: and The results in both equations (V.18) and (V. 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 ± running in the loops. The third term originates from the third sunrise diagram in Fig. 5, while the last three terms -proportional to λ 2 -arise from the eight-shaped diagrams.
b. 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 ± in the IDM as well as the momentum-derivative of the one-loop self-energy of the 125-GeV Higgs boson The conversion of the Higgs VEV is the same as in the SM, as given in eq. (III. 19). Moreover, as the coupling λ 2 only appears at two loops, we do not need any translation for it. Finally, had we not set µ 2 to zero, it would have appeared in the one-loop correction to the Higgs self-couplingsc.f. eqs. (III.10) and (III.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 (MS) relation µ 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 and for the Higgs quartic coupling It is interesting to note that because of WF and VEV renormalisation -which give the second line

D. A Higgs-Singlet Model with Z 2 symmetry
Finally, we turn to the Z 2 -symmetric HSM, introduced in Sec. II C 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 a. MS calculation -Following the same procedure as for the 2HDM and the IDM, we find in terms of MS parameters for the corrections to the Higgs trilinear coupling, and for those to the Higgs quartic coupling. In both equations, the first and second lines 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 ∂ 4 /∂h 4 of V (2) (instead of applying D 4 ), we reproduce the same result as in equation (29) of Ref. [35] for the two-loop corrections to the Higgs quartic coupling (after taking the limit of m h → 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].
b. 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 after neglecting the light scalar masses.
Moreover, as in the 2HDM, the mass parameter µ S appears in the one-loop corrections to the Higgs couplings. Therefore, we also need a finite "OS" counterterm -that we denote δ OS µ 2 S -for it, which we define asμ whereμ S and µ S are the OS-and MS-renormalised versions of the mass parameter. As in the 2HDM, we determine δ OS µ 2 S 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 2 S =μ 2 S + 1 2 λ S v 2 and taking the limitμ S → ∞ while keeping λ S fixed. We find eventually We use the corrections to the Higgs trilinear coupling to determine the finite counterterm δ OS µ 2 S , 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 Using all the above results, we obtain finally the following OS-renormalised expressions for the two-loop corrections to the Higgs self-couplings 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 .

VI. NUMERICAL EXAMPLES
We now turn to the discussion of the numerical behaviour of the BSM corrections computed in Sec. V. 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 Sec. II), we require that unitarity should not be violated. For this, we choose to take as our criterion 9 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.

A. Decoupling limit
A natural first point to study numerically is the decoupling behaviour of the two-loop BSM corrections. In Sec. V, we had discussed the decoupling of BSM effects in terms of analytical 9 Note that one could in principle argue that (one-)loop level perturbative unitarity conditions should be considered as we work at two loops. However, this would open a long discussion, which we prefer to leave for separate work. t β = 1.5 The left part of Fig. 7 illustrates the decoupling of the one-and two-loop corrections to the Higgs trilinear couplingλ hhh in the 2HDM, when expressed in terms of OS parameters. More precisely, the plot shows, as a function ofM , the BSM deviations δR, defined -in the alignment limit -at one and two loops respectively as for an example point with degenerate BSM scalar masses and tan β = 1.5. For each of the values of the difference M 2 Φ −M 2 =λv 2 that we consider -namely (200 GeV) 2 , (300 GeV) 2 , and (400 GeV) 2 -we can indeed observe that the BSM effects decouple rapidly for increasingM , both at the oneloop 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 β dependance at two loops and its impact on the decoupling of BSM corrections. Indeed, as can be seen from equations (V. 15) and (V.16), the two-loop contributions to the Higgs self-couplings involve cot 2 β and cot 2 2β, even in the alignment limit, because these appear in tree-level scalar couplings. While cot 2 β obviously vanishes very quickly with increasing tan β, cot 2 2β grows very fast -like tan 2 β for large tan β. This implies possible enhancements of the two-loop corrections for large values of tan β, only limited by the upper bound that the constraint of perturbative unitarity puts on tan β. We show on the right side of Fig. 7 the magnitude of the two-loop deviations -obtained as δR 2 − δR 1 -as a function ofM . We present results for four values of tan β, tan β = 1.75 being close to the maximal value allowed under the criterion of tree-level perturbative unitarity [81] forM = 0 and M Φ = 400 GeV.
As expected, while δR 1 does not depend on tan β, we can observe significant enhancements of δR 2 when increasing tan β, especially for small values ofM . WhenM grows, the effect of the terms in δ (2)λ hhh proportional to cot 2 2β diminishes because these have higher powers of the suppression factor than some of the other, non-tan β-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 β. when the BSM scalars obtain their masses mostly from the Higgs VEV and quartic couplings, and can therefore not be decoupled. 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 6 Φ dependence of part of the two-loop terms -see equations (V.15), (V.22), and (V.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 β for the 2HDM, λ 2 for the IDM, and λ 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 β = 1, λ 2 = 0, and λ S = 0, up to the maximal value for which the condition of tree-level unitarity is verified all the way to M Φ = 500 GeV -this gives the maximal values tan β = 1.4, λ 2 = 6.0, and λ S = 3.7. Turning first to the 2HDM, we observe that varying tan β in the allowed range does not produce any significant enhancement of the two-loop corrections; in other words the eight-shaped diagrams in V eff only have a limited impact in the 2HDM. Indeed requiring that tree-level unitarity is preserved up to M Φ = 500 GeV and forM = 0 strongly constrains the maximal possible value of tan β -we will return to this point in the next section. In contrast to the 2HDM, in the IDM and the HSM the parameters λ 2 and λ 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  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 → sγ -can also limit severely the allowed values of M H ± and tan β, and for this reason we choose to work in a 2HDM of type I (as mentioned already in Sec. II) where these constraints are weakest.
Indeed in type I, the lower limit on M H ± decreases from approximately 440 GeV for tan β = 1 to below 200 GeV for tan β = 1.5 [99].  [133]. Finally, for the HSM, we have considered a scenario with an additional 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].

C. 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λ BSM hhh can deviate from its SM predictionλ SM hhh , 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 massdegenerate additional scalars. Figure 9 shows the maximal possible size of the BSM deviation at two loops δR 2 in the plane of M Φ and tan β, considering now significantly larger ranges than in the previous section. Once these two parameters are fixed, the value ofM 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   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 λ hhh at future colliders these two-loop contributions could potentially be distinguishable experimentally in the future. For instance, supposing a deviation of 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 introductionwould 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 λ 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 MS calculations, we will encounter IR divergences in the limit m G,G ± → 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 λ 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 λ 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 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 λ 2 for the IDM and λ 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γγ 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.

VIII. 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 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 Z 2 symmetry). We have provided expressions for these corrections both in the 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 scaleM of the Z 2 symmetry in the 2HDM, ensuring the explicit decoupling of the BSM corrections, expressed in terms of OS-renormalised parameters, when taking the limitM → ∞ and using a relation of the form M 2 Φ =M 2 +λv 2 . We have also extended this result to the case of µ 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 (M orμ 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 O(∼ 20%) of the one-loop ones -noting the caveat for the IDM and the HSM that large values of the scalar quartic couplings λ 2 or λ 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 β and intermediate masses -i.e. M Φ between 500 and 900 GeV. Finally, we have carried out a preliminary estimate of the theoretical uncertainty associated with our two-loop computation ofλ hhh , taking this time the example of the IDM, and we found a conservative estimate of about 5% (of the total result). with µ denoting the regularisation scale. We also use the definition where Q is the renormalisation scale, defined as Q 2 = 4πe −γ E µ 2 . At one-loop order, we encounter the two integrals and from their finite (and -independent) part, we obtain the usual Passarino-Veltmann functions [109] as The following limits are also useful Because we neglect all external momenta, we only need one master integral at two loops, namely the sunrise integral (see e.g. [124]) . (A.9) Its finite, -independent, part is obtained as [111] I(x, y, z) ≡ lim where c xxx is a numerical constant defined as Li 2 being the dilogarithm function.

Appendix B: 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.

F F S for degenerate BSM scalar masses
The complete expressions for the derivatives (with the operators D 3 and D 4 ) of the 2HDM effective potential diagrams involving both the top quark and BSM scalars are: When working beyond O(M 4 t ), 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 (V.4) -it receives the following new one-loop terms

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.