#  On two-loop corrections to the Higgs trilinear coupling in models with
extended scalar sectors

Johannes Braathen<EMAIL_ADDRESS>Department of Physics,
Osaka University, Toyonaka, Osaka 560-0043, Japan  Shinya Kanemura
<EMAIL_ADDRESS>Department of Physics, Osaka University,
Toyonaka, Osaka 560-0043, Japan

###### Abstract

We investigate the possible size of two-loop radiative corrections to the
Higgs trilinear coupling $\lambda_{hhh}$ in two types of models with extended
Higgs sectors, namely in a Two-Higgs-Doublet Model (2HDM) and in the Inert
Doublet Model (IDM). We calculate the leading contributions at two loops
arising from the additional (heavy) scalars and the top quark of these
theories in the effective-potential approximation. We include all necessary
conversion shifts in order to obtain expressions both in the $\overline{\rm
MS}$ and on-shell renormalisation schemes, and in particular, we devise a
consistent “on-shell” prescription for the soft-breaking mass of the 2HDM at
the two-loop level. We illustrate our analytical results with numerical
studies of simple aligned scenarios and show that the two-loop corrections to
$\lambda_{hhh}$ remain smaller than their one-loop counterparts, with a
typical size being $10-20\%$ of the one-loop corrections, at least while
perturbative unitarity conditions are fulfilled. As a consequence, the
existence of a large deviation of the Higgs trilinear coupling from the
prediction in the Standard Model, which has been discussed in the literature
at one loop, is not altered significantly.

††preprint: OU-HET-1001

##  I Introduction

Although the Standard Model (SM) particle spectrum has been completed by the
discovery of a 125-GeV Higgs particle at the CERN LHC Chatrchyan:2012xdj ;
Aad:2012tfa , no sign of any new Physics has been found so far, and direct
searches of non-SM particles are currently putting increasingly stringent
bounds on parameter spaces of Beyond-the-Standard-Model (BSM) theories. At the
moment, the measured properties of the Higgs boson appear to be in close
agreement with their SM predictions, which tends to indicate that new Physics
is either heavy or made difficult to observe by some mechanism such as
alignment Gunion:2002zf  – which is defined as the situation where the Higgs
vacuum expectation value (VEV) is colinear in field space with one (often the
lightest) of the CP-even Higgs mass eigenstates. As a consequence, in aligned
scenarios of BSM models, the coupling constants of the 125-GeV Higgs boson are
equal at tree level to those in the SM, and deviations can only arise via
radiative corrections. However, in models with extended Higgs sectors, some of
the couplings of the SM-like Higgs boson can deviate significantly from the SM
case because of non-decoupling loop effects involving the additional scalar
states of the theory, as was found first in Refs. Kanemura:2002vm ;
Kanemura:2004mg . Among these is the Higgs trilinear coupling $\lambda_{hhh}$,
on which we will focus in this letter.

This coupling is especially important because it participates in the
determination of the shape of the Higgs potential, and in turn the type and
strength of the electroweak phase transition (EWPT). In particular, it has
been shown in Refs. Grojean:2004xa ; Kanemura:2004ch  that large –
$\mathcal{O}(20-30\%)$ or more – deviations in $\lambda_{hhh}$ from its SM
prediction are required for the EWPT to be of strong first order, which is
necessary for the scenario of electroweak baryogenesis (EWBG) Sakharov:1967dj
; Kuzmin:1985mm ; Cohen:1993nk  to be successful.

Current experimental limits on the Higgs trilinear coupling, obtained via
searches for Higgs pair production, are at 95% confidence level
$-5.0<\lambda_{hhh}/\lambda_{hhh}^{\text{SM}}<12.1$ from ATLAS Ferrari:2018akh
(see also Aaboud:2018ftw ) and
$-11<\lambda_{hhh}/\lambda_{hhh}^{\text{SM}}<17$ from CMS Sirunyan:2018iwt
(see also Sirunyan:2018two ). As for further measurement prospects, the HL-LHC
with 3$\text{ ab}^{-1}$ of integrated luminosity could reach DiVita:2017vrr
$0.1<\lambda_{hhh}/\lambda_{hhh}^{\text{SM}}<2.3$, while at the HE-LHC (a
possible 27-TeV upgrade of the LHC) it might be possible to achieve the limits
$0.58<\lambda_{hhh}/\lambda_{hhh}^{\text{SM}}<1.45$ using 15$\text{ ab}^{-1}$
of data Homiller:2018dgu ; Cepeda:2019klc . In the case of lepton colliders,
the ILC operating at 250 GeV cannot access directly the Higgs trilinear
coupling Fujii:2017vwa , but its extension to 500 GeV (1 TeV) could measure
$\lambda_{hhh}$ to a precision of 27% (10%) Fujii:2015jha  using all available
datasets; independently CLIC running at 1.4 and 3 TeV could obtain a result to
$\sim 20\%$ precision (at 68% confidence level) DiVita:2017vrr . Finally, at a
possible 100-TeV hadron collider, one could attain a level of accuracy as good
as $\sim 5-7\%$ when using 30$\text{ ab}^{-1}$ of data Goncalves:2018yva ;
Chang:2018uwu .

Radiative corrections to the Higgs trilinear coupling were first investigated
at the one-loop order in the SM and the minimal supersymmetric SM (MSSM) in
Refs. Barger:1991ed ; Hollik:2001px ; Dobado:2002jz . Moreover, one-loop
effects have also been investigated for various (non-supersymmetric) BSM
theories with extended Higgs sectors – namely with additional singlets
Kanemura:2015fra ; Kanemura:2016lkz ; He:2016sqr ; Kanemura:2017wtm , doublets
Kanemura:2002vm ; Kanemura:2004mg ; Kanemura:2015mxa ; Arhrib:2015hoa ;
Kanemura:2016sos ; Kanemura:2017wtm , or triplets Aoki:2012jj  – and most of
these results are now available in the program H-COUP Kanemura:2017gbi . Since
Refs. Kanemura:2002vm ; Kanemura:2004mg  it is known that, in the non-
decoupling regime, the dominant one-loop BSM corrections to the Higgs
trilinear coupling can cause a deviation of $\lambda_{hhh}$ by several tens of
or even a hundred percent from its SM prediction, while still verifying the
criterion of tree-level unitarity Lee:1977eg . After encountering such large
effects at one loop, one may at first ask whether perturbativity is still
preserved. This is indeed the case because the one-loop expressions are not a
perturbation of the tree-level formula, and instead involve new parameters
that only enter the calculation at loop level. However, it remains natural to
enquire about the situation at two loops: $i.e.$, whether new effects as large
as $\mathcal{O}(100\%)$ may add up with the one-loop ones, or whether
contributions at two loops stay smaller than at one loop.

The first study of leading two-loop corrections in a model exhibiting non-
decoupling effects was performed in Ref. Senaha:2018xek  for the case of the
Inert Doublet Model (IDM), and indicated that two-loop corrections enhance the
Higgs trilinear coupling by a few percent and slightly weaken the first-order
EWPT. We should also mention anterior works in the context of supersymmetric
models, motivated by the need for a consistent theoretical determination of
the Higgs mass(es) and trilinear coupling: these are namely Refs.
Brucherseifer:2013qva  and Muhlleitner:2015dua  where the leading
$\mathcal{O}(\alpha_{t}\alpha_{s})$ SUSY-QCD corrections to $\lambda_{hhh}$
were computed in the MSSM and NMSSM respectively, and their effects were found
to be of the order of 10%.

In this work, we continue along this line of research and investigate the
possible size of two-loop corrections to $\lambda_{hhh}$ both in the IDM and
in an aligned scenario of a Two-Higgs-Doublet Model (2HDM), using the
effective-potential method. For the former, we include new scalar diagrams
that were overlooked in Ref. Senaha:2018xek , while for the latter the
expressions that we obtain constitute the first results in the literature for
the 2HDM and in general for two-loop diagrams involving both heavy Higgs
scalars and top quarks. Moreover, we find the need for a careful treatment of
the renormalisation of the soft-breaking mass $M$ of the 2HDM and we therefore
devise a new prescription ensuring explicitly the decoupling of our
expressions in terms of on-shell-renormalised parameters. We will restrict our
attention here to the two-loop BSM contributions to $\lambda_{hhh}$ from the
additional states in the extended Higgs sectors.

##  II Models

We here briefly recall our conventions for the 2HDM and the IDM. For more
complete reviews of these models, see Refs. Gunion:1989we ; Branco:2011iw  for
the 2HDM and Refs. Barbieri:2006dq ; Goudelis:2013uca  for the IDM.

###  II.1 Two-Higgs-Doublet Model

The first type of model that we consider is a CP-conserving 2HDM, defined in
terms of two $SU(2)_{L}$ doublets $\Phi_{1}$, $\Phi_{2}$ of hypercharge $1/2$.
To avoid flavour changing neutral currents that are strongly constrained
experimentally, we impose a $\mathbb{Z}_{2}$ symmetry under which the two
doublets of the theory transform as $\Phi_{1}\to\Phi_{1}$ and
$\Phi_{2}\to-\Phi_{2}$ Glashow:1976nt , but which is softly broken by a mass
term ($m_{3}^{2}$) in the potential. We follow the conventions of Ref.
Kanemura:2004mg  and write the tree-level scalar potential as

$\displaystyle V^{(0)}_{\text{2HDM}}=$ $\displaystyle\
m_{1}^{2}|\Phi_{1}|^{2}+m_{2}^{2}|\Phi_{2}|^{2}-m_{3}^{2}(\Phi_{2}^{\dagger}\Phi_{1}+\Phi_{1}^{\dagger}\Phi_{2})$
(1)
$\displaystyle+\frac{\lambda_{1}}{2}|\Phi_{1}|^{4}+\frac{\lambda_{2}}{2}|\Phi_{2}|^{4}+\lambda_{3}|\Phi_{1}|^{2}|\Phi_{2}|^{2}+\lambda_{4}|\Phi_{2}^{\dagger}\Phi_{1}|^{2}+\frac{\lambda_{5}}{2}\Big{(}(\Phi_{2}^{\dagger}\Phi_{1})^{2}+\text{h.c.}\Big{)}\,.$

Because we assume that CP is conserved, all mass parameters $m_{i}^{2}$ and
quartic coupling constants $\lambda_{i}$ are real. We choose then to expand
the doublets $\Phi_{1}$ and $\Phi_{2}$ as $\Phi_{i}=(\phi_{i}^{+},\
\phi_{i}^{0}/\sqrt{2})$. Depending on the parameters of the Lagrangian, the
neutral components of the doublets may acquire non-zero (real) vacuum
expectation values (VEVs), which are denoted
$v_{i}\equiv\langle\phi_{i}^{0}\rangle$ and verify
$v_{1}^{2}+v_{2}^{2}=v^{2}$, where $v\simeq 246$ GeV.

Assuming that both $v_{1}$ and $v_{2}$ are non-zero, two dimensionful
parameters – typically $m_{1}^{2}$ and $m_{2}^{2}$ – can be eliminated using
the tadpole equations (see $e.g.$ eqs (9)-(10) in Kanemura:2004mg  for their
tree-level expressions). Seven free parameters then remain in the scalar
sector: $m_{3}^{2}$, $\lambda_{i}$ ($i=1-5$) and the ratio of the two VEVs
$v_{2}/v_{1}=\tan\beta$. The latter defines an angle $\beta$ that diagonalises
the charged and CP-odd Higgs mass matrices at tree-level, while for the CP-
even Higgs mass matrix a second mixing angle $\alpha$ needs to be introduced.
We can then obtain the charged and neutral components of the $SU(2)_{L}$
doublets in terms of mass eigenstates as

$\begin{pmatrix}\phi_{1}^{+}\\\
\phi_{2}^{+}\end{pmatrix}=R_{\beta}\begin{pmatrix}G^{+}\\\
H^{+}\end{pmatrix}\,,\qquad\begin{pmatrix}\phi_{1}^{0}\\\
\phi_{2}^{0}\end{pmatrix}=v\begin{pmatrix}c_{\beta}\\\
s_{\beta}\end{pmatrix}+R_{\alpha}\begin{pmatrix}H\\\
h\end{pmatrix}+iR_{\beta}\begin{pmatrix}G\\\ A\end{pmatrix}\,,$ (2)

with $R_{x}\equiv\left(\begin{smallmatrix}\cos x&-\sin x\\\ \sin x&\cos
x\end{smallmatrix}\right)$, and where $h$ and $H$ are CP-even Higgs bosons,
$A$ is a CP-odd Higgs boson, and $H^{+}$ is a charged Higgs boson. In
addition, $G$ and $G^{+}$ are respectively the neutral and charged Goldstone
bosons associated with electroweak symmetry breaking (EWSB).

It is common to replace the Lagrangian mass parameter $m_{3}^{2}$ by the soft-
breaking mass $M^{2}\equiv 2m_{3}^{2}/s_{2\beta}$, and to express the five
quartic couplings in terms of the four scalar mass eigenvalues and the mixing
angle $\alpha$, using tree-level relations given for example in equations
(26)-(30) of Ref. Kanemura:2004mg . In this respect, it is important to
emphasize that the mass eigenvalues in these equations should be interpreted
as the tree-level ones – otherwise radiative corrections would need to be
taken into account to obtain the relation between quartic couplings and loop-
level mass eigenvalues (see for example Refs. Kanemura:2017wtm ;
Kanemura:2015mxa ; Krause:2016oke ; Basler:2016obg ; Braathen:2017izn ;
Braathen:2017jvs ; Basler:2017uxn ).

To ensure compatibility with experimental constraints, we will throughout this
letter consider the so-called alignment limit. This limit is defined by the
requirement that one of the CP-even Higgs mass eigenstates is aligned in field
space with the full Higgs VEV $v$ Gunion:2002zf . Additionally, we want the
SM-like state to be the lightest eigenstate $h$, which in terms of mixing
angles implies $s_{\beta-\alpha}=1$, or equivalently $\alpha=\beta-\pi/2$. The
tree-level couplings of $h$ to other particles are then equal to their SM
values, and in particular $\lambda_{hhh}^{(0)}=3m_{h}^{2}/v^{2}$.

Furthermore, in the alignment limit and for $m_{h}\ll m_{\Phi}$, we can obtain
simple expressions for the field-dependent tree-level masses of the additional
scalars $\Phi=H,A,H^{\pm}$ as

$m_{\Phi}^{2}(h)\simeq M^{2}+\frac{m_{\Phi}^{2}-M^{2}}{v^{2}}(v+h)^{2}\,.$ (3)

Finally, it should be noted that we neglect throughout this letter
contributions from quarks other than the top and from leptons, so there is no
need to specify the type of fermion couplings in our setting. Indeed, at tree-
level, the couplings of the top quark to the scalar sector are the same in all
types of 2HDMs, and its field-dependent mass is
$m_{t}^{2}(h)=y_{t}^{2}s_{\beta}^{2}(v+h)^{2}/2$.

###  II.2 The Inert Doublet Model

The IDM Deshpande:1977rw ; Barbieri:2006dq  is one of the simplest extensions
of the SM, and corresponds to the limit of the 2HDM in which the previously-
mentioned $\mathbb{Z}_{2}$ symmetry is exact after EWSB. This ensures that
there is no mixing between the SM-like doublet $\Phi_{1}$, and the
$\mathbb{Z}_{2}$-odd one $\Phi_{2}$. Furthermore, it allows the model to
accommodate a dark matter candidate – namely the lightest $\mathbb{Z}_{2}$-odd
scalar. Under the gauge and $\mathbb{Z}_{2}$ symmetries, the scalar potential
is given by

$\displaystyle V^{(0)}_{\text{IDM}}=$
$\displaystyle\mu_{1}^{2}|\Phi_{1}|^{2}+\mu_{2}^{2}|\Phi_{2}|^{2}+\frac{\lambda_{1}}{2}|\Phi_{1}^{2}|^{4}+\frac{\lambda_{2}}{2}|\Phi_{2}^{2}|^{4}+\lambda_{3}|\Phi_{1}|^{2}|\Phi_{2}|^{2}+\lambda_{4}|\Phi_{1}^{\dagger}\Phi_{2}|^{2}+\frac{\lambda_{5}}{2}\left((\Phi_{1}^{\dagger}\Phi_{2})^{2}+\text{h.c.}\right)\,.$
(4)

Following Ref. Aoki:2013lhm , we here decompose the two doublets in terms of
mass eigenstates as

$\Phi_{1}=\begin{pmatrix}G^{+}\\\
\frac{1}{\sqrt{2}}(v+h+iG)\end{pmatrix}\,,\qquad\Phi_{2}=\begin{pmatrix}H^{+}\\\
\frac{1}{\sqrt{2}}(H+iA)\end{pmatrix}\,.$ (5)

where we use the same notations as in the 2HDM. Finally, from the above tree-
level potential, we can derive field-dependent masses for the inert ($i.e.$
$\mathbb{Z}_{2}$-odd) scalars $H$, $A$, and $H^{\pm}$ as
$m_{\Phi}^{2}(h)=\mu_{2}^{2}+\lambda_{\Phi}/2(v+h)^{2}$, where
$\lambda_{H,A}=\lambda_{3}+\lambda_{4}\pm\lambda_{5}$ and
$\lambda_{H^{\pm}}=\lambda_{3}$.

##  III Effective-potential calculation of $\lambda_{hhh}$ at two loops

We investigate leading two-loop corrections to the effective Higgs trilinear
coupling in the effective-potential approximation, which is equivalent to
setting external momenta to zero in a diagrammatic calculation. We define our
loop expansion of the effective potential $V_{\text{eff}}$ as

$V_{\text{eff}}\equiv V^{(0)}+\Delta V_{\text{eff}}=V^{(0)}+\kappa
V^{(1)}+\kappa^{2}V^{(2)}+\cdots\,,$ (6)

where $\kappa\equiv 1/(16\pi^{2})$ is the usual loop factor. While they miss
potential threshold effects (shown at one loop for example in Ref.
Kanemura:2004mg ), effective-potential computations are considerably simpler
than diagrammatic ones and are sufficient for a first study of the magnitude
of two-loop corrections. Furthermore, from past experience with scalar mass
calculations we may expect the inclusion of momentum at two loops to give only
subleading effects – see $e.g.$ Martin:2004kr ; Borowka:2014wla ;
Degrassi:2014pfa ; Braathen:2017izn .

Normalising the effective Higgs trilinear coupling as
$\mathcal{L}\supset-\frac{1}{6}\lambda_{hhh}h^{3}$, the radiative corrections
that it receives can be computed by taking derivatives of the effective
potential as

$\lambda_{hhh}=\frac{\partial^{3}V_{\text{eff}}}{\partial
h^{3}}\bigg{|}_{\text{min}}=\lambda_{hhh}^{(0)}+\Delta\lambda_{hhh}=\lambda_{hhh}^{(0)}+\kappa\delta^{(1)}\lambda_{hhh}+\kappa^{2}\delta^{(2)}\lambda_{hhh}+\cdots\,.$
(7)

In the scenarios without mixing in the scalar sector that we consider in
section IV, the tree-level result $\lambda_{hhh}^{(0)}$ can be simply
expressed in terms of
$[M_{h}^{2}]_{V_{\text{eff}}}=\left[-\frac{1}{v}\frac{\partial}{\partial
h}+\frac{\partial^{2}}{\partial
h^{2}}\right]V_{\text{eff}}\big{|}_{\text{min.}}$, the effective-potential (or
curvature) mass of the lightest Higgs boson, as

$\lambda_{hhh}=\frac{3[M_{h}^{2}]_{V_{\text{eff}}}}{v}+\mathcal{D}_{3}\Delta
V_{\text{eff}}\Big{|}_{\text{min}}\quad\text{where}\quad\mathcal{D}_{3}\equiv\frac{\partial^{3}}{\partial
h^{3}}-\frac{3}{v}\left[-\frac{1}{v}\frac{\partial}{\partial
h}+\frac{\partial^{2}}{\partial h^{2}}\right]\,.$ (8)

The above definition of the differential operator
$\operatorname{\mathcal{D}_{3}}$ ensures that tadpole conditions are taken
into account – $i.e.$ the calculation is performed at the minimum of the loop-
corrected potential.

We follow the common choice of performing renormalisation before taking
derivatives of the potential, which allows us to make use of existing results
for the effective potential Martin:2001vx . The renormalised effective
potential is calculated in terms of field-dependent ($\overline{\rm MS}$)
tree-level masses, and therefore the results we find for $\lambda_{hhh}$ are
also expressed in terms of $\overline{\rm MS}$-renormalised parameters. While
theoretically consistent and simple, $\overline{\rm MS}$-scheme calculations
may be plagued by large logarithmic contributions that appear because of the
explicit renormalisation scale dependence, and furthermore it requires the
inclusion of renormalisation group equations (RGEs) for all running
parameters. Therefore, we choose to use an OS scheme instead and express our
results in terms of physical parameters. For this purpose, we translate the
relevant parameters, $i.e.$ all particle masses and the Higgs VEV, from their
$\overline{\rm MS}$ values
$[M_{h}^{2}]_{V_{\text{eff}}},\,m_{H},\,m_{A},\,m_{H^{\pm}},\,m_{t},\,v$ to OS
ones
$M_{h},\,M_{H},\,M_{A},\,M_{H^{\pm}},\,M_{t},\,v_{\text{phys}}=(\sqrt{2}G_{F})^{-1/2}$
($G_{F}$ being the Fermi constant), and we also include finite wave-function
renormalisation (WFR) as

$\hat{\lambda}_{hhh}=\left(\frac{Z_{h}^{\text{OS}}}{Z_{h}^{\overline{\rm
MS}}}\right)^{3/2}\lambda_{hhh}=\left(1+\frac{3}{2}\delta
Z_{h}^{\text{OS}}-\frac{3}{2}\delta Z_{h}^{\overline{\rm
MS}}\right)\lambda_{hhh}=\bigg{(}1+\frac{3}{2}\frac{d}{dp^{2}}\Pi_{hh}(p^{2})\big{|}_{p^{2}=M_{h}^{2}}\bigg{)}\lambda_{hhh}\,,$
(9)

where $\delta Z_{h}^{\text{OS}}$ and $\delta Z_{h}^{\overline{\rm MS}}$ are
the WFR counterterms in the OS and $\overline{\rm MS}$ schemes respectively,
and $\Pi_{hh}(p^{2})$ is the finite part of the Higgs self-energy evaluated at
external momentum equal to $p^{2}$. We recall that the pole and curvature
masses of the Higgs boson are related as (see $e.g.$ Degrassi:2012ry )

$M_{h}^{2}=[M_{h}^{2}]_{V_{\text{eff}}}+\Pi_{hh}(p^{2}=M_{h}^{2})-\Pi_{hh}(p^{2}=0)\,.$
(10)

As our main concern is the size of the dominant two-loop BSM contributions to
$\lambda_{hhh}$ due to the additional scalar states in the 2HDM and the IDM,
we make the further approximation of neglecting contributions from the 125-GeV
Higgs and the would-be Goldstone bosons, at both one- and two-loop orders,
throughout the following. This will not impact our conclusions on the possible
magnitude of two-loop effects, because these effects are maximal for large
BSM-scalar masses, and this is also exactly when the validity of the
approximation of neglecting the lighter-scalar masses is best. Moreover, as we
work here in aligned scenarios of New Physics, corrections involving only
masses of the Goldstone and 125-GeV-Higgs bosons are common with the SM, and
will therefore drop out of the BSM deviations that we will present in the next
section.

Before turning to the two-loop computation and our new results, we briefly
review here the effective-potential calculation of one-loop corrections to the
Higgs trilinear coupling. The dominant terms in the one-loop effective
potential, for both the 2HDM and the IDM, are Jackiw:1974cv

$V^{(1)}=-3m_{t}^{4}(h)\left(\operatorname{\overline{\text{log}}}m_{t}^{2}(h)-\frac{3}{2}\right)+\sum_{\Phi=H,A,H^{\pm}}\frac{n_{\Phi}m_{\Phi}^{4}(h)}{4}\left(\operatorname{\overline{\text{log}}}m_{\Phi}^{2}(h)-\frac{3}{2}\right)\,,$
(11)

where $m_{t}^{2}(h)$ and $m_{\Phi}^{2}(h)$ are the field-dependent masses of
the top quark and of the extra scalars, respectively, and $n_{\Phi}$ is 1 for
$H$ and $A$, and 2 for $H^{\pm}$ – as mentioned earlier, we have neglected
here the SM-like Higgs and Goldstone boson terms.

One can then derive the leading one-loop contributions to $\lambda_{hhh}$ by
using the operator $\operatorname{\mathcal{D}_{3}}$ Kanemura:2002vm ;
Senaha:2018xek

$\operatorname{\mathcal{D}_{3}}V^{(1)}\Big{|}_{\text{min}}=-\frac{48m_{t}^{4}}{v^{3}}+\sum_{\Phi=H,A,H^{\pm}}\frac{4n_{\Phi}m_{\Phi}^{4}}{v^{3}}\bigg{(}1-\frac{\mathcal{M}^{2}}{m_{\Phi}^{2}}\bigg{)}^{3}\,,$
(12)

where $\mathcal{M}$ is either $M$ in the 2HDM or $\mu_{2}$ in the IDM.

Figure 1: Topologies of diagrams with scalars and fermions contributing to the
effective potential at two loops.

Corrections to the two-loop effective potential are obtained by calculating
one-particle-irreducible vacuum bubble diagrams Jackiw:1974cv . For our study,
we expand the two-loop part of $V_{\text{eff}}$ as
$V^{(2)}=V^{(2)}_{SSS}+V^{(2)}_{SS}+V^{(2)}_{FFS}$, where each index $S$ or
$F$ indicates a scalar or Dirac-fermion propagator – the corresponding
diagrams are shown in figure 1. Analytic expressions for each of these terms
can be obtained in the Landau gauge and $\overline{\rm MS}$ scheme for any
renormalisable model using111Note that our notation differs slightly from that
of Ref. Martin:2001vx  because we work here with Dirac fermions, and not Weyl
fermions, therefore our $V^{(2)}_{FFS}$ corresponds to the sum
$V^{(2)}_{FFS}+V^{(2)}_{\bar{F}\bar{F}S}$ in Martin:2001vx . the results of
Ref. Martin:2001vx  (see also Martin:2018emo  for results with a general gauge
fixing). These involve only two loop functions, namely the one-loop function
$A$ and the two-loop sunrise integral $I$, for both of which complete
expressions are given $e.g.$ in Refs. Ford:1992pn ; Martin:2001vx ;
Martin:2003qz ; Degrassi:2009yq , and useful limits of $I$ with one or more
mass arguments equal or vanishing can be found in Refs. Martin:2003qz ;
Braathen:2016cqe .

At two loops, the $\overline{\rm MS}$ to OS scheme conversion requires adding
finite one-loop or two-loop shifts to the parameters that enter at one loop
and tree level respectively. However for the latter, $i.e.$ the Higgs mass
$[M_{h}^{2}]_{V_{\text{eff}}}$ and VEV $v$, the two-loop shifts yield
corrections to $\hat{\lambda}_{hhh}$ proportional to the 125-GeV Higgs mass
and should hence be neglected in our approximation. Similarly, corrections
involving two-loop WFR are also proportional to the 125-GeV Higgs mass and
thus subleading. To summarise, we only need one-loop scheme translations for
the Higgs VEV, the scalar masses, and the top quark mass, as well as one-loop
finite WFR.

Finally, before considering BSM corrections, we should mention also the case
of the SM calculation, performed, for example, in Ref. Senaha:2018xek .
Starting from the two-loop SM effective potential given in Ref. Ford:1992pn ,
we obtain the same result as equation (11) of Senaha:2018xek  in terms of
$\overline{\rm MS}$ parameters. However, when translating that expression to
the OS scheme, both using the results of Ref. Degrassi:2012ry  as well as with
a standalone calculation, we have

$\delta^{(2)}\hat{\lambda}_{hhh}=\frac{72M_{t}^{4}}{v_{\text{phys}}^{3}}\bigg{(}16g_{3}^{2}-\frac{13M_{t}^{2}}{v_{\text{phys}}^{2}}\bigg{)}\,.$
(13)

Our results do not agree with equation (12) of Senaha:2018xek  as we find for
the numerical coefficient of the two-loop $M_{t}^{6}$ term $-936$ instead of
$336$. We have furthermore checked that the numerical values for our OS
expression and the $\overline{\rm MS}$ one evaluated at renormalisation
$Q=M_{t}$ are in excellent agreement.

##  IV Numerical examples

###  IV.1 An aligned scenario with degenerate heavy scalars in the 2HDM

As our first numerical example, we consider a simplified scenario of the 2HDM
where the additional scalars are degenerate in mass. This ensures that our
calculations contain only three mass scales $\tilde{M}$, $M_{\Phi}$, and
$M_{t}$ and allow relatively compact analytical expressions to be obtained.
Furthermore, to avoid complications arising from mixing between $h$ and $H$,
the CP-even mixing angle $\alpha$ is fixed222Note that, in principle, we
should take into account radiative corrections to the alignment condition, but
as was studied $e.g.$ in Braathen:2017izn  these are typically minute, and we
will neglect these effects here.  as $\alpha=\beta-\pi/2$ to ensure alignment.

In the 2HDM, there are with respect to the SM fifteen new diagrams involving
heavy scalars and top quarks that contribute to the effective potential at two
loops, which we can write as $V^{(2)}_{hHH}$, $V^{(2)}_{hAA}$,
$V^{(2)}_{hH^{\pm}H^{\pm}}$, $V^{(2)}_{HHH}$, $V^{(2)}_{HAA}$,
$V^{(2)}_{HH^{\pm}H^{\pm}}$, $V^{(2)}_{HH}$, $V^{(2)}_{HA}$,
$V^{(2)}_{HH^{\pm}}$, $V^{(2)}_{AA}$, $V^{(2)}_{AH^{\pm}}$,
$V^{(2)}_{H^{\pm}H^{\pm}}$, $V^{(2)}_{ttH}$, $V^{(2)}_{ttA}$, and
$V^{(2)}_{tbH^{\pm}}$. Applying the operator $\operatorname{\mathcal{D}_{3}}$
to these effective-potential terms, we obtain

$\displaystyle\delta^{(2)}\lambda_{hhh}=$ $\displaystyle\
\frac{16m_{\Phi}^{4}}{v^{5}}\left(4+9\cot^{2}2\beta\right)\left(1-\frac{M^{2}}{m_{\Phi}^{2}}\right)^{4}\big{[}-2M^{2}-m_{\Phi}^{2}+(M^{2}+2m_{\Phi}^{2})\operatorname{\overline{\text{log}}}m_{\Phi}^{2}\big{]}$
$\displaystyle+\frac{192m_{\Phi}^{6}\cot^{2}2\beta}{v^{5}}\left(1-\frac{M^{2}}{m_{\Phi}^{2}}\right)^{4}\big{[}1+2\operatorname{\overline{\text{log}}}m_{\Phi}^{2}\big{]}$
$\displaystyle+\frac{96m_{\Phi}^{4}m_{t}^{2}\cot^{2}\beta}{v^{5}}\left(1-\frac{M^{2}}{m_{\Phi}^{2}}\right)^{3}\big{[}-1+2\operatorname{\overline{\text{log}}}m_{\Phi}^{2}\big{]}+\mathcal{O}\left(\frac{m_{\Phi}^{2}m_{t}^{4}}{v^{5}}\right)\,.$
(14)

in terms of the $\overline{\rm MS}$-renormalised parameters – $m_{\Phi}$ being
an $\overline{\rm MS}$-scheme degenerate mass for the heavy scalars
$H,\,A,\,H^{\pm}$ – and with
$\operatorname{\overline{\text{log}}}x\equiv\log(x/Q^{2})$, $Q$ being the
renormalisation scale. The complete expression of the third derivative of the
heavy scalar and top quark sunrise $V^{(2)}_{FFS}$ diagrams is rather long, so
for brevity we only write here the leading
$\mathcal{O}(m_{\Phi}^{4}m_{t}^{2}/v^{5})$ term (while we use the complete
result for our following numerical investigation). We performed some
consistency checks of these results by: $(i)$ verifying that the $\log Q^{2}$
dependence of the total result for $\lambda_{hhh}$ is eliminated when
including the running of all parameters appearing at lower orders; $(ii)$
confirming that each of the terms in this $\overline{\rm MS}$ expression
independently decouples when taking the limit $M\to\infty$. The latter can be
understood as each term is proportional to

$(m_{\Phi}^{2})^{n-1}\left(1-\frac{M^{2}}{m_{\Phi}^{2}}\right)^{n}\underset{m_{\Phi}^{2}=M^{2}+\tilde{\lambda}v^{2}}{=}\frac{(\tilde{\lambda}v^{2})^{n}}{M^{2}+\tilde{\lambda}v^{2}}\,,$
(15)

with $n=3$ or $4$, and where $\tilde{\lambda}$ denotes some combination of
Lagrangian scalar quartic couplings. We should expect to observe decoupling of
the BSM corrections when taking $M\to\infty$ while keeping
$\tilde{\lambda}v^{2}$ finite, so that the additional scalar masses go to
infinity without the calculation entering the non-decoupling regime associated
with large scalar quartic couplings, which would cause perturbativity to be
lost. Here, when taking the limit $M\to\infty$ with $\tilde{\lambda}v^{2}$
fixed our expressions do indeed decouple properly, as can be seen from eq.
(15).

Instead of $\overline{\rm MS}$-renormalised parameters, we prefer to work in
terms of physical parameters, and therefore we now convert our expressions to
the OS scheme. For the masses of the top quark and the heavy scalars, this
simply requires shifting the masses in the one-loop corrections given in eq.
(12) by the corresponding self-energies. At this point, we should emphasise
that the scenarios where the heavy scalars have degenerate $\overline{\rm MS}$
masses or OS masses correspond to different points in the parameter space of
the 2HDM, because the radiative corrections that relate $\overline{\rm MS}$
and OS masses are not the same for the different scalars. Keeping this in
mind, we choose however to consider parameter points for which the scalars
$H$, $A$, and $H^{\pm}$ have a common physical mass $M_{\Phi}$, after the
conversion of our results to the OS scheme. Then, for $\tan\beta$, we do not
need to perform any conversion, as this parameter only enters the calculation
at two loops.

Finally, the treatment of $M^{2}$ is more subtle, as we will discuss now. When
working at one-loop order, one may find decoupling of the effects of the
heavy-scalar loops in the OS scheme result in the limit $M\to\infty$ when
using a relation $M_{\Phi}^{2}=M^{2}+\tilde{\lambda}v^{2}$, with $M_{\Phi}$
renormalised in the OS scheme and $M$ in the $\overline{\rm MS}$ scheme
Kanemura:2004mg . However, when going to two-loop order this is not the case
any more, and one needs to relate parameters that appear at one-loop order in
different schemes with a one-loop equation, so as not to miss two-loop order
effects. We have checked that decoupling does occur if we consistently use a
one-loop relation between $M_{\Phi}$ and $M$ in our results – this is
essentially equivalent to using expressions with the heavy scalar masses
renormalised in the $\overline{\rm MS}$ scheme. Nevertheless, this situation
motivates devising an “on-shell” renormalisation333We write here “on-shell”
with inverted commas, as we are not actually relating $\tilde{M}$ to some
physical observable. condition for the soft-breaking mass, which we then
denote $\tilde{M}$, that would make decoupling apparent when using a relation
of the form $M_{\Phi}^{2}=\tilde{M}^{2}+\tilde{\lambda}v_{\text{phys}}^{2}$.
Furthermore, we still have the freedom to choose this renormalisation
condition in such a way that it ensures the complete cancellation of all
$\operatorname{\overline{\text{log}}}m_{\Phi}^{2}$ terms in
$\delta^{(2)}\lambda_{hhh}$. We emphasise here that $\tilde{M}$ is the OS-
renormalised value of the soft-breaking scale of the $\mathbb{Z}_{2}$ symmetry
of the 2HDM, and that therefore it is the parameter that governs the
possibility of decoupling of the additional states in the extended Higgs
sector.

With this choice, we can derive the finite “OS” counterterm
$\delta^{\text{OS}}M^{2}$ for $M^{2}$ – defined so that
$\tilde{M}^{2}=M^{2}+\delta^{\text{OS}}M^{2}$ – and we obtain at the one-loop
order

$\displaystyle\delta^{\text{OS}}M^{2}=-\kappa\bigg{(}$
$\displaystyle\frac{12(M^{2}-m_{\Phi}^{2})M^{2}\cot^{2}2\beta}{v^{2}}[\operatorname{\overline{\text{log}}}m_{\Phi}^{2}-1]$
$\displaystyle+\frac{3M^{2}m_{t}^{2}\cot^{2}\beta}{v^{2}}\Big{[}B_{0}(m_{\Phi}^{2},m_{t}^{2},m_{t}^{2})+B_{0}(m_{\Phi}^{2},0,m_{t}^{2})\Big{]}\bigg{)}\,,$
(16)

where $B_{0}$ is the (finite part of the) usual Passarino-Veltman function
Passarino:1978jh . Our final OS scheme result for
$\delta^{(2)}\hat{\lambda}_{hhh}$ is then

$\displaystyle\delta^{(2)}\hat{\lambda}_{hhh}=$ $\displaystyle\
\frac{48M_{\Phi}^{6}}{v_{\text{phys}}^{5}}\left(1-\frac{\tilde{M}^{2}}{M_{\Phi}^{2}}\right)^{4}\left\\{4+3\cot^{2}2\beta\left[3-\frac{\pi}{\sqrt{3}}\left(\frac{\tilde{M}^{2}}{M_{\Phi}^{2}}+2\right)\right]\right\\}$
$\displaystyle+\frac{576M_{\Phi}^{6}\cot^{2}2\beta}{v_{\text{phys}}^{5}}\left(1-\frac{\tilde{M}^{2}}{M_{\Phi}^{2}}\right)^{4}+\frac{288M_{\Phi}^{4}M_{t}^{2}\cot^{2}\beta}{v_{\text{phys}}^{5}}\left(1-\frac{\tilde{M}^{2}}{M_{\Phi}^{2}}\right)^{3}$
$\displaystyle+\frac{168M_{\Phi}^{4}M_{t}^{2}}{v_{\text{phys}}^{5}}\left(1-\frac{\tilde{M}^{2}}{M_{\Phi}^{2}}\right)^{3}-\frac{48M_{\Phi}^{6}}{v_{\text{phys}}^{5}}\left(1-\frac{\tilde{M}^{2}}{M_{\Phi}^{2}}\right)^{5}+\mathcal{O}\left(\frac{M_{\Phi}^{2}M_{t}^{4}}{v_{\text{phys}}^{5}}\right)\,,$
(17)

where the terms on the third line come from finite WF and VEV renormalisation.

Figure 2: Illustrations of our results for the deviation $\delta R$ of
$\hat{\lambda}_{hhh}$ computed in the 2HDM with respect to its SM prediction,
with $\delta
R\equiv\Delta\hat{\lambda}_{hhh}^{\text{2HDM}}/\hat{\lambda}_{hhh}^{\text{SM}}=\hat{\lambda}_{hhh}^{\text{2HDM}}/\hat{\lambda}_{hhh}^{\text{SM}}-1$.
(Left side): Decoupling behaviour of the BSM contributions to
$\hat{\lambda}_{hhh}$, shown by plotting $\delta R$ at one loop (solid blue
curve) and two loops (red dot-dashed curve) as a function of $\tilde{M}$. The
degenerate pole mass of the additional scalars $M_{\Phi}$ is taken to be
$M_{\Phi}^{2}=\tilde{M}^{2}+\tilde{\lambda}v_{\text{phys}}^{2}$, with
$\tilde{\lambda}v_{\text{phys}}^{2}=(200\text{ GeV})^{2},\,(300\text{
GeV})^{2},\,(400\text{ GeV})^{2}$, and we fix $\tan\beta=1.5$. (Right side):
$\delta R$ computed at one loop (solid blue curve) and two loops (red dot-
dashed curve) as a function of $M_{\Phi}$, for the maximal non-decoupling case
of $\tilde{M}=0$, and with $\tan\beta=1.1$. The light-red shaded region shows
the values of the additional scalar masses currently excluded by experimental
searches. It should be noted that this constraint on $M_{\Phi}$ can be
weakened for $\tilde{M}>0$.

These two-loop corrections indeed decouple explicitly when taking the limit
$\tilde{M}\to\infty$ with $\tilde{\lambda}v_{\text{phys}}^{2}$ fixed. An
example of this is shown in the left side of figure 2, where we plot the
deviation $\delta
R\equiv\hat{\lambda}_{hhh}^{\text{2HDM}}/\hat{\lambda}_{hhh}^{\text{SM}}-1$ of
$\lambda_{hhh}$ calculated in the 2HDM with respect to the SM prediction as a
function of the OS-renormalised $\tilde{M}$, at one- and two-loop orders, for
different fixed values of $\sqrt{M_{\Phi}^{2}-\tilde{M}^{2}}$ and for
$\tan\beta=1.5$ – other physical inputs being taken from the PDG
Tanabashi:2018oca .

A comment about experimental constraints on the 2HDM parameter space should
also be made at this point. Indeed, the allowed values of the charged Higgs
mass and of $\tan\beta$ can be constrained by searches of charged and neutral
scalars at the LHC as well as by results on flavour observables – some
detailled discussions can be found, for example, in Refs. Misiak:2017bgg ;
Arbey:2017gmh . Furthermore, with the Mathematica package SARAH Staub:2008uz ;
Staub:2009bi ; Staub:2010jh ; Staub:2012pb ; Staub:2013tta ; Staub:2015kfa ,
we have created a SPheno Porod:2003um ; Porod:2011nf  spectrum generator for a
type-I 2HDM444Other types of 2HDMs, in particular type II and type Y, are more
severely constrained by flavour observables – see $e.g.$ Ref. Misiak:2017bgg .
, which allowed us to check with HiggsBounds-5.3.2beta Bechtle:2008jh ;
Bechtle:2011sb ; Bechtle:2013wla ; Bechtle:2015pma  that for $\tilde{M}=0$ and
$\tan\beta=1.1$ ($\tan\beta=1.5$), BSM scalar masses above $350\text{ GeV}$
($355\text{ GeV}$) are not currently excluded by experimental searches.

The other interesting limit to consider with our analytical results is the
case of maximal non-decoupling effects that we obtain for $\tilde{M}=0$. We
illustrate the non-decoupling behaviour of the BSM corrections to
$\hat{\lambda}_{hhh}$ in the right side of figure 2, where we show the same
deviation $\delta R$ as in the left-side plot but now as a function of the
heavy scalar (degenerate) pole mass $M_{\Phi}$, in the case of $\tilde{M}=0$
(to enhance as much as possible the non-decoupling effects) and
$\tan\beta=1.1$. From essentially negligible before the one-loop corrections
for $M_{\Phi}\lesssim v_{\text{phys}}$, the two-loop contributions become as
large as 80% for $M_{\Phi}=500\text{ GeV}$ – the one-loop deviation is then
250%. One should note that the value $M_{\Phi}=500\text{ GeV}$ is close to the
upper limit on $M_{\Phi}$ allowed by the criterion of tree-level unitarity –
using Kanemura:1993hm , we find this limit to be $M_{\Phi}\simeq 600\text{
GeV}$ for $\tan\beta=1.1$ and $\tilde{M}=0$. For $M_{\Phi}=400\text{ GeV}$,
well below the bound from perturbative unitarity, the BSM contributions cause
a deviation of $\hat{\lambda}_{hhh}^{\text{2HDM}}$ with respect to its SM
prediction of 101% at one loop, and of a further 22% at two loops. Finally, we
should mention also the dependence on $\tan\beta$ that appears at two loops in
$\hat{\lambda}_{hhh}$, even in the alignment limit that we have considered. In
particular, the effect of $\tan\beta$ is largest in the scalar contributions
to $\delta^{(2)}\hat{\lambda}_{hhh}$ (the first two terms in eq. (IV.1)),
because of their $\cot^{2}2\beta$ dependence, and hence these terms are
greatly enhanced when $\tan\beta$ increases. However, we observe that the
perturbative expansion is not broken – $i.e.$ two-loop corrections to
$\hat{\lambda}_{hhh}$ remain smaller than their one-loop counterparts – at
least while perturbative unitarity conditions are not violated. To illustrate
this, we consider two example points. For the first one, we take
$M_{\Phi}=400\text{ GeV}$ and $\tilde{M}=0$, and the criterion of tree-level
unitarity then implies Kanemura:1993hm  an upper bound $\tan\beta\lesssim
1.7$. With this maximal value of $\tan\beta$, the one- and two-loop deviations
of $\hat{\lambda}_{hhh}$ from BSM contributions are respectively 101% and 34%.
We fix for the second example point $M_{\Phi}=250\text{ GeV}$ and
$\tilde{M}=0$, which gives the bound $\tan\beta\lesssim 2.8$, and in turn we
obtain for the deviation of $\hat{\lambda}_{hhh}^{\text{2HDM}}$ from
$\hat{\lambda}_{hhh}^{\text{SM}}$ at one- and two-loop orders respectively 15%
and 6%. It therefore appears that, under the criterion of tree-level
unitarity, the two-loop BSM corrections to the Higgs trilinear coupling can
become at most $\mathcal{O}(30-40\%)$ of the one-loop ones.

###  IV.2 A dark-matter-inspired scenario in the IDM

The second case that we consider is a scenario of the IDM, already studied in
Ref. Senaha:2018xek , in which the additional inert scalar $H$ is light
($M_{H}\simeq M_{h}/2\ll M_{A,H^{\pm}}$) and becomes a DM candidate. The
leading two-loop corrections to $\lambda_{hhh}$ are then due to the
pseudoscalar and charged Higgs bosons, and in order to maximise the size of
the radiative corrections, we set the mass parameter $\mu_{2}$ to zero
throughout this section.

In this scenario, there are eight BSM diagrams that give contributions to
$V^{(2)}$, namely $V^{(2)}_{hAA}$, $V^{(2)}_{hH^{\pm}H^{\pm}}$,
$V^{(2)}_{AH^{\pm}G^{\pm}}$, $V^{(2)}_{HAG}$, $V^{(2)}_{HH^{\pm}G^{\pm}}$,
$V^{(2)}_{AA}$, $V^{(2)}_{AH^{\pm}}$, and $V^{(2)}_{H^{\pm}H^{\pm}}$. Only the
first two of these were included in Ref. Senaha:2018xek , and we find
agreement between our results for these and equation (16) of Senaha:2018xek .
Taking into account the other of the above diagrams, we present here for the
first time the complete $\mathcal{O}(M_{\Phi}^{6}/v_{\text{phys}}^{5})$ and
$\mathcal{O}(\lambda_{2}M_{\Phi}^{4}/v_{\text{phys}}^{3})$ contributions to
$\delta^{(2)}\hat{\lambda}_{hhh}$ (here by $M_{\Phi}$ we mean $M_{A}$,
$M_{H^{\pm}}$, or some combination of the two). After conversion to the OS
scheme, and inclusion of finite WF and VEV renormalisation, these read

$\displaystyle\delta^{(2)}\hat{\lambda}_{hhh}=$
$\displaystyle\frac{6\lambda_{2}}{v_{\text{phys}}^{3}}\big{(}3M_{A}^{4}+4M_{A}^{2}M_{H^{\pm}}^{2}+8M_{H^{\pm}}^{4}\big{)}+\frac{60(M_{A}^{6}+2M_{H^{\pm}}^{6})}{v_{\text{phys}}^{5}}+\frac{24(M_{A}^{2}-M_{H^{\pm}}^{2})^{2}(M_{A}^{2}+M_{H^{\pm}}^{2})}{v_{\text{phys}}^{5}}$
$\displaystyle+\frac{24M_{t}^{4}(M_{A}^{2}+2M_{H^{\pm}}^{2})}{v_{\text{phys}}^{5}}+\frac{42M_{t}^{2}(M_{A}^{4}+2M_{H^{\pm}}^{4})}{v_{\text{phys}}^{5}}-\frac{2(M_{A}^{4}+2M_{H^{\pm}}^{4})(M_{A}^{2}+2M_{H^{\pm}}^{2})}{v_{\text{phys}}^{5}}\,.$
(18)

We emphasise that, as $\lambda_{2}$ only appears in the calculation of
$\lambda_{hhh}$ at two-loop order, we do not need here to specify a choice of
renormalisation scheme for it, as opposed to the inert scalar masses, which
appear at one-loop. Moreover, we expect from our result for $\tilde{M}$ in the
previous section that the tree-level ($\overline{\rm MS}$) condition
$\mu_{2}=0$ will also hold in the OS scheme (see eq. (IV.1)). Interestingly,
one may notice that although the inert scalars do not couple directly to the
top quark, terms involving both the top and inert-scalar masses do appear at
two loops through the interplay of one-loop scalar contributions to the Higgs
WFR with the one-loop top quark correction to $\lambda_{hhh}$ (and vice
versa).

At this point, we can also quantify the effects of the new diagrams that we
include here for the first time. When expressed in the OS scheme, the
contributions arising from the two sunrise diagrams $V^{(2)}_{hAA}$ and
$V^{(2)}_{hH^{\pm}H^{\pm}}$ read, respectively,
$48M_{A}^{6}/v_{\text{phys}}^{5}$ and $96M_{H^{\pm}}^{6}/v_{\text{phys}}^{5}$.
We then find that the corrections to $\delta^{(2)}\lambda_{hhh}$ obtained from
$V^{(2)}_{HAG}$ and $V^{(2)}_{HH^{\pm}G^{\pm}}$ amount to 25% of the
aforementioned results, while the terms coming from the last sunrise diagram
$V^{(2)}_{AH^{\pm}G^{\pm}}$ only give significant effects when $M_{A}$ and
$M_{H^{\pm}}$ are not degenerate – $e.g.$ for $M_{A}=2M_{H^{\pm}}$ these terms
result in an additional 25% positive correction to $\delta^{(2)}\lambda_{hhh}$
compared to the result of Ref. Senaha:2018xek . Moreover, the corrections from
the other diagrams $V^{(2)}_{AA}$, $V^{(2)}_{AH^{\pm}}$, and
$V^{(2)}_{H^{\pm}H^{\pm}}$ involving $\lambda_{2}$ can also become large, as
we will see now.

Figure 3: Deviation $\delta R$ of the Higgs trilinear coupling calculated in
the IDM ($\hat{\lambda}_{hhh}^{\text{IDM}}$) with respect to the SM
($\hat{\lambda}_{hhh}^{\text{SM}}$) – $i.e.$ $\delta
R=\hat{\lambda}_{hhh}^{\text{IDM}}/\hat{\lambda}_{hhh}^{\text{SM}}-1$ – as a
function of the degenerate pole masses of the pseudoscalar and charged Higgses
$M_{\Phi}=M_{A}=M_{H^{\pm}}$, at one loop (solid blue curve) and two loops
(dashed and dot-dashed red curves). For the two-loop results, the dashed and
dot-dashed correspond to different values of the inert doublet quartic
coupling $\lambda_{2}$, respectively, $\lambda_{2}=0$ and $6$. We recall that
we take $\mu_{2}=0$ (to maximise the non-decoupling effects), and that we
neglect contributions from $h$ and $H$, as we assume $M_{h},M_{H}\ll
M_{\Phi}$.

We illustrate our numerical results in figure 3, where we plot the deviation
$\delta R$ of the Higgs trilinear coupling $\hat{\lambda}_{hhh}^{\text{IDM}}$
with respect to its SM prediction $\hat{\lambda}_{hhh}^{\text{SM}}$, as a
function of the pole masses of the heavy scalars, which we take to be equal –
$i.e.$ $M_{A}=M_{H^{\pm}}=M_{\Phi}$ – both for convenience and to keep the
$\rho$ parameter close to 1. One should note that this implies that the third
term in the above equation (IV.2), which corresponds to the
$V^{(2)}_{AH^{\pm}G^{\pm}}$ diagram, vanishes. The impact of the non-vanishing
sunrise diagrams – $V^{(2)}_{hAA}$, $V^{(2)}_{hH^{\pm}H^{\pm}}$,
$V^{(2)}_{HAG}$, and $V^{(2)}_{HH^{\pm}G^{\pm}}$ – is given by the difference
between the solid blue (one-loop) and dashed red (two-loop, with
$\lambda_{2}=0$) curves in figure 3. We also try to evaluate the maximal
possible size of the contributions proportional to $\lambda_{2}$ – $i.e.$
coming from $V^{(2)}_{AA}$, $V^{(2)}_{AH^{\pm}}$, $V^{(2)}_{H^{\pm}H^{\pm}}$ –
under the constraint of perturbative unitarity, and for this purpose evaluate
$\delta R$ at two loops for $\lambda_{2}=6$; note that we do not take
$\lambda_{2}$ larger because we find the tree-level unitarity conditions to be
violated Kanemura:1993hm ; Akeroyd:2000wc  for $M_{\Phi}=500\text{ GeV}$,
$\mu_{2}=0$ and $\lambda_{2}\simeq 6.5$. We note also that, in this type of
IDM scenario with $H$ as a DM candidate of mass $M_{H}\simeq M_{h}/2$, the
mass range that we consider here for $M_{\Phi}=M_{A}=M_{H^{\pm}}$ is not
constrained by collider and DM searches Dercks:2018wch .

For $M_{\Phi}=400\text{ GeV}$, $\hat{\lambda}_{hhh}^{\text{IDM}}$ deviates at
one loop by about $76\%$ with respect to its SM prediction, while the sunrise
diagrams give a further enhancement of $\sim 18\%$ and the remaining two-loop
diagrams (with $\lambda_{2}$) another $\sim 21\%$. On the one hand one may
notice that, in relative size, the corrections from the sunrise diagrams grow
significantly faster than the one-loop and two-loop $\lambda_{2}$-dependent
ones, which can be understood from their expression proportional to
$M_{\Phi}^{6}$ as opposed to $M_{\Phi}^{4}$ for the others. On the other hand,
the remaining two-loop diagrams can potentially give large contributions for
low $M_{\Phi}$, if $\lambda_{2}$ is large, because of their large
combinatorial factor, but their relative importance diminishes for increasing
$M_{\Phi}$. To summarise, similarly to what we found in the 2HDM, the two-loop
corrections in the IDM remain smaller than the one-loop ones, meaning that the
perturbative expansion is not breaking down, at least as long as perturbative
unitarity is fulfilled.

Before concluding, we point out that while the coupling $\lambda_{2}$ is quite
difficult to access experimentally as it is the quartic coupling between inert
scalars (see the potential in eq. (4)), a precise measurement of
$\lambda_{hhh}$ could allow us to obtain some information about the value of
$\lambda_{2}$. We can expect this observation to hold also for couplings of
the Higgs boson with other particles, such as for example the $h\gamma\gamma$
or $hZZ$ couplings, because the quartic coupling $\lambda_{2}$ should appear
in internal scalar loops therein as well.

##  V Summary

We have discussed the magnitude of the deviation of the Higgs trilinear
coupling $\lambda_{hhh}$ from its SM prediction at two loops in two models
with extended Higgs sectors, namely the 2HDM – for which we have obtained for
the first time leading two-loop corrections to $\lambda_{hhh}$ – and the IDM –
where we have improved on the existing results of Ref. Senaha:2018xek . We
have performed calculations both in the $\overline{\rm MS}$ and the on-shell
schemes using the effective-potential method. In the cases where it was
possible, we have compared our expressions with existing works in the
literature, explaining the origin of some differences in the SM and IDM. We
also devised a new “on-shell” renormalisation prescription for the soft-
breaking scale $\tilde{M}$ to maintain explicitly the decoupling of the two-
loop corrections for
$M_{\Phi}^{2}=\tilde{M}^{2}+\tilde{\lambda}v_{\text{phys}}^{2}$ with
$\tilde{M}\to\infty$ and fixed $\tilde{\lambda}v_{\text{phys}}^{2}$.

In the two models we studied, we found new dependences of $\lambda_{hhh}$ on
parameters – respectively $\tan\beta$ in the aligned 2HDM and $\lambda_{2}$ in
the IDM – entering the calculation only from two loops, and which may cause
large enhancements of $\delta^{(2)}\hat{\lambda}_{hhh}$. However, we have
shown that, provided one considers parameter points within the region of
parameter space allowed under the criterion of tree-level unitarity, the two-
loop corrections to $\lambda_{hhh}$ do not grow out of control and remain
smaller than the one-loop effects. When expressed in terms of OS-scheme
parameters, our new two-loop contributions (moderately) enhance the non-
decoupling effects appearing at one-loop, but we should emphasise that we do
not obtain new large – $i.e.$ $\mathcal{O}(100\%)$ or so – corrections at two
loops. The typical size that we find for the two-loop corrections in the OS
scheme, up to $\mathcal{O}(\sim 20\%)$ of the one-loop corrections, implies,
on the one hand, that higher-order contributions do not change the existence
of the non-decoupling effects observed from one loop. But on the other hand,
it also means that in the perspective of precise measurements of the Higgs
trilinear coupling, a careful theoretical calculation of $\lambda_{hhh}$ –
including radiative corrections beyond one loop – will be necessary.

Details about the calculations and discussions in this letter, and in
particular complete expressions for a 2HDM, will be shown elsewhere BK:2019xxx
.

###### Acknowledgements.

We are indebted to Kodai Sakurai for numerous useful discussions. We would
also like to thank Benjamin Fuks, Mark Goodsell, Kentarou Mawatari, and
Sebastian Paßehr for helpful discussions. This work is, in part, supported by
Grant-in-Aid for Scientific Research on Innovative Areas, the Ministry of
Education, Culture, Sports, Science and Technology, No. 16H06492 and No.
18H04587, and Grant H2020-MSCA-RISE-2014, No. 645722 (Non Minimal Higgs). This
work is also supported in part by JSPS KAKENHI Grant No. A18F180220.

## References

  * (1) CMS Collaboration, S. Chatrchyan et al., Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC. Phys. Lett. B716 (2012) 30–61, arXiv:1207.7235 [hep-ex].
  * (2) ATLAS Collaboration, G. Aad et al., Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC. Phys. Lett. B716 (2012) 1–29, arXiv:1207.7214 [hep-ex].
  * (3) J. F. Gunion and H. E. Haber, The CP conserving two Higgs doublet model: The Approach to the decoupling limit. Phys. Rev. D67 (2003) 075019, arXiv:hep-ph/0207010 [hep-ph].
  * (4) S. Kanemura, S. Kiyoura, Y. Okada, E. Senaha, and C. P. Yuan, New physics effect on the Higgs selfcoupling. Phys. Lett. B558 (2003) 157–164, arXiv:hep-ph/0211308 [hep-ph].
  * (5) S. Kanemura, Y. Okada, E. Senaha, and C. P. Yuan, Higgs coupling constants as a probe of new physics. Phys. Rev. D70 (2004) 115002, arXiv:hep-ph/0408364 [hep-ph].
  * (6) C. Grojean, G. Servant, and J. D. Wells, First-order electroweak phase transition in the standard model with a low cutoff. Phys. Rev. D71 (2005) 036001, arXiv:hep-ph/0407019 [hep-ph].
  * (7) S. Kanemura, Y. Okada, and E. Senaha, Electroweak baryogenesis and quantum corrections to the triple Higgs boson coupling. Phys. Lett. B606 (2005) 361–366, arXiv:hep-ph/0411354 [hep-ph].
  * (8) A. D. Sakharov, Violation of CP Invariance, C asymmetry, and baryon asymmetry of the universe. Pisma Zh. Eksp. Teor. Fiz. 5 (1967) 32–35. [Usp. Fiz. Nauk161,no.5,61(1991)].
  * (9) V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov, On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe. Phys. Lett. 155B (1985) 36.
  * (10) A. G. Cohen, D. B. Kaplan, and A. E. Nelson, Progress in electroweak baryogenesis. Ann. Rev. Nucl. Part. Sci. 43 (1993) 27–70, arXiv:hep-ph/9302210 [hep-ph].
  * (11) ATLAS Collaboration, A. Ferrari, “Searches for Higgs boson pair production with ATLAS,” in 13th Conference on the Intersections of Particle and Nuclear Physics (CIPANP 2018) Palm Springs, California, USA, May 29-June 3, 2018.  2018\.  arXiv:1809.08870 [hep-ex].
  * (12) ATLAS Collaboration, M. Aaboud et al., Search for Higgs boson pair production in the $\gamma\gamma b\bar{b}$ final state with 13 TeV $pp$ collision data collected by the ATLAS experiment. JHEP 11 (2018) 040, arXiv:1807.04873 [hep-ex].
  * (13) CMS Collaboration, A. M. Sirunyan et al., Search for Higgs boson pair production in the $\gamma\gamma\mathrm{b\overline{b}}$ final state in pp collisions at $\sqrt{s}=$ 13 TeV. Phys. Lett. B788 (2019) 7–36, arXiv:1806.00408 [hep-ex].
  * (14) CMS Collaboration, A. M. Sirunyan et al., Combination of searches for Higgs boson pair production in proton-proton collisions at $\sqrt{s}=$ 13 TeV. Submitted to: Phys. Rev. Lett. (2018) , arXiv:1811.09689 [hep-ex].
  * (15) S. Di Vita, G. Durieux, C. Grojean, J. Gu, Z. Liu, G. Panico, M. Riembau, and T. Vantalon, A global view on the Higgs self-coupling at lepton colliders. JHEP 02 (2018) 178, arXiv:1711.03978 [hep-ph].
  * (16) S. Homiller and P. Meade, Measurement of the Triple Higgs Coupling at a HE-LHC. arXiv:1811.02572 [hep-ph].
  * (17) Physics of the HL-LHC Working Group Collaboration, M. Cepeda et al., Higgs Physics at the HL-LHC and HE-LHC. arXiv:1902.00134 [hep-ph].
  * (18) K. Fujii et al., Physics Case for the 250 GeV Stage of the International Linear Collider. arXiv:1710.07621 [hep-ex].
  * (19) K. Fujii et al., Physics Case for the International Linear Collider. arXiv:1506.05992 [hep-ex].
  * (20) D. Gonçalves, T. Han, F. Kling, T. Plehn, and M. Takeuchi, Higgs boson pair production at future hadron colliders: From kinematics to dynamics. Phys. Rev. D97 (2018) no. 11, 113004, arXiv:1802.04319 [hep-ph].
  * (21) J. Chang, K. Cheung, J. S. Lee, C.-T. Lu, and J. Park, Higgs-boson-pair production $H(\rightarrow b\overline{b})H(\rightarrow\gamma\gamma)$ from gluon fusion at the HL-LHC and HL-100 TeV hadron collider. arXiv:1804.07130 [hep-ph].
  * (22) V. D. Barger, M. S. Berger, A. L. Stange, and R. J. N. Phillips, Supersymmetric Higgs boson hadroproduction and decays including radiative corrections. Phys. Rev. D45 (1992) 4128–4147.
  * (23) W. Hollik and S. Penaranda, Yukawa coupling quantum corrections to the selfcouplings of the lightest MSSM Higgs boson. Eur. Phys. J. C23 (2002) 163–172, arXiv:hep-ph/0108245 [hep-ph].
  * (24) A. Dobado, M. J. Herrero, W. Hollik, and S. Penaranda, Selfinteractions of the lightest MSSM Higgs boson in the large pseudoscalar mass limit. Phys. Rev. D66 (2002) 095016, arXiv:hep-ph/0208014 [hep-ph].
  * (25) S. Kanemura, M. Kikuchi, and K. Yagyu, Radiative corrections to the Higgs boson couplings in the model with an additional real singlet scalar field. Nucl. Phys. B907 (2016) 286–322, arXiv:1511.06211 [hep-ph].
  * (26) S. Kanemura, M. Kikuchi, and K. Yagyu, One-loop corrections to the Higgs self-couplings in the singlet extension. Nucl. Phys. B917 (2017) 154–177, arXiv:1608.01582 [hep-ph].
  * (27) S.-P. He and S.-h. Zhu, One-Loop Radiative Correction to the Triple Higgs Coupling in the Higgs Singlet Model. Phys. Lett. B764 (2017) 31–37, arXiv:1607.04497 [hep-ph].
  * (28) S. Kanemura, M. Kikuchi, K. Sakurai, and K. Yagyu, Gauge invariant one-loop corrections to Higgs boson couplings in non-minimal Higgs models. Phys. Rev. D96 (2017) no. 3, 035014, arXiv:1705.05399 [hep-ph].
  * (29) S. Kanemura, M. Kikuchi, and K. Yagyu, Fingerprinting the extended Higgs sector using one-loop corrected Higgs boson couplings and future precision measurements. Nucl. Phys. B896 (2015) 80–137, arXiv:1502.07716 [hep-ph].
  * (30) A. Arhrib, R. Benbrik, J. El Falaki, and A. Jueid, Radiative corrections to the Triple Higgs Coupling in the Inert Higgs Doublet Model. JHEP 12 (2015) 007, arXiv:1507.03630 [hep-ph].
  * (31) S. Kanemura, M. Kikuchi, and K. Sakurai, Testing the dark matter scenario in the inert doublet model by future precision measurements of the Higgs boson couplings. Phys. Rev. D94 (2016) no. 11, 115011, arXiv:1605.08520 [hep-ph].
  * (32) M. Aoki, S. Kanemura, M. Kikuchi, and K. Yagyu, Radiative corrections to the Higgs boson couplings in the triplet model. Phys. Rev. D87 (2013) no. 1, 015012, arXiv:1211.6029 [hep-ph].
  * (33) S. Kanemura, M. Kikuchi, K. Sakurai, and K. Yagyu, H-COUP: a program for one-loop corrected Higgs boson couplings in non-minimal Higgs sectors. Comput. Phys. Commun. 233 (2018) 134–144, arXiv:1710.04603 [hep-ph].
  * (34) B. W. Lee, C. Quigg, and H. B. Thacker, Weak Interactions at Very High-Energies: The Role of the Higgs Boson Mass. Phys. Rev. D16 (1977) 1519.
  * (35) E. Senaha, Radiative Corrections to Triple Higgs Coupling and Electroweak Phase Transition: Beyond One-loop Analysis. arXiv:1811.00336 [hep-ph].
  * (36) M. Brucherseifer, R. Gavin, and M. Spira, Minimal supersymmetric Higgs boson self-couplings: Two-loop $O(\alpha_{t}\alpha_{s})$ corrections. Phys. Rev. D90 (2014) no. 11, 117701, arXiv:1309.3140 [hep-ph].
  * (37) M. Mühlleitner, D. T. Nhung, and H. Ziesche, The order $\mathcal{O}\left({\alpha}_{t}{\alpha}_{s}\right)$ corrections to the trilinear Higgs self-couplings in the complex NMSSM. JHEP 12 (2015) 034, arXiv:1506.03321 [hep-ph].
  * (38) J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson, The Higgs Hunter’s Guide. Front. Phys. 80 (2000) 1–404.
  * (39) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, Theory and phenomenology of two-Higgs-doublet models. Phys. Rept. 516 (2012) 1–102, arXiv:1106.0034 [hep-ph].
  * (40) R. Barbieri, L. J. Hall, and V. S. Rychkov, Improved naturalness with a heavy Higgs: An Alternative road to LHC physics. Phys. Rev. D74 (2006) 015007, arXiv:hep-ph/0603188 [hep-ph].
  * (41) A. Goudelis, B. Herrmann, and O. Stål, Dark matter in the Inert Doublet Model after the discovery of a Higgs-like boson at the LHC. JHEP 09 (2013) 106, arXiv:1303.3010 [hep-ph].
  * (42) S. L. Glashow and S. Weinberg, Natural Conservation Laws for Neutral Currents. Phys. Rev. D15 (1977) 1958.
  * (43) M. Krause, R. Lorenz, M. Muhlleitner, R. Santos, and H. Ziesche, Gauge-independent Renormalization of the 2-Higgs-Doublet Model. JHEP 09 (2016) 143, arXiv:1605.04853 [hep-ph].
  * (44) P. Basler, M. Krause, M. Muhlleitner, J. Wittbrodt, and A. Wlotzka, Strong First Order Electroweak Phase Transition in the CP-Conserving 2HDM Revisited. JHEP 02 (2017) 121, arXiv:1612.04086 [hep-ph].
  * (45) J. Braathen, M. D. Goodsell, and F. Staub, Supersymmetric and non-supersymmetric models without catastrophic Goldstone bosons. Eur. Phys. J. C77 (2017) no. 11, 757, arXiv:1706.05372 [hep-ph].
  * (46) J. Braathen, M. D. Goodsell, M. E. Krauss, T. Opferkuch, and F. Staub, $N$-loop running should be combined with $N$-loop matching. Phys. Rev. D97 (2018) no. 1, 015011, arXiv:1711.08460 [hep-ph].
  * (47) P. Basler, M. Mühlleitner, and J. Wittbrodt, The CP-Violating 2HDM in Light of a Strong First Order Electroweak Phase Transition and Implications for Higgs Pair Production. JHEP 03 (2018) 061, arXiv:1711.04097 [hep-ph].
  * (48) N. G. Deshpande and E. Ma, Pattern of Symmetry Breaking with Two Higgs Doublets. Phys. Rev. D18 (1978) 2574.
  * (49) M. Aoki, S. Kanemura, and H. Yokoya, Reconstruction of Inert Doublet Scalars at the International Linear Collider. Phys. Lett. B725 (2013) 302–309, arXiv:1303.6191 [hep-ph].
  * (50) S. P. Martin, Strong and Yukawa two-loop contributions to Higgs scalar boson self-energies and pole masses in supersymmetry. Phys. Rev. D71 (2005) 016012, arXiv:hep-ph/0405022 [hep-ph].
  * (51) S. Borowka, T. Hahn, S. Heinemeyer, G. Heinrich, and W. Hollik, Momentum-dependent two-loop QCD corrections to the neutral Higgs-boson masses in the MSSM. Eur. Phys. J. C74 (2014) no. 8, 2994, arXiv:1404.7074 [hep-ph].
  * (52) G. Degrassi, S. Di Vita, and P. Slavich, Two-loop QCD corrections to the MSSM Higgs masses beyond the effective-potential approximation. Eur. Phys. J. C75 (2015) no. 2, 61, arXiv:1410.3432 [hep-ph].
  * (53) S. P. Martin, Two loop effective potential for a general renormalizable theory and softly broken supersymmetry. Phys. Rev. D65 (2002) 116003, arXiv:hep-ph/0111209 [hep-ph].
  * (54) G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori, and A. Strumia, Higgs mass and vacuum stability in the Standard Model at NNLO. JHEP 08 (2012) 098, arXiv:1205.6497 [hep-ph].
  * (55) R. Jackiw, Functional evaluation of the effective potential. Phys. Rev. D9 (1974) 1686.
  * (56) S. P. Martin and H. H. Patel, Two-loop effective potential for generalized gauge fixing. Phys. Rev. D98 (2018) no. 7, 076008, arXiv:1808.07615 [hep-ph].
  * (57) C. Ford, I. Jack, and D. R. T. Jones, The Standard model effective potential at two loops. Nucl. Phys. B387 (1992) 373–390, arXiv:hep-ph/0111190 [hep-ph]. [Erratum: Nucl. Phys.B504,551(1997)].
  * (58) S. P. Martin, Evaluation of two loop selfenergy basis integrals using differential equations. Phys. Rev. D68 (2003) 075002, arXiv:hep-ph/0307101 [hep-ph].
  * (59) G. Degrassi and P. Slavich, On the radiative corrections to the neutral Higgs boson masses in the NMSSM. Nucl. Phys. B825 (2010) 119–150, arXiv:0907.4682 [hep-ph].
  * (60) J. Braathen and M. D. Goodsell, Avoiding the Goldstone Boson Catastrophe in general renormalisable field theories at two loops. JHEP 12 (2016) 056, arXiv:1609.06977 [hep-ph].
  * (61) G. Passarino and M. J. G. Veltman, One Loop Corrections for e+ e- Annihilation Into mu+ mu- in the Weinberg Model. Nucl. Phys. B160 (1979) 151–207.
  * (62) Particle Data Group Collaboration, M. Tanabashi et al., Review of Particle Physics. Phys. Rev. D98 (2018) no. 3, 030001.
  * (63) M. Misiak and M. Steinhauser, Weak radiative decays of the B meson and bounds on $M_{H^{\pm}}$ in the Two-Higgs-Doublet Model. Eur. Phys. J. C77 (2017) no. 3, 201, arXiv:1702.04571 [hep-ph].
  * (64) A. Arbey, F. Mahmoudi, O. Stal, and T. Stefaniak, Status of the Charged Higgs Boson in Two Higgs Doublet Models. Eur. Phys. J. C78 (2018) no. 3, 182, arXiv:1706.07414 [hep-ph].
  * (65) F. Staub, SARAH. arXiv:0806.0538 [hep-ph].
  * (66) F. Staub, From Superpotential to Model Files for FeynArts and CalcHep/CompHep. Comput. Phys. Commun. 181 (2010) 1077–1086, arXiv:0909.2863 [hep-ph].
  * (67) F. Staub, Automatic Calculation of supersymmetric Renormalization Group Equations and Self Energies. Comput. Phys. Commun. 182 (2011) 808–833, arXiv:1002.0840 [hep-ph].
  * (68) F. Staub, SARAH 3.2: Dirac Gauginos, UFO output, and more. Comput. Phys. Commun. 184 (2013) 1792–1809, arXiv:1207.0906 [hep-ph].
  * (69) F. Staub, SARAH 4 : A tool for (not only SUSY) model builders. Comput. Phys. Commun. 185 (2014) 1773–1790, arXiv:1309.7223 [hep-ph].
  * (70) F. Staub, Exploring new models in all detail with SARAH. Adv. High Energy Phys. 2015 (2015) 840780, arXiv:1503.04200 [hep-ph].
  * (71) W. Porod, SPheno, a program for calculating supersymmetric spectra, SUSY particle decays and SUSY particle production at e+ e- colliders. Comput. Phys. Commun. 153 (2003) 275–315, arXiv:hep-ph/0301101 [hep-ph].
  * (72) W. Porod and F. Staub, SPheno 3.1: Extensions including flavour, CP-phases and models beyond the MSSM. Comput. Phys. Commun. 183 (2012) 2458–2469, arXiv:1104.1573 [hep-ph].
  * (73) P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein, and K. E. Williams, HiggsBounds: Confronting Arbitrary Higgs Sectors with Exclusion Bounds from LEP and the Tevatron. Comput. Phys. Commun. 181 (2010) 138–167, arXiv:0811.4169 [hep-ph].
  * (74) P. Bechtle, O. Brein, S. Heinemeyer, G. Weiglein, and K. E. Williams, HiggsBounds 2.0.0: Confronting Neutral and Charged Higgs Sector Predictions with Exclusion Bounds from LEP and the Tevatron. Comput. Phys. Commun. 182 (2011) 2605–2631, arXiv:1102.1898 [hep-ph].
  * (75) P. Bechtle, O. Brein, S. Heinemeyer, O. Stål, T. Stefaniak, G. Weiglein, and K. E. Williams, $\mathsf{HiggsBounds}-4$: Improved Tests of Extended Higgs Sectors against Exclusion Bounds from LEP, the Tevatron and the LHC. Eur. Phys. J. C74 (2014) no. 3, 2693, arXiv:1311.0055 [hep-ph].
  * (76) P. Bechtle, S. Heinemeyer, O. Stal, T. Stefaniak, and G. Weiglein, Applying Exclusion Likelihoods from LHC Searches to Extended Higgs Sectors. Eur. Phys. J. C75 (2015) no. 9, 421, arXiv:1507.06706 [hep-ph].
  * (77) S. Kanemura, T. Kubota, and E. Takasugi, Lee-Quigg-Thacker bounds for Higgs boson masses in a two doublet model. Phys. Lett. B313 (1993) 155–160, arXiv:hep-ph/9303263 [hep-ph].
  * (78) A. G. Akeroyd, A. Arhrib, and E.-M. Naimi, Note on tree level unitarity in the general two Higgs doublet model. Phys. Lett. B490 (2000) 119–124, arXiv:hep-ph/0006035 [hep-ph].
  * (79) D. Dercks and T. Robens, Constraining the Inert Doublet Model using Vector Boson Fusion. arXiv:1812.07913 [hep-ph].
  * (80) J. Braathen and S. Kanemura, in preparation.


also at ]Jawaharlal Nehru Centre For Advanced Scientific Research, Jakkur,
Bangalore, India

#  The formation of compact objects at finite temperatures in a dark-matter-
candidate self-gravitating bosonic system

Akhilesh Kumar Verma<EMAIL_ADDRESS>Centre for Condensed Matter Theory,
Department of Physics, Indian Institute of Science, Bangalore 560012, India
Rahul Pandit<EMAIL_ADDRESS>[  Centre for Condensed Matter Theory,
Department of Physics, Indian Institute of Science, Bangalore 560012, India
Marc E. Brachet<EMAIL_ADDRESS>Laboratoire de Physique de l’École
Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université Université
de Paris, F-75005 Paris, France

###### Abstract

We study self-gravitating bosonic systems, candidates for dark-matter halos,
by carrying out a suite of direct numerical simulations (DNSs) designed to
investigate the formation of finite-temperature, compact objects in the three-
dimensional (3D) Fourier-truncated Gross-Pitaevskii-Poisson equation (GPPE).
This truncation allows us to explore the collapse and fluctuations of compact
objects, which form at both zero temperature and finite temperature. We show
that the statistically steady state of the GPPE, in the large-time limit and
for the system sizes we study, can also be obtained efficiently by tuning the
temperature in an auxiliary stochastic Ginzburg-Landau-Poisson equation
(SGLPE). We show that, over a wide range of model parameters, this system
undergoes a thermally driven first-order transition from a collapsed, compact,
Bose-Einstein condensate (BEC) to a tenuous Bose gas without condensation. By
a suitable choice of initial conditions in the GPPE, we also obtain a binary
condensate that comprises a pair of collapsed objects rotating around their
center of mass.

Superfluidity; Quantum fluids; Quantum vortices

Gravitational effects are important on stellar scales; it might also be
possible to mimic such effects in laboratory Bose-Einstein condensates
ODell2000  and thus emulate gravitationally bound, condensed, assemblies of
bosons, which are candidates for dark-matter halos RuffiniBonazzola1969 ;
Ingrosso:1991 ; Jetzer1992 ; Hu2000 . Although many experiments have been
carried out to establish the identity of dark matter, there is still no
unambiguous dark-matter candidate. For many years the front runners have been
weakly interacting massive particles (WIMPs); but their existence has not been
established convincingly (see, e.g., abdallah2018search ; kachulis2018search
and references therein), so investigations of other candidates, e.g., axions,
boson stars, black holes, and superfluids, have experienced a renaissance
apslink ; lawson2019tunable ; Chavanis1 ; Chavanis2 ; Chavanis3 ; Chavanis4 ;
Chavanis5 .

The Gross-Pitaevskii-Poisson equation (GPPE), for a self-gravitating assembly
of weakly interacting bosons, is the natural theoretical model for dark-
matter-candidate bosonic assemblies, both in laboratory and astrophysical
settings. We address finite-temperature ($T>0$) effects in such condensation,
a problem of central importance in this challenging field. Several ways have
been suggested for including $T>0$ effects in the Gross-Pitaevskii (GP) model
PhDTutorial  without gravity; one important way uses the Fourier-truncated GP
model, in which this truncation generates a classical-field model PhDTutorial
; Berloff14 ; Krstulovic11 ; vmrnjp13 . We generalize these studies by using
the Fourier-truncated GPPE to study $T>0$ effects during the gravitational
collapse of a system of self-gravitating bosons. In addition, we define an
algorithm that directly reconstructs the thermalized state of such a system by
using an auxiliary stochastic Ginzburg-Landau-Poisson equation (SGLPE).

We obtain several interesting results by pseudospectral direct numerical
simulations (DNSs) of the truncated GPPE and the SGLPE in three dimensions
(3D). We follow the spatiotemporal evolution of different initial conditions
for the density of bosons. If we start from a very-nearly uniform density in
the truncated GPPE, this system undergoes gravitational collapse and
thermalizes to a Bose-Einstein condensate (BEC), at low $T$. Our SGLPE study
yields a hitherto unanticipated, thermally driven, first-order phase
transition from a collapsed, compact, BEC to a tenuous Bose gas without
condensation, for a wide variety of parameters in the GPPE. Finally, by a
suitable choice of initial conditions in the GPPE, we obtain a binary
condensate that comprises a pair of collapsed objects rotating around their
center of mass.

Figure 1:  Columns $(1)-(3)$: Ten-level contour plots of ${\left|\psi({\bf
x},t)\right|}^{2}$, at representative times: SGLPE ($T=0$) (top row, run
$R1$); SGLPE (second row, run $R2$); GPPE (third row, run $R3$); GPPE (fourth
row, run $R4$); and $256^{3}$ GPPE (fifth row, run $R5$) [the videos V1-V5
(Supplemental Material supplement ) show, respectively, the complete
spatiotemporal evolution for these cases]. Column$(4)$: Plots of the scaled
radius of gyration
$\frac{R}{L}=\frac{1}{L}\sqrt{\frac{\int_{V}\rho(r)r^{2}d{\bf
r}}{\int_{V}\rho(r)d{\bf r}}}$, (blue) and the scaled gravitational energy
$E_{\rm grav}/E_{a}$ (red ) versus the scaled time $t/(\xi/v)$ for the
different runs, where $E_{\rm a}=2^{5}\pi^{4}(G/g^{3})^{1/2}$.

A self-gravitating BEC is described by a complex wave function $\psi({\bf
x},t)$; for weakly interacting bosons, the spatiotemporal evolution of
$\psi({\bf x},t)$ is governed by the GPPE:

$\displaystyle i\hbar\partial_{t}\psi$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{2m}\nabla^{2}\psi+\left[G\Phi+g|\psi|^{2}\right]\psi,$
$\displaystyle\nabla^{2}\Phi$ $\displaystyle=$
$\displaystyle|\psi|^{2}-<|\psi|^{2}>,$ (1)

where $m$ is the mass of the bosons, $n=|\psi|^{2}$ their number density,
$\Phi$ is the gravitational potential field, $G=4\pi G_{\rm N}m^{2}$ ($G_{\rm
N}$ is Newton’s constant), and $g=4\pi a\hbar^{2}/m$, with $a$ the $s$-wave
scattering length. The subtraction of the mean density
$\langle|\psi|^{2}\rangle$ can be justified either by taking into account the
cosmological expansion Peebles ; falco2013  or by defining a Newtonian
cosmological constant KIESSLING2003 . By linearizing Eq. (1) around the
constant $|\psi|^{2}=n_{0}$, we obtain the dispersion relation
$\omega(k)=\sqrt{-Gn_{0}/m+k^{2}gn_{0}/m+k^{4}(\hbar/2m)^{2}}$, which displays
a low-$k$ Jeans instability, for wave numbers $k<k_{\rm
J}=\sqrt{\frac{G}{g}}\left[\left(1+\sqrt{1+\frac{G\hbar^{2}}{mg^{2}n_{0}}}\right)/2\right]^{-1/2}$.
In the absence of gravity ($G=0$), we identify the speed of sound
$c=\sqrt{\frac{gn_{0}}{m}}$ and the coherence length
$\xi=\sqrt{\frac{\hbar^{2}}{2g{n_{0}}m}}$. Units relevant to astrophysics are
discussed in Chavanis5 .

We solve the GPPE (1) by using the $3D$ Fourier pseudospectral method, with
the $2/3$-rule for dealiasing Krstulovic11 ; vmrnjp13 ; Got-Ors : We expand
the $2\pi$ periodic wave function as $\psi(x)=\sum_{{\bf
k}\in\mathbb{Z}^{3}}\hat{\psi}_{\bf k}\exp(i{\bf k}\cdot{\bf x})$ and then we
truncate it spectrally, by setting $\hat{\psi}_{\bf k}\equiv 0$ for $|{\bf
k}|>k_{\rm max}$, with $k_{\rm max}=[{\cal N}/3]$, where $\cal{N}$ is the
resolution and $[\cdot]$ denotes the integer part PhDTutorial ; Berloff14 . If
we introduce the Galerkin projector $\mathcal{P}_{\rm G}$ [in Fourier space
$\mathcal{P}_{\rm G}[\hat{\psi}_{\bf k}]=\theta(k_{\rm max}-|{\bf
k}|)\hat{\psi}_{\bf k}$, with $\theta(\cdot)$ the Heaviside function], the
Fourier-truncated GPPE becomes

$i\hbar\frac{\partial\psi}{\partial t}=\mathcal{P}_{\rm
G}[-\frac{\hbar^{2}}{2m}{\bf\nabla}^{2}\psi+\mathcal{P}_{\rm
G}[(G{\bf\nabla}^{-2}+g)|\psi|^{2}]\psi].$ (2)

Equation (2) conserves exactly the number of particles $N=\int
d^{3}x|\psi|^{2}$ and the energy $E=E_{kq}+E_{int}+E_{G}$, where
$E_{kq}=\frac{\hbar^{2}}{2m}\int d^{3}x|{\bf\nabla}\psi|^{2}$,
$E_{int}=\frac{g}{2}\int d^{3}x[\mathcal{P}_{\rm G}|\psi|^{2}]^{2}$, and
$E_{G}=E_{\rm grav}=\frac{G}{2}\int d^{3}x[\mathcal{P}_{\rm
G}|\psi|^{2}]{\bf\nabla}^{-2}[\mathcal{P}_{\rm G}|\psi|^{2}]$. If we use the
$2/3$-rule for dealiasing, then the momentum ${\bf P}=\frac{i\hbar}{2m}\int
d^{3}x\left(\psi{\bf\nabla}\overline{\psi}-\overline{\psi}{\bf\nabla}\psi\right)$,
where the over-bar denotes complex conjugation, is also conserved Krstulovic11
.

This spectral truncation generates a classical-field model that allows us to
study finite-$T$ effects in the GPPE (Refs. Berloff14 ; Krstulovic11 ;
vmrnjp13  for the GP case). We show that this spectrally truncated GPPE can
describe dynamical effects and, at the same time, yield thermalized states,
which we can obtain both by the thermalization of the long-time GPPE dynamics
or directly, by using the SGLPE

$\displaystyle\hbar\frac{\partial\psi}{\partial t}$ $\displaystyle=$
$\displaystyle\mathcal{P}_{\rm
G}[\frac{\hbar^{2}}{2m}{\bf\nabla}^{2}\psi+\mu\psi-\mathcal{P}_{\rm
G}[(G{\bf\nabla}^{-2}+g)|\psi|^{2}]\psi]$ (3)
$\displaystyle+\sqrt{\frac{2\hbar}{\beta}}\mathcal{P}_{\rm G}\left[\zeta({\bf
x},t)\right],\hskip 8.53581pt$

where the zero-mean, Gaussian white noise $\zeta({\bf x},t)$ has the variance
$\langle\zeta({\bf x},t)\zeta^{*}({\bf
x^{\prime}},t^{\prime})\rangle=\delta(t-t^{\prime})\delta({\bf x}-{\bf
x^{\prime}})$, with $\beta=\frac{1}{k_{B}T}$ the inverse temperature; we tune
the chemical potential $\mu$ at each time step to conserve the total number of
particles $N$. The finite-$T$ SGLPE dynamics does not describe any physical
evolution; but it converges more rapidly, than does the GPPE dynamics, toward
a thermalized state (Refs. Berloff14 ; Krstulovic11 ; vmrnjp13  for the GP
case). Clearly, the SGLPE leads to a state with a given temperature; but the
GPPE yields a state with a given energy.

Our pseudospectral DNSs of the GPPE (1) and SGLPE (3) use a cubical
computational domain that is $(2\pi)^{3}$ periodic; we normalise $\psi$ such
that $N=(2\pi)^{3}$, and use units with $\hbar=1$ and $m=1$. We list the
parameters for different runs in Tables 1 and 2.

$M_{\rm a}=\frac{\hbar}{\sqrt{G_{N}ma}}$ | $\frac{M}{M_{\rm a}}=\frac{M\sqrt{Gg}}{4\pi\hbar^{2}}$
---|---
$R_{\rm a}=\sqrt{\frac{a\hbar^{2}}{m^{3}G_{N}}}$ | $\frac{R}{R_{\rm a}}=\frac{R}{\sqrt{g/G}}$
$a_{\rm Q}=\frac{\hbar^{2}}{G_{N}M^{2}m}$ | $\frac{a}{a_{\rm Q}}=\frac{M^{2}{Gg}}{(4\pi\hbar)^{2}}$
$R_{\rm Q}=\frac{\hbar^{2}}{G_{N}m^{2}M}$ | $\frac{R}{R_{\rm Q}}=\frac{M{GR}}{4\pi\hbar^{2}}$

Table 1:  Dimensionless variables: in the second column of the table we use
quantities from the first column; $G_{\rm N}$ denotes Newton’s constant and
$G=4\pi G_{\rm N}m^{2}$.

Run | R1 | R2 | R3 | R4 | R5
---|---|---|---|---|---
Type | SGLPE ($T=0$) | SGLPE | GPPE | GPPE | GPPE
$\cal{N}$ | $64$ | $64$ | $64$ | $64$ | $256$
$g$ | $50$ | $50$ | $50$ | $5$ | $50$
$G$ | $105$ | $105$ | $105$ | $550$ | $105$

Table 2: This table shows the representative runs for which we give plots in
Fig.$1$.

In columns $(1)-(3)$ of Fig. 1 we show ten-level contour plots of
${\left|\psi({\bf x},t)\right|}^{2}$ to illustrate, at representative times,
the spatial organization of ${\left|\psi({\bf x},t)\right|}^{2}$ that we
obatin via the SGLPE ($T=0$) (top row, run $R1$), the SGLPE (second row, run
$R2$), and three GPPE runs (third row, run $R3$; fourth row, run $R4$; fifth
row, run $R5$); the videos V1-V5 (Supplemental Material supplement ) show,
respectively, the complete spatiotemporal evolution of ${\left|\psi({\bf
x},t)\right|}^{2}$ for these five runs. In column $(4)$ of Fig. 1 we give, for
these runs, plots of the scaled radius of gyration
$R/L=\frac{1}{L}\sqrt{\frac{\int_{V}\rho(r)r^{2}d{\bf r}}{\int_{V}\rho(r)d{\bf
r}}}$ (blue curves) and the scaled gravitational energy $E_{\rm grav}/E_{a}$
(red curves) versus the scaled time $t/(\xi/v)$, where $E_{\rm
a}=2^{5}\pi^{4}(G/g^{3})^{1/2}$. If we tune $T$ in the SGLPE (3), it yields a
statistically steady state whose properties (like $R/L$ and $E_{\rm
grav}/E_{a}$) are close to their counterparts in the thermalized state of the
GPPE (e.g., by comparing rows (2) and (3) in column(4) of Fig. 1 we see that
the final state of the SGLPE has nearly the same energy and radius as the
corresponding statistically steady GPPE state). Furthermore, convergence to
this thermalized state is more rapid in the (canonical) SGLPE than in the
(microcanonical) GPPE. To validate our DNSs, we have checked explicitly
(Supplemental Material supplement ) that, at zero temperature, our results
agree with those of the $T=0$ study of Refs. Chavanis1 ; Chavanis2 , which
yield spherically symmetric ground states ( row $(1)$ of Fig. 1) with radius
$R$ and $N$ bosons. Their ground-state energy can be approximated as
$E(R)=\frac{\hbar^{2}N}{2mR^{2}}+\frac{gN^{2}}{2R^{3}}-\frac{GN^{2}}{2R}$ and
the equilibrium radius follows from $dE/dR|_{R=R_{0}}=0$, whence
$R_{0}=\frac{R_{\rm Q}}{4\pi}(1+\sqrt{1+48\pi^{2}(\frac{R_{\rm a}}{R_{\rm
Q}}})^{2}$, where $R_{\rm Q}=\frac{\hbar^{2}}{G_{\rm N}m^{3}N}$ and $R_{\rm
a}=\sqrt{\frac{a\hbar^{2}}{G_{\rm N}m^{3}}}$. The details of our extensive
DNSs for the GPPE and the SGLPE are given in the Supplemental Material
supplement  (see, especially, Table I there). Figure 2(a) shows plots of of
$R/R_{\rm Q}$ versus $a/a_{\rm Q}$ for three different values of $G$. In Fig.
2(b) we give plots of $M/M_{\rm a}$ versus $R/R_{\rm a}$, from our GPPE DNSs
for three different values of $g$; the trends in these plots are markedly
different than those at $T=0$ (see Fig.(1) in the Supplemental Material
supplement ). In Fig. 2(c) we plot the scaled energies
$E_{kq}/E_{a},\,E_{int}/E_{a},\,E_{\rm G}/E_{\rm a},$ and their total
$E/E_{\rm a}$ versus the scaled temperature $k_{\rm B}T/E_{a}$ from our SGLPE
DNSs; the jumps in these curves at $k_{\rm B}T/E_{\rm a}\simeq 3.25\times
10^{-5}$ suggest a first-order transition from a collapsed BEC to a tenuous,
non-condensed assembly.

To confirm such a transition, we carry out an SGLPE hysteresis study, whose
results are summarised in the plots in the panels of Fig. 3. The left panel of
this figure shows $R/L$ versus $k_{\rm B}T/E_{\rm a}$; we obtain the red and
green curves by, respectively, increasing and decreasing $k_{\rm B}T/E_{\rm
a}$ (often referred to as heating and cooling runs in statistical mechanics).
In these SGLPE runs, we use the final steady-state configuration for
$\psi({\bf x})$, from the previous temperature, as the initial condition at
the next temperature; clearly, there is significant hysteresis at the first-
order transition from the - collapsed BEC to the non-collapsed state. In the
right side panels of Fig. 3, we show ten-level contour plots of
${\left|\psi({\bf x})\right|}^{2}$ and the associated spectra
${\left|\psi(k)\right|}^{2}$ to illustrate, at representative points on
heating and cooling curves in the hysteresis plot, the real-space density
distribution and the $k-$space density spectra ($k_{\rm B}T/E_{\rm
a}=2.7\times 10^{-5}$ and $k_{\rm B}T/E_{\rm a}=3.62\times 10^{-5}$ in the top
panels A and B, respectively, and $k_{\rm B}T/E_{\rm a}=2.3\times{10^{-5}}$
and $k_{\rm B}T/E_{\rm a}=3.16\times 10^{-5}$ in the bottom panels C and D,
respectively). From these density distributions and spectra, we conclude that
our system undergoes a first-order transition from a collapsed BEC to a
tenuous, non-condensed assembly. The panel of figures at the very bottom of
Fig. 3 show that this first-order transition occurs at $g=0$ too. We expect
that this also occurs when $g<0$, which is the appropriate parameter range for
axion stars Chavanis4 ; Chavanis5 ; the $g<0$ case requires a quintic
nonlinearity in Eq. 1 for stability; finite-temperature effects can be studied
for this case by using the methods that we have described above (as we will
show in future work).

Figure 2: The plots of (a) the scaled radius of gyration $(R/R_{\rm Q})$
versus scaled scattering length $(a/a_{\rm Q})$ for the GPPE Runs A1-A30, (b)
the scaled mass $(M/M_{\rm a})$ versus the scaled radius of gyration $R/R_{\rm
a}$ for GPPE Runs B1-B30 on a log-log scale, (c) the scaled energy components
($E_{\rm kq}/E_{\rm a},E_{\rm int}/E_{\rm a},E_{\rm G}/E_{\rm a})$, and the
total energy $(E/E_{\rm a})$ versus the scaled temperature $k_{\rm B}T/E_{\rm
a}$.

Figure 3: Top left panel: plots of the dimensionless radius ($R/L$) versus the
dimensionless temperature ($k_{\rm B}T/E_{a}$), for heating (red) and cooling
(green) runs showing a hysteresis loop. We show ten-level contour plots of
${\left|\psi({\bf x})\right|}^{2}$ and the associated spectra
${\left|\psi(k)\right|}^{2}$ to illustrate, at representative points on
heating and cooling curves in the hysteresis plot, the real-space density
distribution and the $k-$space density spectra ($k_{\rm B}T/E_{\rm
a}=2.7\times 10^{-5}$ and $k_{\rm B}T/E_{\rm a}=3.62\times 10^{-5}$ in panels
A and B of top panels, respectively, and $k_{\rm B}T/E_{\rm
a}=2.3\times{10^{-5}}$ and $k_{\rm B}T/E_{\rm a}=3.16\times 10^{-5}$ in panels
C and D of bottom panels, respectively). The analogs of these plots, for the
case $g=0$, are shown in the panels at the very bottom. In the bottom left
panel we use $|E_{G}|$ at $T=0$ to make the temperature dimensionless.

Do our Fourier-truncated GPPE yield only single collapsed objects? No. We now
show, for the illustrative parameter values $g=20$, $G=1000$, and $128^{3}$
collocation points, that this GPPE can also yield long-lived states with
temporal oscillations. In particular, by using an initial condition with two
rotating spherical compact objects, the truncated GPPE dynamics yields a
rotating binary system, which we depict at representative times in Fig. 4 and
in the Video V6 in the Supplementary Material supplement  (this also gives the
initial data).

Our study goes well beyond earlier studies harkofinite ; schivecosmic  that
employ the Thomas-Fermi approximation and include finite-temperature effects
at the level of a non-interacting Bose gas (with $g=0$). We have shown that
the truncated GPPE can be used effectively to study the graviational collapse
of an assembly of weakly interacting bosons at finite temperature. Our study
of this collapse shows that there is a clear, thermally driven first-order
transition from a non-condensed bosonic gas, at high $T$, to a condensed BEC
phase at low $T$; this transition should be contrasted with the continuous BEC
transition in the weakly interacting Bose gas, which is described by the GPE,
in the absence of gravitation. Furthermore, we have shown that gravitationally
bound, rotating binary objects can be obtained in our GPPE simulations.
Therefore, our work opens up the possibility of carrying out detailed finite-
temperature studies of self-gravitating bosonic systems, which are potentially
relevant for studies of dark-matter candidates like boson stars and axions.

Figure 4: Ten-level contour plots of the $|\psi({\bf x},t)|^{2}$ (a) at
$t=0.018$, (b) at $t=0.025$, and (c) at $t=0.03$, for the initial condition
for $\psi({\bf x},t)$ given in Eq.(1) of the Supplementary Material supplement
, showing rotating binary system (see video V6 in supplemental material).

Acknowledgements: We thank DST, CSIR (India), and the Indo French Center for
Applied Mathematics (IFCAM) for their support.

## References

  * (1) D. O’Dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin, Phys. Rev. Lett. 84, 5687 (2000).
  * (2) R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * (3) G. Ingrosso, D. Grasso, and R. Ruffini, Astron. Astro- phys. 248, 481 (1991).
  * (4) P. Jetzer, Physics Reports 220, 163 (1992).
  * (5) W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000).
  * (6) H. Abdallah et al., Phys. Rev. Lett. 120, 201101 (2018).
  * (7) C. Kachulis et al., Phys. Rev. Lett. 120, 221301v(2018).
  * (8) https://physics.aps.org/articles/v11/48?utm_campaign=weekly&utm_medium=email&utm_source=emailalert
  * (9) M. Lawson et al., Phys. Rev. Lett. 123, 141802 (2019).
  * (10) P.-H. Chavanis, Phys. Rev. D 84, 043531 (2011).
  * (11) P.-H. Chavanis and L. Delfini, Phys. Rev. D 84, 043532 (2011).
  * (12) P.-H. Chavanis, Phy. Rev. D 94, 083007 (2016).
  * (13) P.-H. Chavanis, Phys. Rev. D 98, 023009 (2018).
  * (14) P.-H. Chavanis, Quantum Aspects of Black Holes, eds. X. Calmet et al, (Springer, 2015), pp. $151–194$.
  * (15) N. P. Proukakis and B. Jackson, Journal of Physics B: Atomic, Molecular and Optical Physics 41, 203002 (2008).
  * (16) N. G. Berloff, M. Brachet, and N. P. Proukakis, Proc. Natl. Acad. Sci. U.S.A. 111, 4675 (2014).
  * (17) G. Krstulovic and M. Brachet, Phys. Rev. E 83, 066311 (2011).
  * (18) V. Shukla, M. Brachet, and R. Pandit, New J. Phys. 15, 113025 (2013).
  * (19) Supplemental Material (2019).
  * (20) P. Peebles, The large-scale structure of the universe, 2011 (Princeton University Press, Princeton, NJ, 1980).
  * (21) M. Falco, S. H. Hansen, R. Wojtak, and G. A. Mamon, Mon. Not. R. Astron. Soc. 431, L6 (2013).
  * (22) M. K.-H. Kiessling, Advances in Applied Mathematics 31, 132 (2003).
  * (23) D. Gottlieb, and S. A. Orszag, Numerical Analysis of Spectral Methods (SIAM, Philadelphia, 1977).
  * (24) T. Harko, and J.M. Madarassy, J. Cosmol. Astropart. Phys. 01 (2012) 025.
  * (25) H-Y. Schive, T. Chiueh, and T. Broadhurst, Nat. Phys. 10, 496 (2014) .