\section{Introduction}
The discovery of a Higgs boson of mass 125 GeV at the CERN Large Hadron Collider~\cite{Chatrchyan:2012xdj,Aad:2012tfa} has completed the particle spectrum of the Standard Model (SM) of particle physics, and has established the Higgs sector as the origin of electroweak symmetry breaking (EWSB). Nevertheless, numerous open questions remain unsolved by the SM, and point to the need for Beyond-the-Standard-Model (BSM) physics. Among these issues, one can mention, on the theory side, for instance the hierarchy problem, the lack of a description of quantum gravity, or inflation. Additionally, a number of experimental results have contradicted prior assumptions -- upon which the SM was built -- such as for example (to cite only a few) the need to realise baryogenesis, the evidence for dark matter and dark energy, or the discovery that neutrinos have tiny masses. 

There is no doubt that the Higgs boson plays a central role in new physics, as well as in searches for it. Indeed, many (if not most) problems of the SM can be related to Higgs physics. In order to address these, BSM theories commonly feature a Higgs sector with extended particle content or modified dynamics -- as there is no compelling argument for it to be minimal. Therefore, shedding light on the structure of the Higgs sector is one of the most pressing tasks for particle physicists. Unfortunately, for the time being, all measurements of Higgs properties seem to be in excellent agreement with SM predictions~\cite{Aad:2019mbh,Sirunyan:2018koj}, and this forces us to investigate how new physics can appear to solve problems of the SM while simultaneously evading detection. A first, but unexciting, possibility is decoupling -- $i.e.$ that all new states are well beyond our experimental reach and their contributions to SM(-like) observables are vanishingly small due to the decoupling theorem~\cite{Appelquist:1974tg}. Another more interesting option is alignment without decoupling~\cite{Gunion:2002zf}, which occurs\footnote{Alignment is usually defined in models with several Higgs doublets, and is then the limit in which one of the CP-even mass eigenstates is colinear in field space with the total electroweak vacuum expectation value. } when the couplings of the 125-GeV Higgs boson to fermions and gauge bosons are SM-like at tree level, possibly because of some mechanism or symmetry of the Lagrangian -- see for instance Refs.~\cite{Dev:2014yca,Deshpande:1977rw,Silveira:1985rk,Barbieri:2006dq,Benakli:2018vqz,Benakli:2018vjk,Coyle:2019exn,Aiko:2020atr}. We will return to this latter case in the following.\footnote{Slightly non-alignment cases will be almost completely surveyed and, if nothing appears, will be excluded by the synergy of the HL-LHC and future precision measurements of the Higgs boson couplings such as at the ILC, FCCee, CEPC, CLIC, etc. \cite{Aiko:2020ksl}. }

Among the issues of the SM, one that has received considerable attention is the (gauge) hierarchy problem~\cite{Veltman:1976rt}, or in other words the need for a mechanism to stabilise the mass of the Higgs boson $M_h$ at the electroweak (EW) scale. Because the Higgs boson is a scalar, its mass is in the SM not protected by any symmetry, meaning that it typically receives huge -- quadratically divergent -- radiative corrections. In turn, this implies that a cancellation of tremendous -- and unrealistic -- accuracy would have to occur between tree-level and loop-level contributions to $M_h$. Solutions to the hierarchy problem usually involve some new symmetry or mechanism that protects the Higgs mass. For example, supersymmetry~\cite{Fayet:1977yc} introduces a new symmetry between bosons and fermions, so that scalar masses become protected by the chiral symmetry. Meanwhile, in composite Higgs models~\cite{Kaplan:1983fs,Kaplan:1983sm}, there are no fundamental scalars and the Higgs boson appears as a composite state of elementary fermions. But because new physics has so far eluded discovery, experimental constraints on most popular BSM scenarios have become more stringent (see for instance Refs.~\cite{ATLAS:2020uqv,ATLAS:2020ujn} and references therein) and the need to explore new avenues of research is increasingly evident. 

Another interesting possibility to address the hierarchy problem, which we will consider in this paper, comes from the concept of classical scale invariance (CSI), first proposed by Bardeen in Ref.~\cite{Bardeen:1995kv}. The idea of CSI is to forbid all mass terms in the Lagrangian, so that the theory is scale invariant at the tree ($i.e.$ classical) level. Of course, as this is only imposed at tree level, scale invariance is explicitly broken by radiative corrections, and EWSB can occur radiatively, via the Coleman-Weinberg mechanism~\cite{Coleman:1973jx} -- which was extended to multi-scalar cases by Gildener and Weinberg in Ref.~\cite{Gildener:1976ih}. The Higgs vacuum expectation value (VEV) $v$ and all particle masses (proportional to $v$) are then generated via dimensional transmutation. If the CSI model is taken to be UV complete, or in other words if scale invariance is assumed at the Planck scale, the absence of a tree-level Higgs mass parameter together with the generation of the Higgs mass entirely at loop level can ensure that there is no hierarchy problem -- see for instance Refs.~\cite{Bardeen:1995kv,Foot:2007as,Foot:2007iy,Meissner:2007xv,Iso:2009ss,Iso:2012jn}. On the other hand, models in which scale invariance is imposed at the EW scale can often not be extended up to the Planck scale, because of Landau poles appearing in the running of dimensionless parameters. However, even if a Landau pole appears at some intermediate scale $\Lambda$ and thus requires the existence of some further new physics beyond $\Lambda$, tree-level scalar masses are prohibited up to this scale thanks to the CSI. Therefore there is no hierarchy problem up to the scale $\Lambda$, and a problem only remains between $\Lambda$ and the Planck scale -- arguably a less severe issue than in the SM. The hierarchy problem is thus not entirely solved in such theories, although it can be somewhat alleviated. 

Nonetheless, motivated by a number of phenomenological considerations, we restrict our attention in this paper on the latter type of CSI models, defined with extended scalar sectors at the EW scale. A first interesting feature of such models is that, in the absence of any (BSM) mass term, new states cannot be decoupled, nor very heavy, and hence should be found in a foreseeable future -- if they exist. Additionally, in CSI models, the 125-GeV Higgs boson, which corresponds in the language of Ref.~\cite{Gildener:1976ih} to the ``scalon'' direction along which EWSB is radiatively realised, is automatically aligned at tree level -- as pointed out in Ref.~\cite{Lane:2019dbc} for a CSI variant of the Two-Higgs-Doublet Model (2HDM) -- thereby explaining current experimental results. 

A further discriminative characteristic of CSI models comes from the Higgs trilinear coupling $\lambda_{hhh}$. This coupling is of great importance as it controls the shape of the Higgs potential, and in turn, the nature and strength of the EW phase transition (EWPT).  In particular, one of the requirements to make the scenario of electroweak baryogenesis work~\cite{Sakharov:1967dj,Kuzmin:1985mm,Cohen:1993nk} is for the EWPT to be of strong first order. In turn, a large deviation of the Higgs trilinear coupling from its SM prediction is one of the main signs of a strong first-order EWPT -- for massive models, a deviation of at least 20\% was found to be necessary in Refs.~\cite{Grojean:2004xa,Kanemura:2004ch} -- together with the special spectrum of gravitational waves that would be produced by bubble collisions during the EWPT~\cite{Grojean:2006bp}.


Couplings of the Higgs boson -- and especially its trilinear coupling -- can potentially differ significantly from their SM values in BSM theories with extended scalar sectors, because of \textit{non-decoupling effects} from loop corrections involving the BSM states -- this was found at one loop, first in the 2HDM in Refs.~\cite{Kanemura:2002vm,Kanemura:2004mg}, and later in various BSM models with extra singlets~\cite{Kanemura:2015fra,Kanemura:2016lkz,He:2016sqr,Kanemura:2017wtm}, doublets~\cite{Kanemura:2015mxa,Arhrib:2015hoa,Kanemura:2016sos,Kanemura:2017wtm}, or triplets~\cite{Aoki:2012jj}. The limit when these BSM corrections are maximal is when the BSM scalars acquire large masses exclusively from the Higgs VEV (rather from some BSM mass parameter). Incidently, this is exactly what happens in CSI models, and as a matter of fact, the Higgs trilinear coupling has been found~\cite{Hill:2014mqa,Hashino:2015nxa,Agrawal:2019bpm} to be at one loop $(\lambda_{hhh}^\text{CSI})^{(1)}=5/3(\lambda_{hhh}^\text{SM})^{(0)}$ in \textit{any} CSI model (without mixing), where $(\lambda_{hhh}^\text{SM})^{(0)}$ is the tree-level SM prediction for $\lambda_{hhh}$ -- in other words $(\lambda_{hhh}^\text{CSI})^{(1)}$ deviates by $\sim 67\%$ from $(\lambda_{hhh}^\text{SM})^{(0)}$. The universality of this prediction is a particularly strong and unique property of CSI models, and follows (as we will see in section~\ref{SEC:GW}) from the especially simple form that the effective potential takes at one-loop order along the ``scalon'' direction. This value of $(\lambda_{hhh}^\text{CSI})^{(1)}$ also means that the EWPT is always of strong first order in CSI theories -- this, together with the synergy between the measurements of gravitational waves and of the Higgs trilinear coupling was investigated in Ref.~\cite{Hashino:2016rvx} (for the example of $O(N)$-symmetric CSI models).

One may then naturally wonder what is the current experimental status for the determination of $\lambda_{hhh}$ as well as future perspectives. It turns out that large deviations of $\lambda_{hhh}$ from its SM prediction are still allowed by experiments: indeed, at the present time the most stringent limits on the trilinear coupling come from the ATLAS collaboration using single-Higgs production searches with LHC Run 2 data and give $-3.2 < \lambda_{hhh}/\lambda_{hhh}^\text{SM} < 11.9$ at 95\% confidence level (CL)~\cite{ATL-PHYS-PUB-2019-009} (see also Refs.~\cite{Ferrari:2018akh,Aaboud:2018ftw} and~\cite{Sirunyan:2018iwt,Sirunyan:2018two} respectively for ATLAS and CMS results using double-Higgs production). However, these limits will be greatly improved at future colliders, thus making $\lambda_{hhh}$ an ideal target to search for BSM effects. We refer the interested reader to Ref.~\cite{deBlas:2019rxi} for a detailed discussion of the current prospects for the determination of $\lambda_{hhh}$ at upcoming experiments, and we only summarise here some of the main results. First of all, the high-luminosity upgrade of the LHC (HL-LHC) is expected to reach an accuracy of about $50-60\%$ with $3\text{ ab}^{-1}$ of data~\cite{Cepeda:2019klc} (see also Ref.~\cite{Chang:2019ncg}) while a high-energy extension (HE-LHC) might -- optimistically\footnote{Note that this result depends on the treatment of backgrounds -- see the discussions in Refs.~\cite{Homiller:2018dgu,Goncalves:2018yva} } -- be able to attain an $\mathcal{O}(15\%)$ precision~\cite{Goncalves:2018yva} at 68\% CL. Turning next to potential lepton colliders, the high-energy stages of the ILC can access the value of $\lambda_{hhh}$ respectively to 27\% accuracy at 500 GeV (with $4\text{ ab}^{-1}$ of data) and 10\% at 1 TeV ($8\text{ ab}^{-1}$ of data)~\cite{Fujii:2015jha} -- both values being once again at 68\% CL. Independently, CLIC could achieve a precision of approximately 10\% (at 68\% CL) by combining the data from all phases (up to 3 TeV)~\cite{Abramowicz:2016zbo,Charles:2018vfv,Roloff:2019crr}. Finally, a possible 100-TeV hadron collider would be able to determine $\lambda_{hhh}$ to a level of accuracy of 5-7\% at 68\% CL using $30\text{ ab}^{-1}$ of data~\cite{Goncalves:2018yva,Cepeda:2019klc,Chang:2018uwu}.


The universality of the value of $\lambda_{hhh}$ in CSI models at one loop also raises the question of the impact of two-loop corrections. Indeed, as the form of the effective potential is changed at two loops -- with new squared-logarithmic terms appearing, see $e.g.$ Refs.~\cite{vanderBij:1983bw,Ford:1992pn,Martin:2001vx} -- one may expect the effective Higgs trilinear coupling derived from it to be modified as well. Moreover, two-loop corrections to $\lambda_{hhh}$ have already been investigated in a number of extension of the SM with enlarged scalar sectors,\footnote{For completeness, we must also mention Refs.~\cite{Brucherseifer:2013qva,Muhlleitner:2015dua} in which the leading two-loop $\mathcal{O}(\alpha_s\alpha_t)$ corrections to $\lambda_{hhh}$ are computed, respectively for the Minimal-Supersymmetric-SM and Next-to-Minimal-Supersymmetric-SM, for the purpose of matching the precision at which the Higgs mass is calculated in these models. } namely in the Inert Doublet Model~\cite{Senaha:2018xek,Braathen:2019pxr,Braathen:2019zoh}, in the 2HDM~\cite{Braathen:2019pxr,Braathen:2019zoh}, and in a real-singlet extension of the SM~\cite{Braathen:2019zoh} (this last model actually corresponds to the $N$-scalar theory we will consider in this work for $N=1$). It is therefore also important to compute two-loop corrections to $\lambda_{hhh}$ in the CSI theories, to match the accuracy achieved in non-CSI scenarios and allow precise comparisons, as the computational method differs between the two types of models. 


In this paper, we have performed the first explicit\footnote{ A study including two-loop corrections to the effective potential, and its derivatives, via renormalisation-group improvement were performed in Ref.~\cite{Hill:2014mqa}. } calculation of dominant two-loop corrections to $\lambda_{hhh}$ in theories with CSI: we extend the known one-loop (effective-potential) calculation to two loops, so as to include the same level of corrections as in our previous works on non-CSI models~\cite{Braathen:2019pxr,Braathen:2019zoh}. We investigate in particular two types of BSM scenarios: first, $O(N)$-symmetric CSI models, and second, a CSI variant of the 2HDM~\cite{Lee:2012jn}.
Moreover, we include the renormalisation scheme conversion, from the {\ensuremath{\overline{\rm MS}}}\xspace scheme (in which the effective potential is obtained) to the OS scheme, in order to express our results in terms of physical quantities -- note that we will be referring to the {\ensuremath{\overline{\rm MS}}}\xspace and OS results respectively as $\lambda_{hhh}$ and $\hat\lambda_{hhh}$ to avoid confusions. We also take into account the requirement of generating the 125-GeV mass of the Higgs boson correctly, at the two-loop order. This yields a relation between the masses of all the states in the CSI models, which we use to constrain the BSM parameters (masses and couplings) appearing in our calculations. We also ensure that perturbative unitarity is maintained and that the EW vacuum remains the true vacuum of the effective potential for the parameter points we consider. As mentioned already, the Higgs trilinear coupling has a universal prediction at one loop in all CSI models -- $(\lambda_{hhh}^\text{CSI})^{(1)}=5/3(\lambda_{hhh}^\text{SM})^{(0)}$ -- and our main result in this paper is finding that the inclusion of two-loop corrections breaks this universality. Through our detailed numerical studies, we observe that the different theoretical constraints we impose on the considered scenarios limit the allowed range of values for $\lambda_{hhh}$ at two loops -- the strongest constraint coming from the requirement of correctly reproducing the 125-GeV Higgs mass. Nevertheless, the two-loop Higgs trilinear coupling retains dependence on the parameters of the BSM scalar sectors, and may provide an opportunity to distinguish different BSM scenarios with CSI.


This paper is organised as follows: in section~\ref{SEC:GW} we present general results, applicable to a wide range of CSI models, and we illustrate our calculational scheme with the example of the CSI Standard Model. Next, in sections~\ref{SEC:CSIO(N)} and~\ref{SEC:CSI2HDM}, we study the Higgs trilinear coupling and the relation between masses in scenarios of $N$-scalar models with $O(N)$ global symmetries and also in a CSI variant of the 2HDM. We discuss implications of our present in section~\ref{SEC:Discussion}, before concluding in section~\ref{SEC:Conclusion}. A number of appendices provide details about our computations and some results we deemed too long for the main text. Definitions of one- and two-loop functions are given in appendix~\ref{APP:loopfn}, complementary expressions for the $O(N)$-symmetric models are provided in appendix~\ref{APP:O(N)}, and results in the CSI-2HDM for general masses are presented in appendix~\ref{APP:CSI2HDM}.  Finally, we provide in appendix~\ref{APP:generic} generic results for the coefficients in the effective potential that can be applied for any CSI model without mixing. 

\section{Generalities about classically scale-invariant theories}
\label{SEC:GW}
We begin this paper by recalling known features of CSI models, before deriving some new results applicable for all such models. 


\subsection{Classically scale-invariant theories}
Classically scale-invariant models arise from the assumption that at some energy scale -- which we will in this paper consider to be at or near the EW scale -- the scalar potential exhibits scale invariance, by which we mean that all mass-dimensionful quantities vanish at the classical ($i.e.$ tree) level. Radiative corrections explicitly violate the CSI and allow the EW gauge symmetry to be broken dynamically, giving rise to a mass scale via dimensional transmutation. This mechanism to realise spontaneous EWSB radiatively was first studied for a model of a single scalar by S. Coleman and E. Weinberg~\cite{Coleman:1973jx}. However, it is known that it is impossible to explain the Higgs boson mass in this simplest single-scalar case: neglecting the top-quark effects, the original paper of Coleman and Weinberg~\cite{Coleman:1973jx} predicted a scalar of mass less than 10 GeV, and including the top quark renders the (one-loop) potential unstable. In order to reproduce the 125-GeV Higgs mass in models with CSI, it is therefore necessary to consider extended scalar sectors. The treatment of this general case of multi-scalar theories requires more care as one must also identify the direction in field space along which the EW gauge symmetry is broken. This is the purpose of the method devised by E. Gildener and S. Weinberg~\cite{Gildener:1976ih}, which we will review briefly now. 


\subsection{The Gildener-Weinberg method}
The tree-level scalar potential of a generic CSI theory with scalars $\{\varphi_i\}$ takes the form
\begin{align}
 V^{(0)}(\{\varphi_i\})=\Lambda_{ijkl}\varphi_i\varphi_j\varphi_k\varphi_l\,,
\end{align}
where the $\Lambda_{ijkl}$ are (dimensionless) quartic couplings -- tadpole, mass, and trilinear terms all have a non-zero mass dimension and must therefore vanish. Clearly, this potential does not admit any other minimum than the origin $\{\varphi_i\}=0$ and in order for the EW gauge symmetry to be broken spontaneously using the Coleman-Weinberg mechanism, it is necessary that the tree-level potential vanish along a ray in field space. This particular flat direction corresponds to the SM-like Higgs boson $\varphi_h$ (called \textit{scalon} in Ref.~\cite{Gildener:1976ih}) which carries the EW vacuum expectation value (VEV) $v$ and acquires a mass of 125 GeV at loop level. 

Along the flat direction, the scalar potential is generated at one-loop order, and it can be written as a function of the order parameter ($i.e.$ the Higgs field) $h\equiv\varphi_h-v$ using the supertrace formula~\cite{Jackiw:1974cv}
\begin{align}
\label{EQ:onelooppot_CW}
 V_{\text{eff}}(h)=\frac{1}{64\pi^2}\bigg\{&\tr\left[M_S^4(h)\left(\log\frac{M_S^2(h)}{Q^2}-\frac32\right)\right]-4\tr\left[M_f^4(h)\left(\log\frac{M_f^2(h)}{Q^2}-\frac32\right)\right]\nonumber\\
 &+3\tr\left[M_V^4(h)\left(\log\frac{M_V^2(h)}{Q^2}-\frac56\right)\right]\bigg\}\,,
\end{align}
where $M_{S,f,V}(h)$ denote respectively the scalar, fermion, and gauge-boson field-dependent mass matrices (at tree level), and $Q$ is the renormalisation scale. Note that, this expression is given in the {\ensuremath{\overline{\rm MS}}}\xspace renormalisation scheme. Additionally, compared to the paper of Gildener and Weinberg~\cite{Gildener:1976ih}, we have kept here all the terms of the supertrace formula, rather than just the logarithmic terms. As we will see in the next section, this has no impact on any physical quantity or discussion.

Due to the absence of any mass term in CSI models, the field-dependent masses of all particles in theory -- scalars, fermions, gauge bosons -- can be related to the field-independent masses as
\begin{equation}
\label{EQ:general_mass_form}
 m_X^2(h)=m_X^2\left(1+\frac{h}{v}\right)^2\,.
\end{equation} 
We can exploit this particularity of CSI models to write the one-loop effective potential under the simple form
\begin{equation}
\label{EQ:form_Veff_1l}
 V_{\text{eff}}(h)=\kappa\left[A^{(1)}\ (v+h)^4+B^{(1)}\ (v+h)^4\log\frac{(v+h)^2}{Q^2}\right]\,,
\end{equation}
where $\kappa$ is the loop factor defined in eq.~(\ref{EQ:loopfactor}) and 
\begin{align}
\label{EQ:def_AB}
 A^{(1)}&\equiv \frac{1}{4v^4}\left\{\tr\left[M_S^4\left(\log\frac{M_S^2}{v^2}-\frac32\right)\right]-4\tr\left[M_f^4\left(\log\frac{M_f^2}{v^2}-\frac32\right)\right]+3\tr\left[M_V^4\left(\log\frac{M_V^2}{v^2}-\frac56\right)\right]\right\}\,,\nonumber\\
 B^{(1)}&\equiv \frac{1}{4v^4}\left(\tr\left[M_S^4\right]-4\tr\left[M_f^4\right]+3\tr\left[M_V^4\right]\right)\,.
\end{align}
We will discuss in the following how two-loop corrections modify this form of $V_{\text{eff}}$, however, first we will review the one-loop result for the effective Higgs trilinear coupling that is computed with this potential.  


\subsection{The Higgs trilinear coupling at one-loop order}
We study in this paper corrections to an effective Higgs trilinear coupling, calculated as the third derivative of the effective potential with respect to the Higgs field, and evaluated at its minimum
\begin{equation}
 \lambda_{hhh}\equiv\frac{\partial^3V_{\text{eff}}}{\partial h^3}\bigg|_{h=0}\,.
\end{equation}
In such a computation based on the effective potential, lower derivatives of the potential also play an important role. First of all, from the first derivative of $V_{\text{eff}}$ we obtain the tadpole equation
\begin{align}
\label{EQ:1l_tadpole}
 \frac{\partialV_{\text{eff}}}{\partial h}\bigg|_{h=0}=2\kappa v^3\bigg(2A^{(1)}+B^{(1)}+2B^{(1)}\log\frac{v^2}{Q^2}\bigg)=0\,,
\end{align}
which allows us -- as $v\neq 0$ -- to eliminate\footnote{This confirms our earlier statement: a finite shift to the quantity $A^{(1)}$ plays no role physically, as $A^{(1)}$ is anyway eliminated using the tadpole equation. } the quantity $A^{(1)}$
\begin{align}
 A^{(1)}=-B^{(1)}\left(\frac12+\log\frac{v^2}{Q^2}\right)\,.
\end{align}
Next, the second derivative of $V_{\text{eff}}$ defines the effective-potential, or curvature, mass of the Higgs boson, which we denote $[M_h^2]_{V_{\text{eff}}}$. We have that
\begin{equation}
\label{EQ:1l_Higgscurvmass}
 [M_h^2]_{V_{\text{eff}}}\equiv\frac{\partial^2V_{\text{eff}}}{\partial h^2}\bigg|_{h=0}=2\kappa v^2\left(6A^{(1)}+7B^{(1)}+6B^{(1)}\log\frac{v^2}{Q^2}\right)=8\kappa v^2B^{(1)}\,.
\end{equation}
This equation has a very important consequence, unique to models with classical scale invariance: it relates on the left-hand side the Higgs curvature mass, close to its known pole mass of 125 GeV, to the tree-level masses of all particles in the theory through the quantity $B^{(1)}$ on the right-hand side. If the model we consider contains no extra fermions and gauge bosons but only additional scalars, we obtain~\cite{Gildener:1976ih,Lee:2012jn,Hashino:2015nxa,Lane:2018ycs} using eq.~(\ref{EQ:def_AB})
\begin{equation}
 4v^4B^{(1)}=8\pi^2v^2[M_h^2]_{V_{\text{eff}}}=\tr[M_S^4]-12m_t^4+6m_W^4+3m_Z^4\,,
\end{equation}
or equivalently
\begin{equation}
 \tr[M_S^4]=8\pi^2v^2[M_h^2]_{V_{\text{eff}}}+12m_t^4-6m_W^4-3m_Z^4\,,
\end{equation}
having taken into account the $N_c=3$ colour factor for the top quark and with $m_t$, $m_W$, and $m_Z$ the tree-level masses of the top quark and of the $W$ and $Z$ bosons respectively. Once one considers specific scenarios for the particle contents of BSM CSI theories, it becomes possible to obtain constraints on the masses of the new states -- as studied for example in Refs.~\cite{Hashino:2015nxa,Lane:2019dbc,Brooijmans:2020yij}. We will return to this type of relations among masses in the following, when we consider particular BSM models.  

Returning now to the one-loop potential ($c.f.$ eq.~(\ref{EQ:form_Veff_1l})) we find that the only two parameters on which it depends -- $A^{(1)}$ and $B^{(1)}$ -- are fixed by the tadpole equation and the knowledge of the Higgs mass. This implies that if we take further derivatives of the one-loop potential, the computed quantities are also fixed in terms of the Higgs curvature mass. Indeed, the effective Higgs trilinear coupling is at one-loop order
\begin{align}
 \lambda_{hhh}\equiv\frac{\partial^3V_{\text{eff}}}{\partial h^3}\bigg|_{h=0}=4 \kappa v \left(6 A^{(1)} + 13 B^{(1)} + 6 B^{(1)} \log\frac{v^2}{Q^2}\right)\,,
\end{align}
and using equations~(\ref{EQ:1l_tadpole}) and (\ref{EQ:1l_Higgscurvmass}), we can eliminate $A^{(1)}$ and $B^{(1)}$, and rewrite $\lambda_{hhh}$ in terms of the tree-level result for $\lambda_{hhh}$ in the Standard Model as~\cite{Hashino:2015nxa}
\begin{align}
\label{EQ:1l_res_lambdahhh}
 \lambda_{hhh}=40\kappa vB^{(1)}=\frac{5[M_h^2]_{V_{\text{eff}}}}{v}=\frac53(\lambda_{hhh}^{(0)})^\text{SM}\,.
\end{align} 
Similarly, we can compute a one-loop (effective) quartic coupling $\lambda_{hhhh}$
\begin{align}
\label{EQ:1l_res_lambdahhhh}
 \lambda_{hhhh}=\frac{\partial^4V_{\text{eff}}}{\partial h^4}\bigg|_{h=0}=&\ 4\kappa \left(6 A^{(1)} + 25 B^{(1)} + 6 B^{(1)} \log \frac{v^2}{Q^2}\right)=88\kappa B^{(1)}\nonumber\\
 =&\ \frac{11[M_h^2]_{V_{\text{eff}}}}{v^2}=\frac{11}{3}(\lambda_{hhhh}^{(0)})^\text{SM}\,.
\end{align}
Note that for these two last equations, we made use of the known tree-level results for the Standard Model
\begin{align}
 (\lambda_{hhh}^{(0)})^\text{SM}=\frac{3m_h^2}{v}\,,\qquad\text{and}\qquad  (\lambda_{hhhh}^{(0)})^\text{SM}=\frac{3m_h^2}{v^2}\,.
\end{align}
The predictions in equations~(\ref{EQ:1l_res_lambdahhh}) and (\ref{EQ:1l_res_lambdahhhh}) deviate significantly from the usual SM values -- for the case of $\lambda_{hhh}$ this effect is sufficiently large to be observed already at the HL-LHC. A further crucial aspect of these results is that, in CSI models, the Higgs trilinear and quartic coupling are, at one loop, \textit{completely independent} of the particle content of the theory. However, as we will see now, this universality is lost when two-loop corrections are included in the effective potential. 

\subsection{Two-loop corrections to the Higgs trilinear coupling}
\label{SEC:2llambda_hhh}
At two loops, corrections to the effective potential are obtained by computing one-particle irreducible vacuum bubble diagrams. In addition to logarithmic and non-logarithmic terms -- as are present at one loop -- such two-loop diagrams also contain squared logarithms, meaning that the effective potential takes the form 
\begin{align}
\label{EQ:form_Veff_1+2l}
 V_{\text{eff}}(h)=A\ (v+h)^4+B\ (v+h)^4\log\frac{(v+h)^2}{Q^2}+C\ (v+h)^4\log^2\frac{(v+h)^2}{Q^2}\,,
\end{align}
where the coefficient $C$ appears only at two-loop order. Expanding this expression in terms of one- and two-loop quantities, it becomes
\begin{align}
 V_{\text{eff}}(h)\equiv&\ \kappaV^{(1)}(h)+\kappa^2V^{(2)}(h)+\mathcal{O}(\kappa^3)\nonumber\\
         =&\ \kappa\bigg[A^{(1)}\ (v+h)^4+B^{(1)}\ (v+h)^4\log\frac{(v+h)^2}{Q^2}\bigg]\\
          &+ \kappa^2\bigg[A^{(2)}\ (v+h)^4+B^{(2)}\ (v+h)^4\log\frac{(v+h)^2}{Q^2}+C^{(2)}\ (v+h)^4\log^2\frac{(v+h)^2}{Q^2}\bigg]  +\mathcal{O}(\kappa^3)\nonumber\,.
\end{align}
While $A^{(1)}$ and $B^{(1)}$ are exactly the quantities given in eq.~(\ref{EQ:def_AB}), $A^{(2)}$, $B^{(2)}$, and $C^{(2)}$ must be derived from the two-loop diagrams contributing to the effective potential -- and will depend on the model one is considering.

Before investigating specific models, it is enlighting to repeat the procedure of the previous section, and derive general expressions for the derivatives of $V_{\text{eff}}$, using now the two-loop form in eq.~(\ref{EQ:form_Veff_1+2l}). First, the tadpole equation reads
\begin{equation}
 \frac{\partialV_{\text{eff}}}{\partial h}\bigg|_{h=0}=2 v^3 \bigg[2 A + B + 2 (B + C) \log\frac{v^2}{Q^2} + 2 C \log^2\frac{v^2}{Q^2}\bigg]=0\,,
\end{equation}
from which we can isolate $A$ as
\begin{align}
\label{EQ:A_2l}
 A=-\frac12B-(B+C)\log\frac{v^2}{Q^2}-C\log^2\frac{v^2}{Q^2}\,.
\end{align}
Next, the Higgs curvature mass is, after eliminating $A$
\begin{align}
\label{EQ:2l_Mh}
 [M_h^2]_{V_{\text{eff}}}=\frac{\partial^2V_{\text{eff}}}{\partial h^2}\bigg|_{h=0
 &=8 v^2 \bigg[B + C + 2 C \log\frac{v^2}{Q^2}\bigg]\,,
\end{align}
so that $B$ can be rewritten in terms of $[M_h^2]_{V_{\text{eff}}}$ and $C$
\begin{equation}
\label{EQ:B_2l}
 B=\frac{[M_h^2]_{V_{\text{eff}}}}{8v^2}-C\bigg(1+2\log\frac{v^2}{Q^2}\bigg)\,.
\end{equation}
Equation~(\ref{EQ:2l_Mh}) is also important as it constitutes the two-loop version of the relation between masses in CSI models -- which we first encountered at one loop in the previous section. Once we consider particular BSM models in sections~\ref{SEC:CSIO(N)} and~\ref{SEC:CSI2HDM}, we will see that this relation places stringent constraints on the allowed BSM mass and parameter ranges. 

Now, if we compute the effective Higgs trilinear coupling, we obtain at two loops
\begin{align}
\label{EQ:res_2l_lambdahhh}
 \lambda_{hhh}=\frac{\partial^3V_{\text{eff}}}{\partial h^3}\bigg|_{h=0

 &=\frac{5 [M_h^2]_{V_{\text{eff}}}}{v} + 32 C v\,,
\end{align}
once $A$ and $B$ are replaced using eqs.~(\ref{EQ:A_2l}) and~(\ref{EQ:B_2l}). 
We see here that because of the new term $C$, appearing from two loops in $V_{\text{eff}}$, non-universal, model-dependent corrections arise in the expression of $\lambda_{hhh}$. In practice we will be computing $C$ from the two-loop diagrams contributing to $V_{\text{eff}}$ and use this result to obtain $\lambda_{hhh}$ at two loops. However, we note that we can formally use the last equation to determine $C$ in terms of $[M_h^2]_{V_{\text{eff}}}$ and $\lambda_{hhh}$, as
\begin{align}
 C=\frac{1}{32v}\bigg(\lambda_{hhh}-\frac{5[M_h^2]_{V_{\text{eff}}}}{v}\bigg)\,,
\end{align}
We can finally use this expression of $C$ to compute $\lambda_{hhhh}$ at two loops in terms of $[M_h^2]_{V_{\text{eff}}}$ and $\lambda_{hhh}$ only
\begin{align}
 \lambda_{hhhh}=\frac{\partial^4V_{\text{eff}}}{\partial h^4}\bigg|_{h=0}
 = \frac{11 [M_h^2]_{V_{\text{eff}}}}{v^2}+192 C
 =\frac{6\lambda_{hhh}}{v}-\frac{19 [M_h^2]_{V_{\text{eff}}}}{v^2}\,.
\end{align}
This equation implies that once $\lambda_{hhh}$ is measured, the Higgs quartic coupling can be predicted in a universal way for all CSI theories, up to two loops. This being said, once three-loop corrections are included, we can expect this universality to be lost because of 
$\log^3$ terms appearing in $V_{\text{eff}}$ -- similarly to what happens with $\lambda_{hhh}$ once two-loop effects are taken into account. 


\subsection{Results for the classically scale-invariant Standard Model}
\label{SEC:CSISM}

\begin{figure}
 \centering
 \includegraphics[width=.8\textwidth]{./BKS_SM_2ldiags.pdf}
 \caption{Dominant two-loop diagrams contributing to the SM effective potential. }
 \label{FIG:SM2ldiags}
\end{figure}


We end this section by calculating the leading two-loop contributions to $V_{\text{eff}}$ and $\lambda_{hhh}$ for the classically scale-invariant version of the SM, which will also provide an occasion to present the renormalisation scheme conversion that we perform to obtain expressions in terms of physical quantities. We must emphasise that the CSI-SM is already known not to be a valid theory of Nature as it cannot reproduce the Higgs mass correctly~\cite{Coleman:1973jx}, however, its simplicity makes it the ideal setting to explain our calculational setup. 

The tree-level scalar potential of the CSI-SM does not contain a mass term, and reads simply
\begin{align}
 V^{(0)}=\lambda|\Phi|^4\,,
\end{align}
where the Higgs doublet $\Phi$ is defined as
\begin{equation}
\label{EQ:SMdoublet}
 \Phi=\begin{pmatrix}
  G^+\\ \frac{1}{\sqrt{2}}(v+h+iG)
 \end{pmatrix}\,.
\end{equation} 
$G$ and $G^\pm$ are respectively the neutral and charged would-be Goldstone bosons, while $h$ is the 125-GeV Higgs boson. 
At tree level, the field-dependent masses of the Higgs and Goldstone bosons read
\begin{equation}
 m_h^2(h)=3\lambda (v+h)^2\,,\qquad\text{and}\qquad m_G^2(h)=\lambda (v+h)^2\,,
\end{equation}
while the top-quark mass is (like for the usual SM)
\begin{align}
 m_t(h)=\frac{y_t}{\sqrt{2}}(v+h)\,.
\end{align}
The need for a flat direction in the tree-level potential imposes $\lambda=0$, which is also the result of solving the tree-level tadpole equation. 

The dominant contributions to the CSI-SM two-loop effective potential arise from the diagrams shown in figure~\ref{FIG:SM2ldiags} and involving the $SU(3)_c$ gauge coupling $g_3$ and the top Yukawa coupling $y_t$ -- the two EW gauge couplings $g_2$ and $g_Y$, and the light fermion Yukawa couplings are negligible in comparison. These contributions read (in the {\ensuremath{\overline{\rm MS}}}\xspace scheme)
\begin{align}
\label{EQ:CSISM_Vtwo}
 V^{(2)}(h)=&-4g_3^2m_t^2(h)\bigg(4J(m_t^2(h))-8m_t^2(h)-\frac{6J(m_t^2(h))^2}{m_t^2(h)}\bigg)\,,\nonumber\\
 &+3 y_t^2 \bigg[2 m_t^2(h) I (m_t^2(h), m_t^2(h), 0) + m_t^2(h) I (m_t^2(h), 0, 0) + J(m_t^2(h))^2\bigg]\,,
\end{align} 
where $J$ and $I$ are respectively one- and two-loop functions, for which definitions\footnote{Note also that the function we denote as $J$ is the same as the Passarino-Veltman function $A$. We choose the notation $J$ as $e.g.$ in Refs.~\cite{Ford:1992pn,Martin:2001vx} to avoid confusion with the quantity $A$ in the effective potential -- $c.f.$ eq.~(\ref{EQ:form_Veff_1+2l}). } and useful limits are given in appendix~\ref{APP:loopfn}.

In the previous section, we found that to obtain the two-loop corrections to $\lambda_{hhh}$, we need to compute $C$, $i.e.$ the coefficient of the $\log^2$ term in $V_{\text{eff}}$. For the CSI-SM, after expanding the expression in equation~(\ref{EQ:CSISM_Vtwo}), we find 
\begin{align}
 C^{(2)}=\frac{24g_3^2m_t^4}{v^4}-\frac{9m_t^6}{v^6}\,.
\end{align}
Using then equation~(\ref{EQ:res_2l_lambdahhh}), we obtain for the effective Higgs trilinear coupling in the CSI-SM, to leading two-loop order,
\begin{align}
\label{EQ:CSISM_lambdahhh_2l}
 \lambda_{hhh}=\frac{5[M_h^2]_{V_{\text{eff}}}}{v}+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2m_t^4}{v^3}-\frac{288m_t^6}{v^5}\bigg]\,.
\end{align}
At this point, we should also discuss the choice of renormalisation scheme for the parameters appearing in the results presented in this work. As the effective potential, from which we derive $\lambda_{hhh}$, is calculated in the {\ensuremath{\overline{\rm MS}}}\xspace scheme, the expression in equation~(\ref{EQ:CSISM_lambdahhh_2l}) should be understood as being written in terms of {\ensuremath{\overline{\rm MS}}}\xspace parameters. In the following, we will however choose to convert {\ensuremath{\overline{\rm MS}}}\xspace expressions into on-shell (OS) ones -- $i.e.$ written in terms of physical quantities and including also finite corrections from wave-function renormalisation (WFR). An OS version of the Higgs trilinear coupling can be obtained as
\begin{align}
 \hat\lambda_{hhh}=\left(\frac{Z_h^\text{OS}}{Z_h^{\ensuremath{\overline{\rm MS}}}\xspace}\right)^{\frac{3}{2}}\lambda_{hhh}=\left(1+\frac32\kappa\frac{d\Pi_{hh}^{(1)}}{dp^2}\bigg|_{p^2=0}+\mathcal{O}(\kappa^2)\right)\lambda_{hhh}\,,
\end{align}
where $Z_h^\text{OS}$ and $Z_h^{\ensuremath{\overline{\rm MS}}}\xspace$ are respectively the on-shell and {\ensuremath{\overline{\rm MS}}}\xspace scheme Higgs WFR constants, $\Pi^{(1)}_{hh}$ is the finite part of the one-loop Higgs boson self-energy, and with $\lambda_{hhh}$ expressed in terms of physical quantities. Note that it suffices to compute the WFR constants to one-loop order, because $\lambda_{hhh}$ is generated at one loop.

In the two-loop result for $\lambda_{hhh}$ in eq.~(\ref{EQ:CSISM_lambdahhh_2l}), $m_t$ and $g_3$ only enter at two loops, and therefore they require no scheme conversion here -- this would only produce corrections at three-loop order. However, it is necessary to include a one-loop scheme conversion for the Higgs {\ensuremath{\overline{\rm MS}}}\xspace VEV $v$ and curvature mass $[M_h^2]_{V_{\text{eff}}}$, both appearing from one loop. First, the {\ensuremath{\overline{\rm MS}}}\xspace VEV can be related to the Fermi constant at (leading) one loop as
\begin{align}
\label{EQ:VEVconv}
 v^2=\frac{1}{\sqrt{2}G_F}+\frac{3M_t^2}{16\pi^2}\big(2\llog M_t^2-1\big)\,.
\end{align}
For notational convenience, we also define an OS Higgs VEV as $v_\text{OS}^2\equiv(\sqrt{2} G_F)^{-1}$. Next, we recall that the Higgs pole ($i.e.$ physical) and curvature masses are related by the following formula
\begin{align}
 M_h^2=[M_h^2]_{V_{\text{eff}}}+\Pi_{hh}(p^2=M_h^2)-\Pi_{hh}(p^2=0)\,.
\end{align}
In the CSI-SM, the Higgs mass is generated at loop level, and thus this relation simplifies into
\begin{align}
 M_h^2=[M_h^2]_{V_{\text{eff}}}+\frac{M_h^2}{16\pi^2}\bigg[\frac{2M_t^2}{v_\text{OS}^2}\big(2+3\llog M_t^2\big)\bigg]\,.
\end{align}
Combining all these intermediate results together, we obtain finally
\begin{align}
 \hat\lambda_{hhh}=\frac{5M_h^2}{v_\text{OS}}+\frac{1}{16\pi^2}\frac{35}{2}\frac{M_h^2M_t^2}{v_\text{OS}^3}+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2M_t^4}{v_\text{OS}^3}-\frac{288M_t^6}{v_\text{OS}^5}\bigg]\,.
\end{align}
It is important to note that, because $M_h^2$ is generated at loop level, the second term, which is proportional to $\kappa M_h^2 M_t^2/v_\text{OS}^3$, is formally of \textit{two-loop} order -- only the first term in this expression is of one-loop order. 

Taking from the PDG~\cite{Zyla:2020zbs}, the following values for the physical inputs
\begin{align}
\label{EQ:SMinputs}
 M_h=125.10\text{ GeV}\,,\quad M_t=172.4\text{ GeV}\,, \quad G_F=1.1663787\cdot 10^{-5}\text{ GeV}^{-2}\,, \quad \alpha_s(M_Z)=0.1179\,,
\end{align}
we find numerically
\begin{align}
 \frac{5M_h^2}{v_\text{OS}}=&\ 317.8\text{ GeV}\nonumber\\
 \frac{5M_h^2}{v_\text{OS}}+\frac{1}{16\pi^2}\frac{35}{2}\frac{M_h^2M_t^2}{v_\text{OS}^3}+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2M_t^4}{v_\text{OS}^3}-\frac{288M_t^6}{v_\text{OS}^5}\bigg]=&\ 323.6\text{ GeV}\,,
\end{align} 
or in other words, the two-loop corrections to $\lambda_{hhh}$ yield a positive shift of the order of $2\%$.


\section{$N$-scalar models}
\label{SEC:CSIO(N)}
\subsection{Model definitions}
The first class of BSM CSI theories that we consider are extensions of the SM with $N$ additional real singlet scalar states. We will furthermore assume that these models are endowed with a global $O(N)$ symmetry, under which only the singlet scalars are charged while the SM particles do not carry any charge. Although this assumption greatly simplifies the expressions we obtain -- by reducing the number of different masses and coupling constants -- it is important to note that it does not significantly affect our discussion about corrections to the Higgs trilinear coupling. 
Such $O(N)$-symmetric CSI models have been discussed in many references -- see Refs.~\cite{Endo:2015ifa,Endo:2015nba,Endo:2016koi,Helmboldt:2016mpi,Fujitani:2017gma} -- from various points of view. In particular, the exact $O(N)$ symmetry provides a symmetry to stabilise dark matter candidates in these models. It was found in Ref.~\cite{Endo:2015nba} that there is an upper bound on the value of $N$ if one considers direct dark matter search results. However, in this paper, we choose to study $O(N)$ models exclusively from the point of view of Higgs physics and corrections to $\lambda_{hhh}$, and we do not consider potential dark matter aspects of these models. In other words, we consider here $O(N)$ as a phenomenological symmetry -- which greatly simplifies the study of Higgs physics in $N$-scalar theories -- rather than a symmetry to stabilise dark matter. 

We can write the scalar potential of $O(N)$-symmetric models with CSI as
\begin{align}
\label{EQ:O(N)_pot}
 V^{(0)}=&\ \lambda |\Phi|^4+\lambda_{\Phi S}\vec{S}^2|\Phi|^2+\frac{1}{4}\lambda_S(\vec{S}^2)^2\,,
\end{align}
where $\Phi$ is the SM-like Higgs doublet -- same as that shown in eq.~(\ref{EQ:SMdoublet}) -- and  $\vec{S}=(S_1\,, S_2\,, \cdots S_N)$ is a real scalar, singlet under the SM gauge groups and belonging to an $N$-dimensional representation of the global $O(N)$ symmetry group -- another way to see this is to think of $N$ real scalar singlets. Only the neutral component of the Higgs doublet acquires a non-zero VEV -- the singlets do not as this would break the global $O(N)$ symmetry. 

The scalar masses at tree-level are 
\begin{align}
 m_h^2(h)=3\lambda (v+h)^2\,,\qquad m_G^2(h)=\lambda (v+h)^2\,,\qquad m_{S_i}^2(h)=\lambda_{\Phi S} (v+h)^2\equiv m_S^2(h)\,.
\end{align} 
and the tree-level tadpole condition simply gives $\lambda=0$. Finally, it follows immediately from this equation that we can replace $\lambda_{\Phi S}$ by the tree-level BSM scalar mass $m_S$, as $\lambda_{\Phi S}=m_S^2/v^2$.

\subsection{Dominant two-loop corrections to the Higgs trilinear coupling}
\label{SEC:O(N)_lambdahhh}
\begin{figure}
 \centering
 \includegraphics[width=.6\textwidth]{./BKS_ON_2ldiags.pdf}
 \caption{Leading two-loop BSM diagrams contributing to the effective potential of the $O(N)$-symmetric CSI model -- in addition to the leading SM diagrams of figure~\ref{FIG:SM2ldiags}. }
 \label{FIG:O(N)2ldiags}
\end{figure}
We can obtain for the leading BSM contributions to $V^{(2)}$ by considering the new diagrams involving $\lambda_{\Phi S}$ and $\lambda_S$, shown in figure~\ref{FIG:O(N)2ldiags}. Their expressions are
\begin{align}
 V^{(2)}=&\ V^{(2)}_{hSS}+V^{(2)}_{SS}\,,\nonumber\\
 \text{with}\qquadV^{(2)}_{hSS}=&-N\lambda_{\Phi S}^2 (v+h)^2 I(m_S^2(h),m_S^2(h),0)\,,\nonumber\\
 V^{(2)}_{SS}=&\ \frac14N(N+2)\lambda_SJ(m_S^2(h))^2\,.
\end{align}
The squared-logarithm terms in these two-loop contributions are, in turn
\begin{align}
 C^{(2)}_{hSS}=&\ \frac{Nm_S^6}{v^6} \,,\nonumber\\
 C^{(2)}_{SS}=&\ \frac14N(N+2)\frac{\lambda_Sm_S^4}{v^4}\,.
\end{align}
Using once again equation~(\ref{EQ:res_2l_lambdahhh}), we can obtain the leading two-loop order {\ensuremath{\overline{\rm MS}}}\xspace result for $\lambda_{hhh}$
\begin{align}
 \lambda_{hhh}=\frac{5[M_h^2]_{V_{\text{eff}}}}{v}+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2m_t^4}{v^3}-\frac{288m_t^6}{v^5}+\frac{32Nm_S^6}{v^5}+8N(N+2)\frac{\lambda_Sm_S^4}{v^3}\bigg]\,.
\end{align}
Among the terms in the brackets, the first two (those involving $m_t$) are the SM-like contributions, while the last two are the new BSM effects. 
Note that, like in the SM case, $m_S$ and $\lambda_S$ only enter the expression of $\lambda_{hhh}$ at two loops, and hence we do not need to specify their renormalisation scheme. If we want to express the one-loop result in terms of $M_h^2$ and $v_\text{OS}$ we must however include scheme conversions for both $[M_h^2]_{V_{\text{eff}}}$ and $v$ and take into account WFR effects -- as discussed already in section~\ref{SEC:CSISM}. The only difference we must consider arises from the new momentum-dependent term in the Higgs self-energy in the $O(N)$-symmetric model, which reads
\begin{align}
 \Pi^{(1)}_{hh}(p^2)\supset-\frac{2Nm_S^4}{v^2}B_0(p^2,m_S^2,m_S^2)\,,
\end{align}
where the one-loop Passarino-Veltman function~\cite{Passarino:1978jh} $B_0$ is defined in appendix~\ref{APP:loopfn}. 
Adding this additional term, we find the OS expression for $\hat\lambda_{hhh}$ at leading two loops to be
\begin{align}
 \hat\lambda_{hhh}=&\ \frac{5M_h^2}{v_\text{OS}}+\frac{1}{16\pi^2}\frac{5M_h^2}{v_\text{OS}^3}\bigg[\frac{7}{2}M_t^2-\frac{1}{6}NM_S^2\bigg]\nonumber\\
 &+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2M_t^4}{v_\text{OS}^3}-\frac{288M_t^6}{v_\text{OS}^5}+\frac{32NM_S^6}{v_\text{OS}^5}+8N(N+2)\frac{\lambda_SM_S^4}{v_\text{OS}^3}\bigg]\,.
\end{align}
Once again, we emphasise that the terms in the first pair of brackets -- proportional to $\kappa M_h^2/v_\text{OS}^3$ -- are formally of two-loop order. 

In principle the above result for $\hat\lambda_{hhh}$ involves 3 BSM parameters, namely $M_S$, $\lambda_S$, and $N$. However, the classical scale invariance produces a further constraint on the BSM parameters, by relating the masses of the particles in the theory, as we will see in the following section.

\subsection{Mass relation}
\label{SEC:O(N)_massrel}
When investigating derivatives of the effective potential in general CSI models in section~\ref{SEC:2llambda_hhh}, we had found that the curvature mass of the Higgs boson is related to the $B$ and $C$ coefficients in $V_{\text{eff}}$. In terms of one- and two-loop quantities, this relation reads
\begin{align}
\label{EQ:2l_massrelation}
 8\pi^2[M_h^2]_{V_{\text{eff}}}=4v^2B^{(1)}+\frac{4v^2}{16\pi^2}\left[B^{(2)}+C^{(2)}\left(1+2\log\frac{v^2}{Q^2}\right)\right]\,.
\end{align}
As we have computed the effective potential in $O(N)$-symmetric CSI models, it is a simple matter to derive all the necessary ingredients to apply this relation to this class of models. In addition to $C^{(2)}$, which was given in the previous section, we also need the logarithmic terms at one and two loops -- $B^{(1)}$ and $B^{(2)}$. Including both SM-like and BSM terms, these are
\begin{align}
 B^{(1)}=&\ \frac{1}{4v^4}\left(Nm_S^4-12m_t^4+6m_W^4+3m_Z^4\right)\,,\nonumber\\
 B^{(2)}=&\ \frac{4g_3^2m_t^4}{v^4}\left[12\log\frac{m_t^2}{v^2}-16\right]+\frac{6m_t^6}{v^6}\left[-3\log\frac{m_t^2}{v^2}+8\right]\nonumber\\
 &+\frac{Nm_S^6}{v^6} \bigg[2\log \frac{m_S^2}{v^2}-4\bigg]+\frac12N(N+2)\frac{\lambda_Sm_S^4}{v^4}\bigg[\log \frac{m_S^2}{v^2}-1\bigg]\,.
\end{align}
Inserting these into equation~(\ref{EQ:2l_massrelation}), we obtain in terms of {\ensuremath{\overline{\rm MS}}}\xspace-renormalised parameters
\begin{align}
 8\pi^2v^2[M_h^2]_{V_{\text{eff}}}=&\ Nm_S^4-12m_t^4+6m_W^4+3m_Z^4\nonumber\\
 &+\frac{4m_t^4}{16\pi^2}\bigg[g_3^2\left(48\log\frac{m_t^2}{Q^2}-40\right)+\frac{m_t^2}{v^2}\left(-18\log\frac{m_t^2}{Q^2}+39\right)\bigg]\nonumber\\
 &+\frac{4m_S^4}{16\pi^2}\bigg[\frac{Nm_S^2}{v^2} \bigg(2\log \frac{m_S^2}{Q^2}-3\bigg)+\frac12N(N+2)\lambda_S\bigg(\log \frac{m_S^2}{Q^2}-\frac12\bigg)\bigg]\,.
\end{align}
The first line is simply the one-loop result, while the second and third lines are respectively the SM-like and BSM two-loop corrections to the relation between masses. 

In order to convert this result in terms of physical masses for the SM and BSM states, we require the Higgs and BSM scalar self-energies -- provided in appendix~\ref{APP:CSIO(N)} -- as well as equation~(\ref{EQ:VEVconv}) for the Higgs VEV. We obtain finally
\begin{align}
\label{EQ:O(N)_2l_massrelationOS}
 \frac{4\sqrt{2}\pi^2}{G_F}M_h^2=&\ NM_S^4-12M_t^4+6M_W^4+3M_Z^4\\
 &+M_h^2\bigg[\frac72 M_t^2-\frac16NM_S^2\bigg]+\frac{3M_t^4}{8\pi^2}\bigg[16g_3^2-\frac{6M_t^2}{v_\text{OS}^2}\bigg]+\frac{NM_S^4}{16\pi^2}\bigg[(N+2)\lambda_S+\frac{4M_S^2}{v_\text{OS}^2}\bigg]\,.\nonumber
\end{align}
If one considers this relation from the point of view of computing $M_S$ in terms of $\lambda_S$, $N$, and the known SM inputs, the two-loop corrections decrease the extracted value of $M_S$. However, in our numerical investigations, we will instead use this relation to compute $\lambda_S$ as a function of $M_S$ and $N$.   


\subsection{Numerical study}
At this point, we have derived all the analytical results necessary to study the possible magnitude of the leading two-loop corrections to the Higgs trilinear coupling in the $O(N)$-symmetric CSI model. Rather than examining the absolute value of $\hat\lambda_{hhh}$, we prefer to investigate by how much it deviates from the prediction in the (non-CSI) Standard Model, at the same leading two-loop order. Therefore, we present in the following results for the BSM deviation $\delta R$, defined as 
\begin{align}
\label{EQ:deltaR_def}
 \delta R&\equiv \frac{\hat\lambda_{hhh}^\text{BSM}-\hat\lambda_{hhh}^\text{SM}}{\hat\lambda_{hhh}^\text{SM}}\,,
\end{align}
where $\hat\lambda_{hhh}$ is the leading two-loop result for the Higgs trilinear coupling in the OS scheme and the SM -- its analytical expression can be found for instance in Ref.~\cite{Braathen:2019zoh} (see also Ref.~\cite{Senaha:2018xek,Braathen:2019pxr}). Taking the same input values as in eq.~(\ref{EQ:SMinputs}), we obtain numerically $\hat\lambda_{hhh}^\text{SM}\simeq 177.4\text{ GeV}$.

\begin{figure}
 \centering
 \includegraphics[width=.7\textwidth]{./BKS_ON_deltaRnonCSI_noconstraint.pdf}
 \caption{BSM deviation $\delta R$ -- defined in eq.~(\ref{EQ:deltaR_def}) -- of the Higgs trilinear coupling $\hat\lambda_{hhh}$ computed in the $O(N)$-symmetric CSI model with respect to its prediction in the usual ($i.e.$ non-CSI) SM. The blue (green) curves show the results at two loops for $N=1$ ($N=4$). Solid and dashed curves are plotted respectively for $\lambda_S=0$ and $\lambda_S=2$. The red dot-dashed line shows the one-loop result -- $i.e.$ comparing the one-loop CSI result of $5M_h^2/v_\text{OS}$ to the one-loop (non-CSI) SM result. The black dotted line corresponds to the comparison of the one-loop CSI result for $\lambda_{hhh}$ to the SM result at tree level -- $i.e.$ the 67\% deviation mentioned in the introduction. }
 \label{FIG:O(N)_deltaRnonCSI1}
\end{figure}

We present in figure~\ref{FIG:O(N)_deltaRnonCSI1} our results for $\delta R$ using only the expressions from section~\ref{SEC:O(N)_lambdahhh} -- in other words, we do not impose the relation between masses from eq.~(\ref{EQ:O(N)_2l_massrelationOS}) yet. The red dot-dashed line shows the one-loop result, comparing the one-loop CSI result $(\hat\lambda_{hhh})^{(1)}=5M_h^2/v_\text{OS}$ to the (leading) one-loop expression in the SM -- see for instance Ref.~\cite{Kanemura:2004mg}. We should point out here that this result differs from the $\sim 67\%$ mentioned in the introduction, and shown with the black dotted curve: we have $(\delta R)^{(1)}\simeq 82\%$. This is because the value of $\sim 67\%$ is obtained by comparing the one-loop $(\hat\lambda_{hhh}^\text{CSI})^{(1)}$ to the tree-level SM result for $(\hat\lambda_{hhh}^\text{SM})^{(0)}$, but once we include the one-loop contribution from the top quark in the SM expression we obtain a smaller value for $(\hat\lambda_{hhh}^\text{SM})^{(1)}<(\hat\lambda_{hhh}^\text{SM})^{(0)}$, and hence a larger deviation. We must emphasise that with our definition of $(\delta R)^{(1)}$ we are comparing \textit{one-loop effective-potential results}, which do not include any dependence on external momenta. It should however be noted that, as studied in Ref.~\cite{Kanemura:2004mg}, the momentum dependence in the SM result can be quite significant, and could modify significantly the value of the BSM deviation. In Ref.~\cite{Fujitani:2017gma}, momentum-dependent corrections were also included in the study of the Higgs trilinear coupling in $O(N)$-symmetric CSI models, and were found to be non-negligible.  
Nevertheless, we still choose to define $(\delta R)^{(1)}$ in terms of effective Higgs trilinear couplings, as this allows simpler comparisons with our new two-loop results, calculated with the effective potential. The result we find for $(\delta R)^{(1)}=(\hat\lambda_{hhh}^\text{CSI})^{(1)}/(\hat\lambda_{hhh}^\text{SM})^{(1)}-1$ remains totally independent of the particle content of the CSI model. 

Next, the blue and green curves in figure~\ref{FIG:O(N)_deltaRnonCSI1} are the values of $\delta R$ at two loops, for different choices of $N$ and $\lambda_S$: blue and green curves correspond respectively to $N=1$ and $N=4$; solid curves are made for $\lambda_S=0$ while the dashed ones are for $\lambda_S=2$. For the lowest values of $M_S$ -- $i.e.$ where the BSM effects are the smallest -- the two-loop deviation $(\delta R)^{(2)}$ is only slightly larger than at one loop, however, as one increases $M_S$, the BSM contributions grow rapidly. This growth is stronger for larger values of $N$ and $\lambda_S$, and, \textit{crucially}, the universality of the BSM shift found at one loop is entirely lost at two loops. Looking at concrete numerical values of $(\delta R)^{(2)}$, for example for $M_S=500\text{ GeV}$, we have on the one hand for $N=1$ and $\lambda_S=0$ a deviation of $\sim 94\%$, while on the other hand for $N=4$ and $\lambda_S=2$ we find $(\delta R)^{(2)}\sim 168\%$. Anticipating a little on the discussion in the next pages, we note that all parameter points in this figure~\ref{FIG:O(N)_deltaRnonCSI1} fulfill the criterion of tree-level perturbative unitarity~\cite{Lee:1977eg,Hashino:2016rvx}. 

\begin{figure}[h]
 \includegraphics[width=\textwidth]{./BKS_ON_compareCSInonCSI.pdf}
 \caption{Comparison of the BSM deviations $\delta R$ computed in $O(N)$-symmetric models, with (solid curves) and without (dashed curves) CSI, as a function of the pole mass $M_S$ of the BSM scalars. Results for $N=1$ and $N=4$ are shown respectively in the left- and right-hand plots. Red curves are one-loop values, blue and green curves are two-loop results for $\lambda_S=0$ (respectively for $N=1$ and $N=4$), and black and brown curves are two-loop results for $\lambda_S=2$ (respectively for $N=1$ and $N=4$).}
 \label{FIG:O(N)_compareCSInonCSI}
\end{figure}

Another interesting comparison we can perform is to investigate the difference in the BSM deviation in $\hat\lambda_{hhh}$ between CSI and non-CSI versions of $O(N)$-symmetric models. For this purpose, we derive in appendix~\ref{APP:nonCSIO(N)} the dominant two-loop corrections to $\hat\lambda_{hhh}$ in a non-CSI variant of the $O(N)$-symmetric model, and for generic values of $N$ -- this calculation follows the procedure developed in Refs.~\cite{Braathen:2019pxr,Braathen:2019zoh} (note also that the case of $N=1$ was already treated in Ref.~\cite{Braathen:2019zoh}, under the name of ``Higgs singlet model''). We present in figure~\ref{FIG:O(N)_compareCSInonCSI} a comparison of our results for the BSM deviation $\delta R$ computed with respect to the SM result in both CSI and non-CSI $O(N)$-symmetric models, for $N=1$ (left side) and for $N=4$ (right side). Note that, as there is no mass term in the CSI $O(N)$-symmetric model, we must set $\tilde\mu_S=0$ in the result for the non-CSI version of the model in order to be able to compare both variants consistently.  
The solid curves in fig.~\ref{FIG:O(N)_compareCSInonCSI} are obtained for the versions of the models with CSI, while the results for the non-CSI variants are given by the dashed curves. Additionally, in both plots, the red curves show one-loop values, while the blue and black (green and brown) curves are the two-loop results for $N=1$ ($N=4$) for $\lambda_S=0$ and $\lambda_S=2$. 

One clearly notices that CSI and non-CSI variants of this BSM model behave quite differently. First, for low BSM masses, the BSM deviation in the non-CSI case is minute, while in the CSI case the deviation is already large ($80\%$ or so) -- this is because of the difference between the normal SM and CSI SM. Then, when $M_S$ increases, the non-CSI corrections grow much faster than their CSI counterparts, as in the non-CSI models the one-loop corrections are proportional to $M_S^4$ -- $c.f.$ equation~(\ref{EQ:nonCSIO(N)_lambdahhh_1l}). 

\begin{figure}[h]
 \includegraphics[width=\textwidth]{./BKS_ON_unitarityandmassrelation.pdf}
 \caption{Regions of the $M_S$ and $\lambda_S$ parameter space of the $O(N)$-symmetric CSI  model allowed (light green) and excluded (light red) by the requirement of tree-level perturbative unitarity -- following Ref.~\cite{Hashino:2016rvx}. Additionally, the blue and red curves give the values of $\lambda_S$ computed at respectively one and two loops using equation~(\ref{EQ:O(N)_2l_massrelationOS}) as a function of $M_S$. The left pane is for $N=1$, while the right one is for $N=4$. }
 \label{FIG:O(N)_unitaritycontours}
\end{figure} 

Until this point, we have not taken into account the mass relation -- as derived in section~\ref{SEC:O(N)_massrel} -- nor the theoretical constraint from unitarity. For the latter, we choose to use as our criterion tree-level\footnote{We want to constrain the values of $\lambda_S$ and $M_S$, both of which appearing in $\hat\lambda_{hhh}$ only at two-loop order. Therefore, the difference from including higher-order (instead of tree-level) perturbative unitarity constraints should be formally a three-loop effect, and thus within the expected theoretical uncertainty of our result.} perturbative unitarity~\cite{Lee:1977eg}, and we can for the $O(N)$-symmetric model employ the results given in Ref.~\cite{Hashino:2016rvx}, in particular equation (22) therein. We show in figure~\ref{FIG:O(N)_unitaritycontours} the regions of parameter space --  in the $M_S$-$\lambda_S$ plane -- allowed (light green) and excluded (light red) under this criterion of tree-level perturbative unitarity, for both $N=1$ (left) and $N=4$ (right). Moreover, we superimpose on top of this plane the curves obtained when computing $\lambda_S$ as a function of $M_S$ (and $N$) using the mass relation we derived in eq.~(\ref{EQ:O(N)_2l_massrelationOS}) -- the blue line is for $N=1$ while the red one is for $N=4$. We recall also that negative values of $\lambda_S$ are forbidden, otherwise the potential would not be bounded from below. Although a wide range of BSM masses seems to be allowed by unitarity, once one takes into account the relation between masses in the CSI model -- or in other words once one imposes that $M_h=125\text{ GeV}$ -- only a narrow interval remains. Concretely, we find for $N=1$ that $509\text{ GeV}\lesssim M_S\lesssim 524\text{ GeV}$, and for $N=4$ that $365\text{ GeV}\lesssim M_S\lesssim 375\text{ GeV}$. In fact, at one-loop order -- see the first line of eq.~(\ref{EQ:O(N)_2l_massrelationOS}) --  only a single value of $M_S$ is allowed for a given $N$, however once two-loop contributions are included the presence of a new parameter $\lambda_S$ allows moderate variations of $M_S$. Conversely, as it only appears from two loops, $\lambda_S$ is not constrained very severely: indeed values as large as $\lambda_S=7.4$ (3.7) are possible for $N=1$ ($N=4$). 

Another important criterion that we use to verify the validity of the parameter points we consider is the true-vacuum condition~\cite{Politzer:1978ic} -- $i.e.$ the requirement that the EW vacuum is the true minimum of the (two-loop) potential and corresponds to a lower value of $V_{\text{eff}}$ than the origin of the potential (at $\varphi_h=v+h=0$). In CSI theories, as all states acquire their masses entirely from their coupling to the Higgs boson, all masses vanish at the origin $v+h=0$, and hence $V_{\text{eff}}(\varphi_h=0)=0$ and we only need to make sure that $V_{\text{eff}}(\varphi_h=v)$ is negative. The effective potential depends not only on the field $\varphi_h$, but also on the renormalisation scale $Q$ and therefore we examine the value of $V_{\text{eff}}(\varphi_h=v)$ for $Q=v$ and $Q=M_S$ -- the two choices that we believe to be the most natural. We find for $N=1$ (where $M_S\sim 515\text{ GeV}$) that the true-vacuum condition is fulfilled for $Q=M_S$ but not for $Q=v$, while for $N=4$ (where $M_S\sim 370\text{ GeV})$ the condition is verified for both values of $Q$. As it is unclear which value of $Q$ would be optimal (both values being around the EW scale anyway), the case of $N=1$ is not conclusively excluded by this criterion. 

\begin{figure}[h]
 \centering
 \includegraphics[width=\textwidth]{./BKS_ON_deltaRnonCSI.pdf}
 \caption{BSM deviation $\delta R$ -- defined in eq.~(\ref{EQ:deltaR_def}) -- of the Higgs trilinear coupling $\hat\lambda_{hhh}$ computed at two loops in the $O(N)$-symmetric CSI model with respect to its usual ($i.e.$ non-CSI) SM prediction. The blue curve (left side) is the two-loop result obtained for $N=1$, while the green curve (right side) is the two-loop result for $N=4$. Moreover, for this figure, $\lambda_S$ is calculated as a function of $M_S$ (for $N=1$ on the left side and $N=4$ on the right side) using equation~(\ref{EQ:O(N)_2l_massrelationOS}). The red dot-dashed line shows the one-loop result for $\delta R$, while the black dotted line is the comparison of the one-loop CSI result with the tree-level SM one. Finally, the light-red shaded regions are excluded by tree-level perturbative unitarity.}
 \label{FIG:O(N)_deltaRnonCSI_conslambdaS}
\end{figure}

Finally, we present in figure~\ref{FIG:O(N)_deltaRnonCSI_conslambdaS} the BSM deviations $\delta R$ we obtain once we take into account the mass relation and the constraint from perturbative unitarity. As for previous figures, the left and right panes are for $N=1$ and $N=4$ respectively, and the red dot-dashed lines show the one-loop result for $\delta R$ -- independent of $M_S$ and $N$. The blue and green curves give the two-loop results for $N=1$ and $N=4$. We emphasise that there is now only a single result for $\delta R$ at two loops because once $M_S$ and $N$ are fixed, the only remaining parameter -- $\lambda_S$ -- is fixed via equation~(\ref{EQ:O(N)_2l_massrelationOS}). The ranges of $M_S$ for these plots correspond to $\lambda_S$ varying between $0$ and $4\pi$, but the light-red shaded regions are excluded by the criterion of tree-level perturbative unitarity, as discussed above. 

Taking into consideration all these theoretical constraints, we find that the two-loop BSM deviation to $\hat\lambda_{hhh}$ must be within the intervals from $\sim 99\%$ to $\sim114\%$ and from $\sim 93\%$ to $\sim 109\%$ respectively for $N=1$ and $N=4$. This means that while the universality found at one loop for $(\delta R)^{(1)}$ is lost at two loops, only a limited range of values of $(\delta R)^{(2)}$ is possible. Moreover, our computation finds the Higgs trilinear coupling in CSI $O(N)$-symmetric model to be about twice as large as in the SM -- $(\delta R)^{(2)}\sim 100\%\pm 10\%$ -- and is larger than what is obtained at one loop, which makes the BSM deviation easier to access in experimental searches.

Before ending our discussion of the CSI $O(N)$-symmetric models, it is also important to comment on the experimental status of the parameter points we have considered. We have verified using \texttt{HiggsBounds}~\cite{Bechtle:2008jh,Bechtle:2011sb,Bechtle:2013wla,Bechtle:2015pma,Bechtle:2020pkv} that the points allowed by under the theoretical criteria discussed above are also allowed from the point of view of collider searches.  
The input values for \texttt{HiggsBounds} are produced by using the \textsc{Mathematica} package \texttt{SARAH}~\cite{Staub:2008uz,Staub:2009bi,Staub:2010jh,Staub:2012pb,Staub:2013tta} to create a \texttt{SPheno}-based~\cite{Porod:2003um,Porod:2011nf} spectrum generator for the CSI $O(N)$-symmetric models, with $N=1$ as well as $N=4$. 



\section{CSI Two-Higgs-Doublet Model}
\label{SEC:CSI2HDM}

\subsection{Model definition}
We now turn to the case of a CSI version of the Two-Higgs-Doublet Model, as devised in Ref.~\cite{Lee:2012jn} (see also Ref.~\cite{Funakubo:1993jg}). We will be assuming here that CP is conserved in the Higgs sector, and thus we can write the tree-level potential of the model in terms of two $SU(2)_L$ doublets $\Phi_1$ and $\Phi_2$ of hypercharge $1/2$ as
\begin{align}
 V^{(0)}=&\ \frac12\lambda_1|\Phi_1|^4+\frac12\lambda_2|\Phi_2|^4+\lambda_3|\Phi_1|^2|\Phi_2|^2+\lambda_4|\Phi_1^\dagger\Phi_2|^2+\frac12\lambda_5\left((\Phi_1^\dagger\Phi_2)^2+(\Phi_2^\dagger\Phi_1)^2\right)\,,
\end{align}
where all the quartic couplings $\lambda_i$ are real. 
To prevent the occurrence of tree-level flavour-changing neutral currents, we impose a discrete $\mathbb{Z}_2$ symmetry~\cite{Glashow:1976nt,Paschos:1976ay} under which the scalar doublets transform as  
\begin{equation}
 \Phi_1\xrightarrow{\mathbb{Z}_2} \Phi_1\,,\qquad\Phi_2\xrightarrow{\mathbb{Z}_2}-\Phi_2\,,
\end{equation}
while the fermions remain unaffected. 

We expand the two doublets as 
\begin{align}
 \Phi_i=\begin{pmatrix}
  w_i^+ \\ \frac{1}{\sqrt{2}}(v_i+h_i +iz_i)
 \end{pmatrix}\,,\qquad i=1,2\,,
\end{align}
where $v_i$ ($i=1,2$) are the vacuum expectation values of the neutral components of the doublets. These can be taken to be real -- owing to the assumption of CP conservation -- and furthermore verify the relation $v_1^2+v_2^2=v^2\simeq(246\text{ GeV})^2$. From the ratio of $v_1$ and $v_2$, we can also define the usual quantity $\tan\beta\equiv v_2/v_1$. 

Requiring both VEVs to be non-vanishing, we find for the tree-level tadpole equations
\begin{align}
 \lambda_1 v_1^2+\lambda_{345} v_2^2&=0\,,\nonumber\\
 \lambda_2 v_2^2+\lambda_{345} v_1^2&=0\,,
\end{align}
with the shorthand notation $\lambda_{345}\equiv\lambda_3+\lambda_4+\lambda_5$. 
These equations can be solved as
\begin{align}
 \lambda_1=-\lambda_{345}\tan^2\beta\,,\quad\text{and}\quad\lambda_2=-\lambda_{345}\cot^2\beta\,,
\end{align}
or equivalently as
\begin{align}
 \frac{\lambda_1}{\lambda_2}=\tan^4\beta\,,\quad\text{and}\quad \lambda_1\lambda_2=\lambda_{345}^2\,.
\end{align}

Next, mass eigenstates are obtained by rotating the gauge eigenstate fields with the angle $\beta$
\begin{align}
 \begin{pmatrix}
  w_1^\pm\\ w_2^\pm
 \end{pmatrix}=R(\beta)\begin{pmatrix}
  G^\pm\\H^\pm
 \end{pmatrix}\,,\quad
 \begin{pmatrix}
  z_1\\ z_2
 \end{pmatrix}=R(\beta)\begin{pmatrix}
  G\\A
 \end{pmatrix}\,,\quad
 \begin{pmatrix}
  h_1\\ h_2
 \end{pmatrix}=R(\beta)\begin{pmatrix}
  h\\H
 \end{pmatrix}\,,
\end{align} 
having defined
\begin{align}
 R(x)\equiv\begin{pmatrix}
  \cos x & -\sin x\\
  \sin x &  \cos x
 \end{pmatrix}\,.
\end{align}
Among these fields, $h,H$ are CP-even Higgs boson, $A$ is a CP-odd Higgs boson, and $H^\pm$ is the charged Higgs boson. As in the SM, $G$ and $G^\pm$ are respectively the neutral and charged would-be Goldstone bosons.

In this new basis, the tree-level scalar mass eigenvalues can be found to be 
\begin{align}
\label{EQ:CSI2HDM_treemasses}
 m_G^2&=m_{G^\pm}^2=m_h^2=0\,,\nonumber\\
 m_H^2=-\lambda_{345}v^2\,,\qquad m_A^2&=-\lambda_5v^2\,,\qquad m_{H^\pm}^2=-\frac12(\lambda_4+\lambda_5)v^2\,.
\end{align}
Interestingly, and to the difference of the usual (non-CSI) 2HDM, the CP-even mass matrix is already diagonal after the rotation of angle $\beta$. In other words, no additional rotation of angle $\alpha-\beta$ is necessary, and the CSI-2HDM is naturally \textit{aligned} (at tree level). Moreover, finding $m_h=0$ also tells us that $h$ is indeed the flat direction along which we want to work in what follows.

We also follow the common choice of trading the free Lagrangian quartic couplings $\lambda_3,\,\lambda_4,\,\lambda_5$ for the tree-level masses of the scalars. From eq.~(\ref{EQ:CSI2HDM_treemasses}), one has
\begin{align}
 \lambda_3=-\frac{m_H^2-2m_{H^\pm}^2}{v^2}\,,\qquad
 \lambda_4=\frac{m_A^2-2m_{H^\pm}^2}{v^2}\,,\qquad
 \lambda_5=-\frac{m_A^2}{v^2}\,.
\end{align}
We should note also that in the following, we neglect loop-induced deviations from the alignment found at tree level. Such effects have been studied for instance in Ref.~\cite{Lane:2019dbc} in the CSI 2HDM (see also Ref.~\cite{Braathen:2017izn} for the usual 2HDM) and were found to be relatively small for most of the parameter space, so that we consider them subleading before the two-loop corrections from the BSM scalar states. This enables us to write the effective potential under the form of equation~(\ref{EQ:form_Veff_1+2l}) and thus to follow the calculational procedure described in section~\ref{SEC:GW}. 

As in the usual 2HDM, the two scalar doublets in the CSI-2HDM can couple to quarks and leptons. Under the $\mathbb{Z}_2$ symmetry, four types of coupling assignments are possible, corresponding to the four types of 2HDMs~\cite{Barger:1989fj,Grossman:1994jb,Aoki:2009ha}. We note however that, for the major part of the discussion in this section, we do not need to specify a type as we only consider effects from the top quark -- only when verifying experimental constraints do we assume interactions of type I in order to have less severe limits from flavour physics than for instance for types II and Y~\cite{Misiak:2017bgg}.

\subsection{Dominant two-loop corrections to the Higgs trilinear coupling}
At one loop, the result for $\lambda_{hhh}$ is model-independent and was discussed in section~\ref{SEC:GW}. We can now consider the dominant two-loop contributions to the effective potential, and in turn to the Higgs trilinear coupling. We choose in this section to present results for the case of degenerate BSM scalar masses $M_H=M_A=M_{H^\pm}\equiv M_\Phi$, because the expressions for general values of the BSM scalar masses are quite long and cumbersome, and moreover because this will not affect the physics we discuss here. We provide the complete expressions for $V^{(2)}$, $\lambda_{hhh}$, and $\hat\lambda_{hhh}$ in appendix~\ref{APP:CSI2HDM_lambdahhh}. We should also note that as the BSM scalar masses will only appear in the Higgs trilinear coupling starting from the two-loop order, it does not make a difference whether we use pole or tree-level masses for the BSM scalars -- the difference in $\lambda_{hhh}$ from using one or the other is formally of three-loop order.

\begin{figure}[h]
 \centering
 \includegraphics[width=.8\textwidth]{./BKS_CSI2HDM_2ldiags.pdf}
 \caption{New dominant two-loop BSM diagrams contributing to the effective potential of the CSI-2HDM -- in addition to the leading SM diagrams of figure~\ref{FIG:SM2ldiags}. The diagrams in the second line all vanish in the limit of degenerate BSM scalar masses. }
 \label{FIG:CSI2HDM2ldiags}
\end{figure}

The BSM diagrams contributing to leading order to $V^{(2)}$ are shown in figure~\ref{FIG:CSI2HDM2ldiags}. In the limit of degenerate masses, their expressions read
\begin{align}
 V^{(2)}_{SSS}(h)=&-\frac{4m_\Phi^4(v+h)^2}{v^4}I(m_\Phi^2(h),m_\Phi^2(h),0)-\frac{6m_\Phi^4\cot^22\beta(v+h)^2}{v^4}I(m_\Phi^2(h),m_\Phi^2(h),m_\Phi^2(h))\,,\nonumber\\
 V^{(2)}_{SS}(h)=&\ \frac{12m_\Phi^2\cot^22\beta}{v^2}J(m_\Phi^2(h))^2\,,\nonumber\\
 V^{(2)}_{FFS}(h)=&\ \frac{6 m_t^2\cot^2\beta}{v^2}\Big[J(m_t^2(h))^2 - 3 J(m_\Phi^2(h)) J(m_t^2(h)) + (m_t^2(h) - m_\Phi^2(h)) I(0, m_\Phi^2(h),m_t^2(h))\nonumber\\
 &\hspace{2cm} + (2m_t^2(h) - m_\Phi^2(h)) I(m_\Phi^2(h), m_t^2(h), m_t^2(h)) \Big]\,.
\end{align}
From these equations we can extract the BSM logarithm-squared terms, and including also the SM-like contributions, we find in total for $C^{(2)}$
\begin{align}
\label{EQ:CSI2HDM_C2}
 C^{(2)}=&\ \frac{24g_3^2m_t^4}{v^4}-\frac{9m_t^6}{v^6}+\frac{m_\Phi^6}{v^6}(4+21\cot^22\beta)+ \frac{6 m_t^2\cot^2\beta}{v^6}\bigg(m_\Phi^4- 3m_\Phi^2m_t^2 - \frac32m_t^4 \bigg)\,.
\end{align}
Applying equation~(\ref{EQ:res_2l_lambdahhh}), we obtain the dominant two-loop corrections to $\lambda_{hhh}$, involving $g_3$, $m_t$, $m_\Phi$, and $\tan\beta$ as
\begin{align}
\label{EQ:CSI2HDM_2lres}
 \lambda_{hhh}=&\ \frac{5[M_h^2]_{V_{\text{eff}}}}{v}+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2m_t^4}{v^3}-\frac{288m_t^6}{v^5}+\frac{32m_\Phi^6}{v^5}(4+21\cot^22\beta)\nonumber\\
 &\hspace{3.5cm}+\frac{192m_t^2m_\Phi^4\cot^2\beta}{v^5}\left(1-\frac{3m_t^2}{m_\Phi^2}-\frac{3m_t^4}{2m_\Phi^4}\right)\bigg]\,.
\end{align}

Lastly, we can -- like in the previous sections -- express the one-loop result $5[M_h^2]_{V_{\text{eff}}}/v$ in terms of the Higgs pole mass and the physical Higgs VEV, and take into account finite WFR effects, to obtain an OS result. The only 2HDM-specific corrections come from the momentum-dependent BSM scalar contributions to the Higgs self-energy. These can be found to be
\begin{align}
 \Pi_{hh}^{(1)}(p^2)\supset-\sum_{\mathclap{\Phi=H,A,H^\pm}}\frac{2n_\Phi m_\Phi^4}{v^2}B_0(p^2,m_\Phi^2,m_\Phi^2)\,,
\end{align} 
where $n_{H,A}=1$ and $n_{H^\pm}=2$. From this, we have
\begin{align}
 \frac{d\Pi_{hh}^{(1)}}{dp^2}\bigg|_{p^2=0}\supset-\frac13\sum_{\mathclap{\Phi=H,A,H^\pm}}\frac{n_\Phi m_\Phi^2}{v^2}\,. 
\end{align}
Finally, we obtain in terms of OS-renormalised parameters
\begin{align}
\label{EQ:CSI2HDM_lambdahhh_OS}
 \hat\lambda_{hhh}=&\ \frac{5M_h^2}{v_\text{OS}}+\frac{1}{16\pi^2}\frac{5M_h^2}{v_\text{OS}}\bigg[\frac{7}{2}\frac{M_t^2}{v_\text{OS}^2}-\frac{2}{3}\frac{M_\Phi^2}{v_\text{OS}^2}\bigg]\nonumber\\
 &+\frac{1}{(16\pi^2)^2}\bigg[\frac{768g_3^2M_t^4}{v_\text{OS}^3}-\frac{288M_t^6}{v_\text{OS}^5}+\frac{32M_\Phi^6}{v_\text{OS}^5}(4+21\cot^22\beta)\nonumber\\
 &\hspace{2cm}+\frac{192M_t^2M_\Phi^4\cot^2\beta}{v_\text{OS}^5}\left(1-\frac{3M_t^2}{M_\Phi^2}-\frac{3M_t^4}{2M_\Phi^4}\right)\bigg]\,.
\end{align}



\subsection{Mass relation}

Next, we also consider the mass relation for the case of the CSI-2HDM, which will serve to relate the BSM scalar masses $M_H,\,M_A,\,M_{H^\pm}$ and $\tan\beta$. As in the previous section, we present here only expressions in the limit of degenerate scalar masses, and leave the complete results for appendix~\ref{APP:CSI2HDM_massrelation}. To derive corrections to $\lambda_{hhh}$, we have already used the expression of $C^{(2)}$ given in eq.~(\ref{EQ:CSI2HDM_C2}). We must now also derive results for $B$ at one- and two-loop orders. First, at one loop it is straightforward to calculate
\begin{align}
 B^{(1)}=&\ \frac{1}{4v^4}\big(4m_\Phi^4-12m_t^4+6m_W^4+3m_Z^4\big)\,.
\end{align}

At two loops, we have
\begin{align}
 B^{(2)}=&\ \frac{4g_3^2m_t^4}{v^4}\left[12\log\frac{m_t^2}{v^2}-16\right]+\frac{6m_t^6}{v^6}\left[-3\log\frac{m_t^2}{v^2}+8\right]+B^{(2)}_{SS}+B^{(2)}_{SSS}+B^{(2)}_{F\bar{F}S}\,,
\end{align}
where the first two terms are SM-like contributions and the latter three -- $B^{(2)}_{SS},\,B^{(2)}_{SSS},\,B^{(2)}_{F\bar{F}S}$ -- are the BSM ones. These can be found to be in the limit of degenerate BSM masses
\begin{align}
 B^{(2)}_{SS}=&\ \frac{24m_\Phi^6\cot^22\beta}{v^6}\left(\log\frac{m_\Phi^2}{v^2}-1\right)\,,\nonumber\\
 B^{(2)}_{SSS}=&\ \frac{2m_\Phi^6}{v^6}(4+9\cot^22\beta)\left(\log \frac{m_\Phi^2}{v^2}-2\right)\,,\nonumber\\
 B^{(2)}_{F\bar{F}S}=&\ \frac{6m_t^2\cot^2\beta}{v^6}\bigg\{ 2m_\Phi^4\left(\log\frac{m_\Phi^2}{v^2}-2\right)-6m_t^2m_\Phi^2\left(\log\frac{m_\Phi^2}{v^2}-1\right)-m_t^4\left(3\log\frac{m_t^2}{v^2}-8\right)\bigg\}\,.
\end{align}


Applying equation~(\ref{EQ:2l_massrelation}), the mass relation at two loops is, in terms of {\ensuremath{\overline{\rm MS}}}\xspace-renormalised quantities,
\begin{align}
\label{EQ:CSI2HDM_massrelation_msbar}
 8\pi^2v^2[M_h^2]_{V_{\text{eff}}}=&\ 4m_\Phi^4-12m_t^4+6m_W^4+3m_Z^4\nonumber\\
 &+\frac{m_t^4}{4\pi^2}\bigg[g_3^2\left(48\log\frac{m_t^2}{Q^2}-40\right)+\frac{m_t^2}{v^2}\left(-18\log\frac{m_t^2}{Q^2}+39\right)\bigg]\nonumber\\
 &+\frac{1}{4\pi^2}\Bigg\{\frac{24m_\Phi^6\cot^22\beta}{v^2}\left(\log\frac{m_\Phi^2}{Q^2}-\frac12\right)+\frac{2m_\Phi^6}{v^2}[4+9\cot^22\beta]\left(\log\frac{m_\Phi^2}{Q^2}-\frac32\right)\nonumber\\
 &\hspace{1.5cm}+\frac{6m_t^2\cot^2\beta}{v^2}\bigg[2m_\Phi^4\left(\log\frac{m_\Phi^2}{Q^2}-\frac32\right)-6m_t^2m_\Phi^2\left(2\log\frac{m_\Phi^2}{Q^2}-\frac12\right)\nonumber\\
 &\hspace{4cm}-m_t^4\left(3\log\frac{m_t^2}{Q^2}-\frac{13}{2}\right)\bigg]\Bigg\}\,.
\end{align} 
The first line in this equation is the one-loop result, while the additional lines are the two-loop corrections. Among these, the first (square) brackets are the SM-like contributions, and the remains terms (in the curly brackets) are the BSM terms, proper to the CSI-2HDM.  

To convert equation~(\ref{EQ:CSI2HDM_massrelation_msbar}) from the {\ensuremath{\overline{\rm MS}}}\xspace to the OS scheme, we require as previously the correction to the Higgs VEV -- $c.f$ equation~(\ref{EQ:VEVconv}) and the top-quark and Higgs-boson one-loop self-energies, as well as the self-energies for the BSM scalars $H,\,A,\,H^\pm$ (also at one loop). All necessary expressions are provided in appendix~\ref{APP:CSI2HDM_massrelation}. With these, we are able to obtain the following result, in terms only of physical parameters
\begin{align}
\label{EQ:CSI2HDM_2l_massrelationOS}
 \frac{4\sqrt{2}\pi^2}{G_F}M_h^2=&\ 4M_\Phi^4-12M_t^4+6M_W^4+3M_Z^4+M_h^2\left[\frac72M_t^2-\frac23M_\Phi^2\right]+\frac{3M_t^4}{8\pi^2}\left[16g_3^2-\frac{6M_t^2}{v_\text{OS}^2}\right]\nonumber\\
 &+\frac{M_\Phi^6}{4\pi^2v_\text{OS}^2}\bigg[4+3\bigg(7-2\sqrt{3}\pi\bigg)\cot^22\beta\bigg]+\frac{3M_t^2\cot^2\beta}{2\pi^2v_\text{OS}^2} \Bigg\{- M_\Phi^4 - 3 M_t^2 M_\Phi^2 + \frac{13}{2} M_t^4\nonumber\\
 &\hspace{.5cm}- 
 2 M_t^2 (M_\Phi^2 - 2 M_t^2) \mathfrak{Re}\bigg[f_B\left(\frac{M_\Phi^2}{2M_t^2} - \sqrt{\frac{M_\Phi^2}{4M_t^4}-1} \right) + f_B\left(\frac{M_\Phi^2}{2M_t^2} + \sqrt{\frac{M_\Phi^2}{4M_t^4}-1} \right)\bigg] \nonumber\\
 &\hspace{.5cm}- M_\Phi^2 (M_\Phi^2 - 2 M_t^2) \mathfrak{Re}\bigg[f_B\left(\frac12 - \sqrt{\frac14 - \frac{M_t^2}{M_\Phi^2}}\right) + f_B\left(\frac12 + \sqrt{\frac14 - \frac{M_t^2}{M_\Phi^2}}\right)\bigg]\Bigg\}\,,
\end{align}
where the function $f_B$ is defined in equation~(\ref{EQ:def_fB}), and $\mathfrak{Re}[x]$ denotes the real part of a given quantity $x$.

\subsection{Numerical study}

\begin{figure}[h]
 \centering
 \includegraphics[width=.8\textwidth]{./BKS_CSI2HDM_deltaRnonCSI.pdf}
 \caption{BSM deviation $\delta R$ -- defined in eq.~(\ref{EQ:deltaR_def}) -- of the Higgs trilinear coupling $\hat\lambda_{hhh}$ computed at one loop (dot-dashed line) and at two loops (solid and dashed curves) in the CSI-2HDM with respect to its SM prediction as a function of the degenerate mass of the BSM scalars $M_\Phi=M_H=M_A=M_{H^\pm}$. The different grey, teal, blue, and green two-loop curves correspond respectively to $\tan\beta$ equal to 1, 1.1, 1.2, and 1.4. As in the previous section, the black dotted line is obtained by comparing the one-loop CSI result for $\lambda_{hhh}$ to the SM tree-level result.}
 \label{FIG:CSI2HDM_deltaRnonCSI1}
\end{figure}

We have now derived all the results we require to study numerically the leading two-loop corrections to the Higgs trilinear coupling in the CSI-2HDM. We follow here the same outline as we did for the $O(N)$-symmetric model, and we provide results for the same type of BSM deviation $\delta R$ -- as defined in eq.~(\ref{EQ:deltaR_def}) -- comparing the prediction for the Higgs trilinear coupling in the CSI-2HDM to that in the usual (non-CSI) SM. Moreover, we will consider the three BSM scalars $H,\,A,\,H^\pm$ to be degenerate in mass -- this ensures that the $\rho$ parameter remains close to one, and furthermore this conveniently reduces the number of BSM mass scales to one. 

We present first in figure~\ref{FIG:CSI2HDM_deltaRnonCSI1} the deviation $\delta R$ at one loop (red line) and two loops (grey, teal, blue, and green curves) for different values of $\tan\beta$. We employ here the expression in equation~(\ref{EQ:CSI2HDM_lambdahhh_OS}), without for now taking into account the relation between masses. As we had seen for the $N$-scalar model, while at one loop the BSM deviation in $\hat\lambda_{hhh}$ has a universal value of $\delta R\sim 82\%$, the two-loop corrections introduce new dependences on both $M_\Phi$ and $\tan\beta$, thereby spoiling the universality of $\delta R$. We consider in this plot different values of $\tan\beta$, ranging from $\tan\beta=1$ to $\tan\beta=1.4$, which we find to be the largest value of $\tan\beta$ for which tree-level perturbative unitarity~\cite{Lee:1977eg,Kanemura:1993hm} is maintained all the way to $M_\Phi=500\text{ GeV}$ -- see figure~\ref{FIG:CSI2HDM_unitaritycontours} and the discussion in the following. Considering for instance $M_\Phi=500\text{ GeV}$, we find $(\delta R)^{(2)}=134\%$ for $\tan\beta=1$, and $(\delta R)^{(2)}=162\%$ for $\tan\beta=1.4$, in stark contrast to the one-loop result. 

At this point, we should comment briefly on our interpretation of the quantity $\tan\beta$. While $\tan\beta$ is originally defined as the ratio of the VEVs of the two scalar doublets in the model, it appears in two-loop corrections to the Higgs trilinear coupling via the dependence of four-point scalar interactions at tree level and of the coupling between BSM scalars and the top quark. It should be emphasised however that because $\tan\beta$ appears only from the two-loop order, any scheme difference in how its value could be extracted from some experimental data (in the event a 2HDM Higgs sector is discovered) only has an impact from three-loop order and beyond -- $i.e.$ below the level of accuracy to which we are working. 


\begin{figure}[h]
 \includegraphics[width=.8\textwidth]{./BKS_2HDM_compareCSInonCSI.pdf}
 \caption{Comparison of the BSM deviations $\delta R$ computed in versions of the 2HDM, with (solid curves) and without (dashed curves) CSI, as a function of the degenerate pole mass $M_\Phi$ of the BSM scalars. Red curves present one-loop values, while grey and green curves are two-loop results for $\tan\beta=1$ and $\tan\beta=1.4$ respectively.}
 \label{FIG:2HDM_compareCSInonCSI}
\end{figure}

Next, it is also interesting to compare predictions for the Higgs trilinear coupling in variants of the 2HDM with or without CSI. For this purpose, we present in figure~\ref{FIG:2HDM_compareCSInonCSI} the BSM deviations in $\hat\lambda_{hhh}$ computed with respect to the usual-SM prediction in both variants of the 2HDM. For the usual 2HDM, we use the analytic expression for $\hat\lambda_{hhh}$ derived in Ref.~\cite{Braathen:2019pxr,Braathen:2019zoh} -- see specifically equation~(5.15) in Ref.~\cite{Braathen:2019zoh} -- in which we set the mass term $\tilde M$ (representing the scale of the soft breaking of the $\mathbb{Z}_2$ symmetry acting on the Higgs doublets) to zero to allow the comparison with the CSI scenario. 

We can observe that the Higgs trilinear coupling behaves very differently in CSI and non-CSI versions of the 2HDM. The main discrepancy occurs at one-loop level: indeed, while in the CSI-2HDM $(\delta R)^{(1)}$ is constant, in the usual (non-CSI) 2HDM, one-loop corrections involving the BSM scalars are minute in the low-mass range, but grow rapidly (as $M_\Phi^4$) and reach $\simeq 240\%$ for $M_\Phi=500\text{ GeV}$. At two loops, the radiative corrections to $\hat\lambda_{hhh}$ behave quite similarly in the CSI and non-CSI scenarios -- the main difference arising from effects due to the scheme conversion of the one-loop BSM corrections in the non-CSI 2HDM case ($c.f.$ Ref.~\cite{Braathen:2019zoh}).  

\begin{figure}[h]
 \includegraphics[width=\textwidth]{./BKS_CSI2HDM_unitarityandmassrelation.pdf}\
 \caption{Regions of the $\tan\beta$ and $M_\Phi$ parameter space of the CSI-2HDM allowed (light green) and excluded (light red) by the requirement of tree-level perturbative unitarity~\cite{Kanemura:1993hm}. Additionally, the red and blue curves give the values of $M_\Phi$ computed at respectively one and two loops using equation~(\ref{EQ:CSI2HDM_2l_massrelationOS}) as a function of $\tan\beta$. The right-hand side plot is an enlargement of the left-hand one. }
 \label{FIG:CSI2HDM_unitaritycontours}
\end{figure} 

At this point, we should discuss the theoretical constraints on the BSM parameters -- $M_\Phi$ and $\tan\beta$ -- coming from unitarity and from the mass relation in eq.~(\ref{EQ:CSI2HDM_2l_massrelationOS}). For the former, we choose to take tree-level perturbative unitarity~\cite{Lee:1977eg} as our criterion, and we employ results from Ref.~\cite{Kanemura:1993hm} (see also Ref.~\cite{Akeroyd:2000wc}). Figure~\ref{FIG:CSI2HDM_unitaritycontours} shows the allowed (light green) and excluded (light red) regions of the CSI-2HDM parameter space, in the $M_\Phi-\tan\beta$ plane. Additionally, the values of $M_\Phi$ that are extracted from the equation~(\ref{EQ:CSI2HDM_2l_massrelationOS}) are given by the dashed red (at one loop) and solid blue (at two loops) curves. On the one hand, at one loop $\tan\beta$ does not appear in the mass relation and hence at a constant value is found for $M_\Phi$. On the other hand, at two loops, the result for $M_\Phi$ obtained from eq.~(\ref{EQ:CSI2HDM_2l_massrelationOS}) becomes a function of $\tan\beta$. Along the two-loop (blue) curve, we find that the maximal possible value of $\tan\beta$ is 1.855, which corresponds to a maximal value of the degenerate mass of the BSM scalars $M_\Phi\leq 382.2 \text{ GeV}$. As for the $O(N)$-symmetric models, we have also verified that the true-vacuum condition is respected in the CSI-2HDM: we have found that the value of the two-loop potential is, at the EW minimum, indeed lower than the value of the potential at the origin, and this is true both for $Q=v$ and $Q=M_\Phi\sim 378\text{ GeV}$. 

\begin{figure}[h]
 \centering
 \includegraphics[width=.8\textwidth]{./BKS_CSI2HDM_deltaRnonCSI_tanbcons.pdf}
 \caption{BSM deviation $\delta R$ -- defined in eq.~(\ref{EQ:deltaR_def}) -- of the Higgs trilinear coupling $\hat\lambda_{hhh}$ computed at two loops in the CSI-2HDM with respect to its SM prediction as a function of the degenerate mass of the BSM scalars $M_\Phi=M_H=M_A=M_{H^\pm}$. The value of $\tan\beta$ in the point in this figure is computed as a function of $M_\Phi$ using equation~(\ref{EQ:CSI2HDM_2l_massrelationOS}). Finally, the maximal value of $M_\Phi$ (equivalently the maximal value of $\tan\beta$) in this plot is constrained by the requirement of tree-level perturbative unitarity -- $c.f.$ figure~\ref{FIG:CSI2HDM_unitaritycontours}. The red dot-dashed and black dotted lines show the comparison of the one-loop CSI value of $\lambda_{hhh}$ with respectively the one-loop effective-potential and tree-level results in the SM.} 
 \label{FIG:CSI2HDM_deltaRnonCSIconstanbeta}
\end{figure}

Having considered the constraints on $M_\Phi$ and $\tan\beta$, we can finally investigate the theoretically-allowed values of the Higgs trilinear coupling. Figure~\ref{FIG:CSI2HDM_deltaRnonCSIconstanbeta} therefore presents our results for the deviation $\delta R$ at one loop (red dot-dashed line) and at two loops (blue solid curve), once $\tan\beta$ is extracted as a function of $M_\Phi$ from the mass relation in eq.~(\ref{EQ:CSI2HDM_2l_massrelationOS}). As can be seen in figure~\ref{FIG:CSI2HDM_unitaritycontours}, the possible values of $M_\Phi$ are all within the narrow interval from $373.9\text{ GeV}$ to $382.1\text{ GeV}$. We find that $\delta R$ varies between $89.9\%$ at the lowest, and $112.5\%$ at most. In other words, as for the $O(N)$-symmetric model earlier, we obtain a larger deviation of $\hat\lambda_{hhh}$ from its deviation at two loops, and the total BSM deviation is (approximately) in the range $100\pm10\%$. 

We have verified that the corresponding parameter points are not excluded by experimental searches for BSM scalars. Assuming Yukawa interactions of type I (to benefit from less severe constraints from flavour physics than for types II or Y, see $e.g.$ Ref.~\cite{Misiak:2017bgg}), we create a \texttt{SPheno}-based spectrum generator for the CSI-2HDM using \texttt{SARAH}, which we in turn used to provide inputs for \texttt{HiggsBounds} to verify the experimental status of the parameter choices we took. We ensured that the Lagrangian inputs provided to \texttt{SPheno} yielded the corrected spectrum of scalar masses -- a 125-GeV SM-like Higgs boson and BSM scalars at 374 GeV, with an allowed error of 1 GeV -- using the \textsc{Mathematica} package \texttt{SSP}~\cite{Staub:2011dp}. 

\section{Discussion}
\label{SEC:Discussion}
The computations presented in this paper allow studying the Higgs trilinear coupling in models with classical scale invariance to the same level of accuracy as for non-CSI extensions of the SM -- as done for instance in Refs.~\cite{Braathen:2019pxr,Braathen:2019zoh,Senaha:2018xek} (and appendix~\ref{APP:nonCSIO(N)} of this work). Classical scale invariance is an attractive concept for model building, and may relate to the solution to the hierarchy problem. It is therefore of paramount importance to find new ways to investigate CSI models, in complement of collider and dark matter searches -- see $e.g.$ Refs~\cite{Endo:2015ifa,Endo:2015nba,Endo:2016koi,Helmboldt:2016mpi,Fujitani:2017gma,Lane:2019dbc,Brooijmans:2020yij}.  

Due to the simple form of the potential at one loop, containing only two free quantities -- recall equation~(\ref{EQ:form_Veff_1l}) -- $\hat\lambda_{hhh}$ is found to be predicted universally at the one-loop level, in \textit{all} CSI models. This universality of $(\hat\lambda_{hhh})^{(1)}$ is lifted by the inclusion of two-loop effects, which thus allows distinguishing different scenarios of CSI models. We found that the new corrections at two loops are quite significant, leading to a further 5-30\% deviation of $\hat\lambda_{hhh}$ from the (non-CSI) SM prediction, compared to the one-loop CSI result. We should emphasise that this is by no means a problem from the point of view of the validity of the perturbative loop expansion: indeed the two-loop corrections to $\hat\lambda_{hhh}$ involve new parameters that are not present at one loop -- namely $M_S,\,\lambda_S,\,N$ for the $O(N)$-symmetric models and $M_H,\,M_A,\,M_{H^\pm},\,\tan\beta$ in the CSI-2HDM. Furthermore, as there are numerous arguments in favour of extended scalar sectors, it is interesting to note that the two-loop contributions to the Higgs trilinear coupling receives from scalars are always positive, which makes it deviate further from the SM value so that it is more easily accessible in collider searches than what is expected from the one-loop result. 

We have also paid particular attention to the relation among masses arising in theories with CSI, which is often over-looked in the literature. More than being simply bounded from above, the masses of BSM states in CSI models are actually strongly constrained by the requirement of reproducing the correct 125-GeV mass of the Higgs boson -- which is entirely generated at loop level. In the scenarios considered in the main text of this paper, we imposed additional constraints -- a global $O(N)$ symmetry for the $N$-scalar model, and mass degeneracy of the BSM scalars in the CSI-2HDM -- which resulted in severe bounds on the allowed BSM masses.\footnote{Admittedly, if these restrictions were to be loosened (as $e.g.$ in appendix~\ref{APP:CSI2HDM}), there would be more freedom to vary the BSM scalar masses. Nevertheless, the masses would remain related so as to ensure the correctness of the Higgs mass. } In turn, this allows us to find predictions for $\hat\lambda_{hhh}$ within a narrow range: in both types of BSM models, we obtain $90\%\lesssim (\delta R)^{(2)}\lesssim 115\%$. Looking back at figures~\ref{FIG:O(N)_compareCSInonCSI} and~\ref{FIG:2HDM_compareCSInonCSI}, it is worth noting that the values obtained for $(\delta R)^{(2)}$ for the allowed BSM mass ranges are very close to those found for the same masses in the non-CSI counterparts of the scenarios. Taking the example of the 2HDM, while in the CSI version of the model we have $90\%\lesssim (\delta R)^{(2)}\lesssim 110\%$, in the non-CSI variant and taking the non-decoupling limit ($i.e.$ $\tilde M=0$) the value of $(\delta R)^{(2)}$ computed at $M_\Phi\simeq 374\text{ GeV}$ varies from $\sim 90\%$ (for $\tan\beta=1$) to $105\%$ (for $\tan\beta=1.9$, close to the upper limit from perturbative unitarity). This implies that the magnitude of a potential BSM deviation in the Higgs trilinear coupling may not provide a clear indication of whether the associated new physics exhibits classical scale invariance -- instead, a better criterion might come from whether masses and couplings fulfill the mass relations imposed by CSI.
 
Rather than considering only a single type of observable, a powerful strategy to investigate CSI theories, and to discriminate theories with or without CSI, would be to study the correlation between several observables. In particular, synergies between measurements of the Higgs trilinear coupling and direct searches of new particles at colliders may prove crucial to ascertain whether Nature exhibits classical scale invariance or not. Indeed, we should recall that due to the relation among masses in CSI models, the new BSM states cannot be arbitrarily heavy (even in non-degenerate scenarios) and will be within reach of experiments in a foreseeable future (see for instance the discussions for the CSI-2HDM in Refs.~\cite{Lane:2019dbc,Brooijmans:2020yij}). In the eventuality that a new scalar state is found in direct searches, determining the value of the Higgs trilinear coupling and verifying whether the relation among masses is fulfilled or not can potentially provide clear evidence of the CSI nature of the underlying theory. Furthermore, if one considers a model with charged scalars ($e.g.$ $H^\pm$ in the CSI-2HDM), sizeable one-loop corrections to the Higgs-to-two-photon ($h\gamma\gamma$) coupling can also be expected, as studied in Ref.~\cite{Hashino:2015nxa}. Such effects could also help confirm the existence (or absence) of CSI in BSM models. In this context, depending on the achieved accuracy of the determination of the $h\gamma\gamma$ coupling at future colliders -- possibly down to percent level, see $e.g.$ Refs.~\cite{Fujii:2015jha,Fujii:2017vwa} for prospects at the ILC -- it may become important (if not unavoidable) to compute theory predictions for this coupling at two-loop level in order to match the precision of experimental results and allow consistent comparisons between theory and experiment. We leave this endeavour for future work.   


A caveat that should be mentioned is that the scenarios considered in this work require relatively large scalar quartic couplings -- once again, this comes from the need to generate the correct mass of the Higgs boson, in the absence of any BSM mass term in the Lagrangian. Consequently, under renormalisation-group running, we can expect these couplings to grow rapidly and encounter Landau poles well before the Planck scale -- in a similar way as what occurs in the non-decoupling limit of (non-CSI) extensions of the SM. This is consistent with our expectation that the assumption of CSI at the electroweak scale can a priori not be extended up to the Planck scale, and therefore does not solve the hierarchy problem by itself.


Finally, we should also point out the importance of computing the Higgs trilinear coupling to high precision when studying the behaviour of a number of BSM phenomena occuring in models with extended scalar sectors -- in particular the possibility of a strong-first order electroweak phase transition and the spectrum of gravitational waves that may be produced during the former. In complementarity with collider searches, the future measurement of primordial gravitational waves (GW) at the LISA and DECIGO space-based interferometers will provide a new way to probe the value of the Higgs trilinear coupling and the shape of the Higgs potential. The synergy between the two -- collider and GW measurements -- was studied at one loop for $O(N)$-symmetric CSI models in Ref.~\cite{Hashino:2016rvx} (see also Refs.~\cite{Kakizaki:2015wua,Hashino:2016xoj,Hashino:2018wee} for similar studies in non-CSI theories).
Given the large deviation from the SM in the one-loop predictions for $\lambda_{hhh}$, one expects in CSI models the EWPT to be of strong first order ($c.f.$ Ref~\cite{Kanemura:2004ch}). However, two-loop corrections to the Higgs trilinear coupling certainly affect the strength of the EWPT to some extent, although for the time being this type of effects have only been studied for a non-CSI IDM in Ref.~\cite{Senaha:2018xek} -- the EWPT was then found to be slightly weaker once two-loop corrections to $\hat\lambda_{hhh}$ where included. In turn, this also modifies the spectrum of gravitational waves produced during the EWPT if it is of strong first-order nature. There is therefore strong motivation to consider two-loop corrections to the Higgs trilinear coupling at finite temperatures in CSI theories. 


\section{Summary}
\label{SEC:Conclusion}
In this paper, we have performed the first explicit calculation of leading two-loop corrections to the Higgs trilinear coupling in models with classical scale invariance, using the effective-potential approximation. For the wide range of CSI models in which the scalon/Higgs direction does not mix with other states -- so that we can compute $V_{\text{eff}}$ as in eq.~(\ref{EQ:form_Veff_1+2l}) -- we have derived general two-loop results relating the coefficients of the effective potential to the corrections to the Higgs trilinear coupling. Importantly, we find that while at one loop the prediction for $\hat\lambda_{hhh}$ is the same in all CSI theories, this is not any more the case at two loops. 

We focused our investigations on two particular types of CSI theories: first an $N$-scalar model with a global $O(N)$ symmetry (where we took either $N=1$ or $N=4$ for numerical applications), and next a CSI-2HDM. For both models, we computed the leading corrections at two loops involving the BSM scalars -- and in the CSI-2HDM, the top quark as well. Because these expressions are derived from $V^{(2)}$, which is computed in the {\ensuremath{\overline{\rm MS}}}\xspace scheme, results are obtained at first in the {\ensuremath{\overline{\rm MS}}}\xspace scheme. Therefore, we have included the necessary scheme conversion to rewrite the corrections to the Higgs trilinear coupling in terms of physical parameters (pole masses and physical Higgs VEV) and to take into account the effects from finite WFR. 
Once we turn to the numerical study of these models, we observe that the two-loop corrections clearly lift the degeneracy in $\hat\lambda_{hhh}$ among different BSM masses and parameters. Furthermore, the results for $\hat\lambda_{hhh}$ in CSI and non-CSI variants of the same model strongly differ in their theoretical behaviours, mainly because of the difference at one loop -- constant (mass-independent) value in the CSI case, and corrections proportional to the fourth power of the BSM mass(es) in the non-CSI case. 

However, when we also consider the constraints on the allowed ranges of the BSM parameters coming from the criterion of tree-level perturbative unitarity and from the mass relation in CSI theories, the possible values of $\hat\lambda_{hhh}$ become quite restricted, namely $90\%\lesssim (\delta R)^{(2)}\lesssim 110-115\%$. In other words, while the Higgs trilinear coupling at two loops is not universally predicted in all CSI models (without mixing), its allowed values are still severely limited by the classical scale invariance. Nevertheless, the additional \textit{positive} deviation of order $5-30\%$ of the Higgs trilinear coupling from its SM prediction at two loops can make the BSM discrepancy easier to find at future experiments, and is of the same order of magnitude as the expected accuracy of the determination of $\lambda_{hhh}$ at future colliders (as discussed in the introduction). Therefore, two-loop calculations will in the future prove necessary for CSI theories in order to consistently compare theoretical predictions with experimental results. 

While it may be difficult to distinguish between CSI scenarios and non-CSI scenarios (in the non-decoupling limit) simply from the size of the Higgs trilinear coupling, our work demonstrates the important role of the mass relation that strongly constrains masses and parameters in CSI theories, and may prove to be one of the most powerful tools to probe the CSI nature of a potential BSM discovery. Additionally, as we have mentioned in the previous section, the synergy of the measurement of the spectrum of primordial gravitational waves produced during the EWPT together with that of the Higgs trilinear coupling can also help distinguish scenarios with or without CSI~\cite{Hashino:2016rvx}. 


Finally, in appendix~\ref{APP:generic} we have provided generic {\ensuremath{\overline{\rm MS}}}\xspace results for the coefficients $B^{(2)}$ and $C^{(2)}$ (see eq.~(\ref{EQ:form_Veff_1+2l})) applicable to all CSI models without mixing, and which served to cross-check some of our calculations. We hope these can be of use for the community to study further CSI scenarios. 

\section*{Acknowledgements} 
This work is supported by JSPS, Grant-in-Aid for Scientific Research, No. 16H06492, 18F18022, 18F18321 and 20H00160. This work is also partly supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306.

\vfill 
\pagebreak
