# Topology of Hawking-Page transition in Born-Infeld AdS black holes

Pavan Kumar Yerra1,2 Chandrasekhar Bhamidipati3 and Sudipta Mukherji1,2 1
Institute of Physics, Sachivalaya Marg, Bhubaneswar, Odisha, 751005, India  2
Homi Bhabha National Institute, Training School Complex, Anushakti Nagar,
Mumbai, 400085, India  3 School of Basic Sciences, Indian Institute of
Technology Bhubaneswar, Bhubaneswar, Odisha, 752050, India
pk11@iitbbs.ac.in,chandrasekhar<EMAIL_ADDRESS>
###### Abstract

Black holes in anti de Sitter (AdS) spacetimes undergo phase transitions which
typically lead to the existence of critical points, that can be classified
using topological techniques. Availing the Bragg-Williams construction which
provides an off-shell free energy formalism, we compute the topological charge
of the Hawking-Page (HP) transition for Einstein-Born-Infeld black holes in
anti de Sitter spacetime and match the result with the confinement-
deconfinement transition in the dual gauge theory, which turn out to be in
perfect agreement.

##  1 Introduction

Thermodynamics of Black holes and their phase transitions, particularly in
anti de Sitter space times (AdS) have continued to provide interesting
insights [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13], with a special place for
the Hawking-Page transition [14], which can be interpreted as a confinement-
deconfinement transition in the dual gauge/gravity setting [15, 16, 17]. Mean
field approximation turns out to be a reliable approach to study phase
transitions, not just in black holes, but for general thermodynamic systems
[18, 19, 20, 21, 22]. A particularly elegant approach is the one by Bragg-
Williams, where one starts by constructing an off shell free energy given in
terms of an order parameter. Minimising this free energy then gives the
required phases of the thermodynamic system. This approach has been applied
successfully to black holes in AdS as well [22, 23].

Novel topological techniques have be pioneered to classify critical points of
black holes111see [24, 25, 26, 27, 28, 29] for other developments which were
initiated earlier from the study of light-rings [30, 31, 32]. The method and
the results have now been generalised to several interesting situations [33,
34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51]. It is
now understood that topological charges can be assigned not just to first and
second order phase transition points, but to the black hole solutions
themselves. In particular, it has now been possible to compute the topological
charge of the confinement-deconfinement transition using an effective
potential in the dual gauge theory and match the result with the charge of
Hawking-Page transition point in the bulk [52]. For the case of charged black
hole in AdS, it is found that the two topological charges are same, even
though the order parameters in the bulk and boundary are quite different,
namely, horizon radius $r_{+}$ for the Hawking-Page transition point, and a
charge parameter $Q$ for the confinement-deconfinement transition in gauge
theory. There are further checks possible, knowing that the saddle points of
the off shell Bragg-Williams free energy give the equilibrium phases of the
thermodynamic system [18, 53, 22, 54]. Using the inputs from [32], one can
assign topological charges to stable and unstable phases of the system, which
turn out to have opposite values [36], both in the bulk and in the boundary
[52, 55]. The aim of this short article is to show that the computation of
topological charges of Hawking-Page transition point and equilibrium phases of
black holes done in [52], can be extended to the case of Born-Infeld AdS black
holes [56] in the grand canonical ensemble.

Rest of the article is structured as follows. In section-(2), we present the
Bragg-Williams off-shell free energy for Born-Infeld AdS black holes and
locate the Hawking-Page transition point. Subsection-(2.1) contains the
computation of topological charge of the HP point. Section-(3) contains the
effective potential of the boundary gauge theory and the computation of
topological charge of confinement-deconfinement transition point. We conclude
in section-(4) with some open questions.

##  2 Bulk

We start with the Bragg-Williams off-shell free energy $f$ of Born-Infeld AdS
black holes in grand canonical ensemble (i.e., fixed potential $\mu$) in
$(n+1)$-dimensional spacetime, given by [56, 53]:

$f(r_{+},T,\mu)=M-TS-\mu Q,$ (2.1)

where, the mass $M$, entropy $S$, potential $\mu$ and charge $Q$ of the black
hole are given in terms of the horizon radius $r_{+}$ as:

$\displaystyle M$ $\displaystyle=$
$\displaystyle\omega_{n-1}(n-1)\bigg{\\{}r_{+}^{n-2}+\frac{r_{+}^{n}}{l^{2}}+\frac{4\beta^{2}r_{+}^{n}}{n(n-1)}\left(1-\sqrt{1+\frac{(n-2)(n-1)q^{2}}{2\beta^{2}r_{+}^{2n-2}}}\right)$
$\displaystyle\,\,\,+\frac{2(n-1)q^{2}}{nr_{+}^{n-2}}\,_{2}F_{1}\big{[}\frac{n-2}{2(n-1)},\frac{1}{2},\frac{3n-4}{2(n-1)},-\frac{(n-1)(n-2)q^{2}}{2\beta^{2}r_{+}^{2n-2}}\big{]}\bigg{\\}},$
$\displaystyle S$ $\displaystyle=$ $\displaystyle
4\pi\omega_{n-1}r_{+}^{n-1},$ $\displaystyle\mu$ $\displaystyle=$
$\displaystyle\frac{q}{cr_{+}^{n-2}}\,_{2}F_{1}\big{[}\frac{n-2}{2(n-1)},\frac{1}{2},\frac{3n-4}{2(n-1)},-\frac{(n-1)(n-2)q^{2}}{2\beta^{2}r_{+}^{2n-2}}\big{]},$
$\displaystyle Q$ $\displaystyle=$ $\displaystyle
2\sqrt{2(n-1)(n-2)}\,\omega_{n-1}\,q,$ (2.2)

with $\omega_{n-1}$ being the volume of a unit $(n-1)$ sphere,
$c=\sqrt{\frac{2(n-2)}{n-1}}$, $l$ is the AdS length, $\beta$ is the Born-
Infeld coupling, ${}_{2}F_{1}$ is the hypergeometric function, and $q$ is an
integration constant related to the charge $Q$.

(a)

(b)

Figure 1: (a) Bragg-Williams free energy $f$ is plotted against $r_{+}$ (using
$x$ as a parameter) at different temperatures $T$ for a fixed $\mu$.
Temperature of the curves increases from top to bottom. Hawking-Page
transition happens for the temperature $T_{\rm HP}$ (dashed red curve) at
$r_{\rm HP}$ (red dot). (b) The curve $T_{\rm 0}$, as a function of $r_{+}$
(using $x$ as a parameter), shows the HP transition point at its minima (red
dot). Here, $(T_{\rm HP},r_{\rm HP},x_{\rm HP})=(0.464,0.972,0.243)$. Plots
are displayed for $n=4,\mu=0.2,\omega_{3}=l=\beta=1$.

Treating horizon radius $r_{+}$ as order parameter, $T$ and $\mu$ as the
external parameters, the behaviour of the free energy $f$ (plotted in the Fig.
1a), shows the Hawking-Page (HP) transition, which happens when:

$f=0,\,\text{and},\,\frac{\partial f}{\partial r_{+}}=0.$ (2.3)

Alternatively (which will be useful for later purpose), one can also find the
Hawking-Page transition point using the first condition of eqn. 2.3, which
gives,

$f=0\implies T=T_{\rm 0}^{\phantom{0}}(r_{+},\mu).$ (2.4)

The Hawking-Page transition point can be located at the minima of the curve
$T_{\rm 0}^{\phantom{0}}$, as shown in the Fig. 1b. We note here that, due to
the non-linear relation between $\mu$ and $q$ as in eqn. 2, we study the free
energy $f$ by defining a parameter $x$ as [53]:

$x=\frac{q}{r_{+}^{n-1}}.$ (2.5)

The horizon radius, $r_{+}$ can now be rewritten as

$r_{+}=\frac{\mu
c}{x\,_{2}F_{1}\big{[}\frac{n-2}{2(n-1)},\frac{1}{2},\frac{3n-4}{2(n-1)},-\frac{(n-1)(n-2)x^{2}}{2\beta^{2}}\big{]}}.$
(2.6)

The HP point can be defined now using:

$f=0,\,\text{and},\,\frac{\partial f}{\partial x}=0,$ (2.7)

or, from the minima of the curve $T_{0}(x,\mu)$.

###  2.1  Topological charge of Hawking-Page transition point

In order to assign the topological charge to HP transition point, we employ
the temperature $T_{\rm 0}^{\phantom{0}}(x,\mu)$ to define the vector filed
$\phi(\phi^{x},\phi^{\theta})$, in the following way:

$\displaystyle\phi^{x}$ $\displaystyle=$
$\displaystyle(\partial_{x}\Phi)_{\theta,\mu}\ $ (2.8)
$\displaystyle\phi^{\theta}$ $\displaystyle=$
$\displaystyle(\partial_{\theta}\Phi)_{x,\mu}\ $ (2.9)

where, $\Phi=\frac{1}{\text{sin}\theta}T_{\rm 0}^{\phantom{0}}(x,\mu)$. This
vector field $\phi$ vanishes exactly at the Hawking-Page transition point,
which can be seen clearly from the normalised vector field
$(\frac{\phi^{x}}{||\phi||},\frac{\phi^{\theta}}{||\phi||})$ plot in the Fig.
2a.

(a)

(b)

Figure 2: (a) The normalised vector field in the $\theta-x$ plane, shows the
presence of Hawking-Page (HP) transition point (black dot, located at $x_{\rm
HP}=0.243$) at $\phi=0$. Contour $C$ encloses the HP transition point. (b)
$\Omega$ vs $\vartheta$ for contour $C$. Plots are displayed for
$n=4,\mu=0.2,\omega_{3}=l=\beta=1$.

The behaviour of the deflection angle $\Omega(\vartheta)$, for the given
contour $C$ in the Fig. 2a, is as shown in Fig 2b, from where we find that the
topological charge associated with the HP transition point would be
$Q_{t}\big{|}_{\rm HP}^{\phantom{HP}}=\frac{1}{2\pi}\Omega(2\pi)=+1$.

##  3 Boundary gauge theory

The effective potential in the gauge theory dual to Born-Infeld AdS black
holes (in grand canonical ensemble) in $n=4$ is [53]:

$\displaystyle W$ $\displaystyle=$
$\displaystyle\omega_{3}\,N_{c}^{2}\bigg{\\{}-\frac{Tr_{+}^{3}}{2\pi}-Q\mu+\frac{3}{8\pi^{2}}\bigg{(}\frac{r_{+}^{4}}{l^{2}}+r_{+}^{2}+\frac{\beta^{2}r_{+}^{4}}{3}\Big{(}1-\sqrt{1+\frac{4\pi^{4}Q^{2}}{\beta^{2}r_{+}^{6}}}\Big{)}$
(3.1)
$\displaystyle+\frac{2\pi^{4}Q^{2}\,_{2}F_{1}\big{[}\frac{1}{3},\frac{1}{2},\frac{4}{3},-\frac{4\pi^{4}Q^{2}}{\beta^{2}r_{+}^{6}}\big{]}}{r_{+}^{2}}\bigg{)}\bigg{\\}}.$

where, $N_{c}$ stands for number of colours. In order to study the phase
structure, again we use the parameter $x$ defined in eqn. (2.5), and carry out
a parametric plot of $W$ against $Q$, the order parameter in the boundary
theory. The resulting phase structure [Fig. 3a] shows a confinement-
deconfinement transition at critical temperature $T_{c}$, which is exactly at
HP temperature $T_{\rm HP}$.

(a)

(b)

Figure 3: For boundary gauge theory: (a) The effective potential $W$ is
plotted against $Q$ (using $x$ as a parameter) at different temperatures $T$
for a fixed $\mu$. Temperature of the curves increases from top to bottom.
Confinement-deconfinement transition happens for the temperature $T_{\rm c}$
(dashed red curve) at $Q_{\rm c}$ (red dot). (b) The curve $T_{\rm 0}$, as a
function of $Q$ (using $x$ as a parameter), shows the confinement-
deconfinement transition point at its minima (red dot). Here, $(T_{\rm
c},Q_{\rm c},x_{\rm c})=(0.464,0.0196,0.243)$. Plots are displayed for
$n=4,\mu=0.2,N_{c}=\omega_{3}=l=\beta=1$.

The critical temperature for the confinement-deconfinement transition can be
obtained on satisfying the following two conditions simultaneously:

$W=0,\,\text{and},\,\frac{\partial W}{\partial Q}=0.$ (3.2)

Alternatively, as shown in the Fig. 3b, one can also find $T_{c}$ from the
minima of the curve $T_{0}$, using $W=0$. We note here that, using the
parameter $x$, the confinement-deconfinement point can again be located at the
minima of the curve $T_{0}(x,\mu)$.

(a)

(b)

Figure 4: For boundary gauge theory: (a) The vector field $\phi$ vanishes at
the confinement-deconfinement transition point (black dot, located at $x_{\rm
c}$). Contour $C$ encloses the transition point. (b) $\Omega$ vs $\vartheta$
for contour $C$. Plots are displayed for
$n=4,\mu=0.2,N_{c}=\omega_{3}=l=\beta=1$.

Now, using $T_{0}(x,\mu)$, we define the vector filed
$\phi(\phi^{x},\phi^{\theta})$, as:

$\displaystyle\phi^{x}$ $\displaystyle=$
$\displaystyle(\partial_{x}\Phi)_{\theta,\mu}\ $ (3.3)
$\displaystyle\phi^{\theta}$ $\displaystyle=$
$\displaystyle(\partial_{\theta}\Phi)_{x,\mu}\ $ (3.4)

where, $\Phi=\frac{1}{\text{sin}\theta}T_{0}(x,\mu)$. This vector field $\phi$
vanishes exactly at the confinement-deconfinement transition point, as can be
seen from the Fig. 4a. Then, the topological charge corresponding to this
transition point turns out to be $Q_{t}\big{|}_{\rm
c}^{\phantom{c}}=\frac{1}{2\pi}\Omega(2\pi)=+1$ (given by the contour $C$ in
Fig. 4b). This charge exactly matches the one obtained from bulk calculation.

##  4 Conclusions

Following the topological classification of phase transition points of black
holes discussed in [32, 33, 34, 35, 36, 52], we computed the topological
charge of the Hawking-Page transition point in Einstein-Born-Infeld black
holes in AdS using an off shell Bragg-Williams free energy approach and found
it to be $+1$. It was shown that the same value appears from an analogous
computation performed using an effective potential in the boundary gauge
theory. Topological ideas from other directions are being advanced to
understand critical points with interesting connections to quantum error
correction [57, 58]. It seems plausible to attempt a topological
classification of critical points in matrix models and possibly general
thermodynamic systems [55]. It is expected that classification will reveal
important information about the nature of phase transitions in a wide range of
theories.

C.B. thanks IIT Bhubaneswar for financial support, and the DST (SERB),
Government of India, through the Mathematical Research Impact Centric Support
(MATRICS) grant no. MTR/2020/000135. C.B. also thanks Renann Lipinski
Jusinskas, the Institute of Physics (FZU) of the Czech Academy of Sciences &
CEICO, Prague, for warm hospitality, and the organisers of the XII.
International Symposium on Quantum Theory and Symmetries (QTS12), Czech
Technical University, Prague for a nice conference.

## References

  * [1] Bekenstein J D 1973 Phys. Rev. D 7 2333–2346
  * [2] Bardeen J M, Carter B and Hawking S W 1973 Commun. Math. Phys. 31 161–170
  * [3] Hawking S W 1975 Commun. Math. Phys. 43 199–220 [Erratum: Commun.Math.Phys. 46, 206 (1976)]
  * [4] Chamblin A, Emparan R, Johnson C V and Myers R C 1999 Phys. Rev. D 60 064018 (Preprint hep-th/9902170)
  * [5] Caldarelli M M, Cognola G and Klemm D 2000 Class. Quant. Grav. 17 399–420 (Preprint hep-th/9908022)
  * [6] Kastor D, Ray S and Traschen J 2009 Class. Quant. Grav. 26 195011 (Preprint 0904.2765)
  * [7] Cvetic M, Gibbons G W, Kubiznak D and Pope C N 2011 Phys. Rev. D84 024037 (Preprint 1012.2888)
  * [8] Dolan B P 2011 Class. Quant. Grav. 28 235017 (Preprint 1106.6260)
  * [9] Karch A and Robinson B 2015 JHEP 12 073 (Preprint 1510.02472)
  * [10] Kubiznak D and Mann R B 2012 JHEP 07 033 (Preprint 1205.0559)
  * [11] Kubiznak D, Mann R B and Teo M 2017 Class. Quant. Grav. 34 063001 (Preprint 1608.06147)
  * [12] Banerjee R and Roychowdhury D 2012 Phys. Rev. D 85 044040 (Preprint 1111.0147)
  * [13] Lu J X, Roy S and Xiao Z 2011 JHEP 01 133 (Preprint 1010.2068)
  * [14] Hawking S W and Page D N 1983 Commun. Math. Phys. 87 577
  * [15] Maldacena J M 1998 Adv. Theor. Math. Phys. 2 231–252 (Preprint hep-th/9711200)
  * [16] Gubser S S, Klebanov I R and Polyakov A M 1998 Phys. Lett. B 428 105–114 (Preprint hep-th/9802109)
  * [17] Witten E 1998 Adv. Theor. Math. Phys. 2 253–291 (Preprint hep-th/9802150)
  * [18] Chaikin P M and Lubensky T C 1995 Principles of Condensed Matter Physics (Cambridge: Cambridge University Press)
  * [19] Bragg W L and Williams E J 1934 Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 145 699–730
  * [20] Bragg W L and Williams E J 1935 Proceedings of the Royal Society of London. Series A - Mathematical and Physical Sciences 151 540–566
  * [21] Kubo R Statistical Mechanics: An Advanced Course with Problems and Solutions Ryōgo Kubo ISBN 9780720400991
  * [22] Banerjee S, Chakrabarti S K, Mukherji S and Panda B 2011 Int. J. Mod. Phys. A 26 3469–3489 (Preprint 1012.3256)
  * [23] Chandrasekhar B, Mukherji S, Sahay A and Sarkar S 2012 JHEP 05 004 (Preprint 1202.4059)
  * [24] Cunha P V P, Berti E and Herdeiro C A R 2017 Phys. Rev. Lett. 119 251102 (Preprint 1708.04211)
  * [25] Cunha P V P and Herdeiro C A R 2020 Phys. Rev. Lett. 124 181101 (Preprint 2003.06445)
  * [26] Guo M and Gao S 2021 Phys. Rev. D 103 104031 (Preprint 2011.02211)
  * [27] Wei S W 2020 Phys. Rev. D 102 064039 (Preprint 2006.02112)
  * [28] Junior H C D L, Yang J Z, Crispino L C B, Cunha P V P and Herdeiro C A R 2022 Phys. Rev. D 105 064070 (Preprint 2112.10802)
  * [29] Wu S P and Wei S W 2023 Phys. Rev. D 108 104041 (Preprint 2307.14003)
  * [30] Duan Y S 1984 SLAC-PUB-3301
  * [31] Duan Y S and Ge M L 1979 Sci. Sin. 9 1072
  * [32] Wei S W and Liu Y X 2022 Phys. Rev. D 105 104003 (Preprint 2112.01706)
  * [33] Yerra P K and Bhamidipati C 2022 Phys. Rev. D 105 104053 (Preprint 2202.10288)
  * [34] Ahmed M B, Kubiznak D and Mann R B 2023 Phys. Rev. D 107 046013 (Preprint 2207.02147)
  * [35] Yerra P K and Bhamidipati C 2022 Phys. Lett. B 835 137591 (Preprint 2207.10612)
  * [36] Wei S W, Liu Y X and Mann R B 2022 Phys. Rev. Lett. 129 191101 (Preprint 2208.01932)
  * [37] Barzi F, El Moumni H and Masmar K 2023 (Preprint 2309.14069)
  * [38] Gogoi N J and Phukon P 2023 Phys. Rev. D 107 106009
  * [39] Alipour M R, Afshar M A S, Noori Gashti S and Sadeghi J 2023 (Preprint 2305.05595)
  * [40] Gogoi N J and Phukon P 2023 Phys. Rev. D 108 066016 (Preprint 2304.05695)
  * [41] Wu D 2023 Phys. Rev. D 107 024024 (Preprint 2211.15151)
  * [42] Fang C, Jiang J and Zhang M 2023 JHEP 01 102 (Preprint 2211.15534)
  * [43] Fan Z Y 2023 Phys. Rev. D 107 044026 (Preprint %****␣conference.bbl␣Line␣175␣****2211.12957)
  * [44] Liu C and Wang J 2023 Phys. Rev. D 107 064023 (Preprint 2211.05524)
  * [45] Bai N C, Li L and Tao J 2023 Phys. Rev. D 107 064015 (Preprint 2208.10177)
  * [46] Zhang M Y, Chen H, Hassanabadi H, Long Z W and Yang H 2023 Eur. Phys. J. C 83 773 (Preprint 2305.15674)
  * [47] Hung T N and Nam C H 2023 Eur. Phys. J. C 83 582 (Preprint 2305.15910)
  * [48] Wang Y S, Xu Z M and Wu B 2023 (Preprint 2307.01569)
  * [49] Chen D, He Y, Tao J and Yang W 2023 (Preprint 2311.11606)
  * [50] Sadeghi J, Afshar M A S, Noori Gashti S and Alipour M R 2024 Annals Phys. 460 169569 (Preprint 2312.04325)
  * [51] Shahzad M U, Mehmood A, Sharif S and Övgün A 2023 Annals Phys. 458 169486
  * [52] Yerra P K, Bhamidipati C and Mukherji S 2022 Phys. Rev. D 106 064059 (Preprint 2208.06388)
  * [53] Banerjee S 2010 Phys. Rev. D 82 106008 (Preprint 1009.1780)
  * [54] Nayak B P 2008 Prayas Students’ Journal of Physics
  * [55] Yerra P K, Bhamidipati C and Mukherji S 2023 (Preprint 2304.14988)
  * [56] Dey T K 2004 Phys. Lett. B 595 484–490 (Preprint hep-th/0406169)
  * [57] Almheiri A, Dong X and Harlow D 2015 JHEP 04 163 (Preprint 1411.7041)
  * [58] Bao N, Cao C and Zhu G 2022 Phys. Rev. D 106 046009 (Preprint 2202.04710)


# Gravitational waves from a scotogenic two-loop neutrino mass model

Cesar Bonilla  A. E. Cárcamo Hernández  João Gonçalves  Vishnudath K. N
António P. Morais  and Roman Pasechnik

###### Abstract

We propose a framework to account for neutrino masses at the two-loop level.
This mechanism introduces new scalars and Majorana fermions to the Standard
Model. It is assumed the existence of a global $U(1)\times\mathcal{Z}_{2}$
symmetry which after partial breaking provides the stability of the dark
matter candidates of the theory. The rich structure of the potential allows
for the possibility of first-order phase transitions (FOPTs) in the early
Universe which can lead to the generation of primordial gravitational waves as
one of the potentially observable signatures of this model. Taking into
account relevant constraints from lepton flavour violation, neutrino physics
as well as the trilinear Higgs couplings at next-to-leading order accuracy, we
have found a wide range of possible FOPTs which are strong enough to be probed
at the proposed gravitational-wave interferometer experiments such as LISA.

CERN-TH-2023-071

##  1 Introduction

The Standard Model (SM) is a highly successful theory that describes the
strong, electromagnetic, and weak interactions whose predictions have been
experimentally verified at the LHC with the highest degree of accuracy.
However it has several unaddressed issues such as, the observed pattern of SM
fermion masses and mixing angles, the number of SM fermion families, the
measured amount of Dark Matter (DM) relic density and baryon asymmetry
observed in the Universe, among others. Experiments with solar, atmospheric,
and reactor neutrinos have brought evidence of neutrino oscillations as a
signature of their nonzero masses and mixing. Several extensions of the SM
have been constructed in order to explain the tiny masses of the active
neutrinos; see for instance Ref. [1] for a review.

The most economical way to generate the small neutrino masses considering the
SM gauge symmetry is by adding two right-handed Majorana neutrinos that mix
with the light active neutrinos thus triggering a type-I seesaw mechanism [2,
3, 4, 5, 6, 7, 8]. In this mechanism, either the right-handed Majorana
neutrinos have to be extremely heavy, with masses of order the Grand
Unification scale, or they can be around the TeV scale implying very tiny
Dirac Yukawa couplings, in order to successfully reproduce the neutrino data.
In both scenarios, the mixing between the active and sterile neutrinos is very
small, making the detection at the collider experiments difficult and strongly
suppressed charged lepton flavor violating (LFV) signatures, several orders of
magnitude below the experimental sensitivity. This makes models with tree-
level type-I seesaw realizations difficult to test at the particle physics
experiments.

Alternatively, radiative seesaw models are examples of interesting and
testable extensions of the SM naturally explaining tiny neutrino masses and
the absence of a strong hierarchy between the neutrino mixing angles. In most
radiative seesaw models, neutrino masses arise at a one-loop level. This
implies that in order to successfully reproduce the experimental neutrino
data, one has to rely either on very small neutrino Yukawa couplings (of order
the electron Yukawa coupling) or on an unnaturally small mass splitting
between the CP-even and CP-odd components of the neutral scalar mediators.
Two-loop neutrino mass models have been proposed in the literature [9, 10, 11,
12, 13, 14, 15, 16, 17] to naturally provide further suppressed active
neutrino masses than those based on one-loop level radiative seesaw
mechanisms. See [18, 19] for comprehensive studies of one- and two-loop
neutrino mass models.

Our goal here is to explore a particularly attractive new scenario for
neutrino mass generation at the two-loop level based upon an additional global
$U(1)\times\mathcal{Z}_{2}$ symmetry in the lepton sector, featuring also
inert scalars and heavy Majorana neutrinos. Among the most promising
phenomenological signatures of this model are the emergence of strong first-
order phase transitions (FOPTs) that can be studied via a possible detection
of a stochastic gravitational wave background (SGWB) as well as a sizeable
modification in the trilinear Higgs coupling that can be probed at the LHC and
future colliders. In this work, we study intricate correlations among these
signatures and the properties of the neutrino sector.

The paper is organised as follows. In Sec. 2, we introduce our model with a
focus on the Yukawa and scalars sectors; in Sec. 3 we discuss the
phenomenological implications of the new additional states in the context of
neutrino physics and charged LFV observables; in Sec. 4 we introduce the main
building blocks for the finite-temperature description of our model and how it
can be correlated with SGWB observables, also showcasing our numerical
results; finally, we finish with Sec. 5 where we take our conclusions.

##  2 Theoretical structure of the model

###  2.1 Particle content and charge assignments

In this work, we augment the SM gauge groups with an extended global abelian
symmetry $\mathrm{U(1)}\times\mathcal{Z}_{2}$, where the $\mathrm{U(1)}$ is
spontaneously broken while the $\mathcal{Z}_{2}$ is preserved. Besides the
standard particle content of the SM, additional fields are present in the
spectrum, including an inert SU(2) doublet scalar field, $\eta$, and two
scalar singlets, $\varphi$ and $\sigma$. In our setup, $\sigma$ acquires a
nonzero vacuum expectation value (VEV) while $\varphi$ does not, such that the
potential DM candidates can arise from the neutral components of $\eta$ or
from $\varphi$. Additional Majorana fermions are also included into this
scheme, namely, two $N_{R}$ fields and two $\Psi_{R}$. While the electroweak
(EW) gauge group does not distinguish $N_{R}$ and $\Psi_{R}$, the extended
symmetry group treats both differently.

The charge assignments of the fields are shown in Tab. 1 for the scalar fields
and in Tab. 2 for matter fields. As shown in Tab. 1, the SU(2) scalar doublet
$\eta$ and the scalar singlet $\varphi$ have non-trivial $\mathcal{Z}_{2}$
charges. Since these two scalars do no acquire VEVs and the charge assignments
do not let neutrinos to couple to the SM Higgs, light active neutrinos do not
acquire masses at tree level. The specific
$\mathrm{U(1)}\times\mathcal{Z}_{2}$ charge assignments of the
$N_{k,\mathrm{R}}$ Majorana neutrinos and $\sigma$ scalar field enables the
radiative generation of neutrino mass and the corresponding mixing only at
two-loop level that can be noticed from a typical topology in Fig. 1. Namely,
this is possible through the virtual exchange of the neutral components of the
inert SU(2) scalar doublet $\eta$ fields as well as of the real and imaginary
parts of the inert singlet scalar $\varphi$, with the loop closing through the
$\sigma$ VEV and the VEV of the active Higgs doublet $\Phi$. For simplicity of
the analysis, we assume that the charged lepton sector is purely diagonal,
such that the only contribution to the Pontecorvo–Maki–Nakagawa–Sakata (PMNS)
mixing matrix comes exclusively from the loop topology shown in Fig. 1.

| $\mathrm{SU(3)_{C}}$ | $\mathrm{SU(2)_{L}}$ | $\mathrm{U(1)_{Y}}$ | $\mathrm{U(1)_{X}}$ | $\mathcal{Z}_{2}$
---|---|---|---|---|---
$\Phi$ | $\bm{1}$ | $\bm{2}$ | $1/2$ | $0$ | $1$
$\sigma$ | $\bm{1}$ | $\bm{1}$ | $0$ | $-1/2$ | $1$
$\eta$ | $\bm{1}$ | $\bm{2}$ | $1/2$ | $0$ | $-1$
$\varphi$ | $\bm{1}$ | $\bm{1}$ | $0$ | $-1/2$ | $-1$
Table 1: Charge assignments under the SM gauge group ($\mathrm{SU(3)_{C}}\times\mathrm{SU(2)_{L}}\times\mathrm{U(1)_{Y}}$) and the supplemental global symmetry $\mathrm{U(1)_{X}}\times\mathcal{Z}_{2}$ for the new scalar fields ($\sigma$, $\eta$ and $\varphi$) and the Higgs doublet ($\Phi$). | $\mathrm{SU(3)_{C}}$ | $\mathrm{SU(2)_{L}}$ | $\mathrm{U(1)_{Y}}$ | $\mathrm{U(1)_{X}}$ | $\mathcal{Z}_{2}$
---|---|---|---|---|---
$\ell_{i,\mathrm{L}}$ | $\bm{1}$ | $\bm{2}$ | $-1/2$ | $1/2$ | $1$
$\ell_{i,\mathrm{R}}$ | $\bm{1}$ | $\bm{1}$ | $-1$ | $1/2$ | $1$
$N_{k,\mathrm{R}}$ | $\bm{1}$ | $\bm{1}$ | $0$ | $1/2$ | $-1$
$\Psi_{k,\mathrm{R}}$ | $\bm{1}$ | $\bm{1}$ | $0$ | $0$ | $1$
Table 2: Charge assignments under the SM gauge group
($\mathrm{SU(3)_{C}}\times\mathrm{SU(2)_{L}}\times\mathrm{U(1)_{Y}}$) and the
supplemental global symmetry $\mathrm{U(1)_{X}}\times\mathcal{Z}_{2}$ for the
new matter fields ($N_{R}$, $\Psi_{R}$) and the SM-like lepton doublet
($\ell_{i,\mathrm{L}}$) and singlet ($\ell_{i,\mathrm{R}}$). Here, $i=1,2,3$
and $k=1,2$.

###  2.2 Yukawa sector and scalar potential

Based on the charge assignments of Tabs. 1 and 2, the most general and
renormalisable Yukawa Lagrangian reads as

$\displaystyle-\mathcal{L}_{Y}^{(l)}=$
$\displaystyle(y_{\ell})_{ij}\leavevmode\nobreak\
\overline{\ell}_{i,\mathrm{L}}\Phi\ell_{j,\mathrm{R}}+(y_{N})_{ik}\leavevmode\nobreak\
\overline{\ell}_{i,\mathrm{L}}\tilde{\eta}N_{k,\mathrm{R}}+(y_{\Omega})_{nk}\leavevmode\nobreak\
\overline{N}_{n,\mathrm{R}}\varphi^{\ast}\Psi_{k,\mathrm{R}}^{c}+(m_{\Psi})_{nm}\leavevmode\nobreak\
\overline{\Psi}_{n,\mathrm{R}}\Psi_{m,\mathrm{R}}^{c}+\mathrm{h.c.}\,,$ (2.1)

where we have defined $\tilde{\eta}=i\sigma_{2}\eta^{\dagger}$ and the
superscript $c$ in the $\Psi$ field indicates charge conjugation. We adopt the
Einstein summation convention where repeated indices indicate sum over them,
with $i,j=1,2,3$ and $m,n,k=1,2$. Here, $y_{\ell}$ is a $3\times 3$ matrix,
$y_{\Omega}$ and $m_{\psi}$ are $2\times 2$ matrices and $y_{N}$ is a $2\times
3$ matrix. The neutrino mass matrix can then be induced by the combination of
the $y_{N}$ and $y_{\Omega}$ Yukawa matrices.

Figure 1: Two-loop Feynman diagram which is responsible for the generation of
the neutrino masses and mixing. Here, $\nu$ are the neutrinos, $\eta_{R,I}$
are the real and imaginary parts of the neutral component of the $\eta$
doublet, $\varphi_{R,I}$ are the real and imaginary parts of the $\varphi$
singlet, $v$ and $v_{\sigma}$ are the VEVs of the $\Phi$ doublet the $\sigma$
singlet fields, respectively.

The complete gauge invariant scalar potential in this model is given as,

$\displaystyle V=$
$\displaystyle-\mu_{\Phi}^{2}(\Phi^{\dagger}\Phi)+\mu_{\eta}^{2}(\eta^{\dagger}\eta)+\mu_{\varphi}^{2}(\varphi^{\ast}\varphi)-\mu_{\sigma}^{2}(\sigma^{\ast}\sigma)-\mu_{sb}^{2}\left(\sigma^{2}+h.c\right)+\lambda_{1}(\Phi^{\dagger}\Phi)(\Phi^{\dagger}\Phi)+\lambda_{2}(\eta^{\dagger}\eta)(\eta^{\dagger}\eta)$
$\displaystyle+\lambda_{3}(\varphi^{\ast}\varphi)(\varphi^{\ast}\varphi)+\lambda_{4}(\sigma^{\ast}\sigma)(\sigma^{\ast}\sigma)+\lambda_{5}(\Phi^{\dagger}\Phi)(\eta^{\dagger}\eta)+\lambda_{6}(\Phi^{\dagger}\eta)(\eta^{\dagger}\Phi)+\frac{\lambda_{7}}{2}\left[(\Phi^{\dagger}\eta)^{2}+\mathrm{h.c.}\right]$
$\displaystyle+\lambda_{8}(\Phi^{\dagger}\Phi)(\sigma^{\ast}\sigma)+\lambda_{9}(\Phi^{\dagger}\Phi)(\varphi^{\ast}\varphi)+\lambda_{10}(\eta^{\dagger}\eta)(\sigma^{\ast}\sigma)+\lambda_{11}(\eta^{\dagger}\eta)(\varphi^{\ast}\varphi)+\lambda_{12}(\varphi^{\ast}\varphi)(\sigma^{\ast}\sigma)$
$\displaystyle+\frac{\lambda_{13}}{2}\left[\varphi^{2}\left(\sigma^{\ast}\right)^{2}+\mathrm{h.c.}\right]\,,$

where we take all parameters of the potential to be real, i.e. CP symmetry is
conserved in the considered scenario. Notice that the quartic couplings
$\lambda_{7}$ and $\lambda_{13}$ are crucial for generating neutrino masses
through out the two-loop level mechanism. The nonzero scalar VEVs are denoted
as $\langle\Phi\rangle=v\sim 246\leavevmode\nobreak\ \mathrm{GeV}$ and
$\langle\sigma\rangle=v_{\sigma}$. Also, the term $\mu_{sb}^{2}$ that breaks
the $\mathrm{U(1)_{X}}$ symmetry softly is incorporated in order to make the
CP-odd scalar associated with $\sigma$ massive. Once the scalars $\Phi$ and
$\sigma$ acquire VEVs and the EW symmetry is broken, the model will have nine
physical scalars – two charged scalars $(\eta^{\pm})$, four neutral CP-even
scalars ($\eta_{R},\,\varphi_{R},\,h_{1}$ and $h_{2}$) and three neutral CP-
odd scalars ($\eta_{I},\,\varphi_{I}$ and $\sigma_{I}$). It is to be noted
that the scalars associated with $\eta$ and $\varphi$ mix only among
themselves and so do the scalar fields coming from $\Phi$ and $\sigma$111Let
us note that similar scalar potentials can be found in scenarios that explain
Dirac neutrino masses at the one-loop level [20]. A detailed investigation of
such scenarios, following the lines of the present work, will be dedicated in
a subsequent paper..

###  2.3 Scalar mass spectrum

After employing the minimization conditions along $v$ and $v_{\sigma}$, the
scalar mass matrices can be written as

$\displaystyle M^{2}_{\eta^{\pm}}=\begin{pmatrix}0&0&0&0\\\ 0&0&0&0\\\
0&0&\mu_{\eta}^{2}+\dfrac{v_{\sigma}^{2}\lambda_{10}+v^{2}\lambda_{5}}{2}&0\\\
0&0&0&\mu_{\eta}^{2}+\dfrac{v_{\sigma}^{2}\lambda_{10}+v^{2}\lambda_{5}}{2}\end{pmatrix}\,,$
(2.2) $\displaystyle M^{2}_{\mathrm{CP-
even}}=\begin{pmatrix}2v^{2}\lambda_{1}&vv_{\sigma}\lambda_{8}&0&0\\\
vv_{\sigma}\lambda_{8}&2v_{\sigma}^{2}\lambda_{4}&0&0\\\
0&0&\mu_{\eta}^{2}+\dfrac{v_{\sigma}^{2}\lambda_{10}+v^{2}\lambda_{5}+v^{2}\lambda_{6}+\lambda_{7}v^{2}}{2}&0\\\
0&0&0&\mu_{\varphi}^{2}+\dfrac{v_{\sigma}^{2}\lambda_{12}+v_{\sigma}^{2}\lambda_{13}+v^{2}\lambda_{9}}{2}\end{pmatrix}\,,$
$\displaystyle M^{2}_{\mathrm{CP-odd}}=\begin{pmatrix}0&0&0&0\\\
0&4\mu_{sb}^{2}&0&0\\\
0&0&\mu^{2}_{\eta}+\dfrac{v_{\sigma}^{2}\lambda_{10}+v^{2}\lambda_{5}+v^{2}\lambda_{6}-v^{2}\lambda_{7}}{2}&0\\\
0&0&0&\mu_{\varphi}^{2}+\dfrac{v_{\sigma}^{2}\lambda_{12}-v_{\sigma}^{2}\lambda_{13}+v^{2}\lambda_{9}}{2}\end{pmatrix}\,,$

where the zero-eigenvalues of the charged and CP-odd matrices are absorbed by
the longitudinal components of the $W^{\pm}$ and $\mathrm{Z^{0}}$ gauge
bosons. To avoid constraints coming from the Higgs sector, we work closely to
the alignment limit, such that $vv_{\sigma}\lambda_{8}\ll 1$. The above
matrices can be analytically diagonalised, which results in the mass
eigenvalues

$\displaystyle\mathrm{Charged\leavevmode\nobreak\
scalars}:\begin{cases}m_{\eta^{\pm}}^{2}=\mu_{\eta}^{2}+\dfrac{1}{2}\left(\lambda_{5}v^{2}+\lambda_{10}v_{\sigma}^{2}\right)\,.\end{cases}$
(2.3) $\displaystyle\mathrm{CP-even\leavevmode\nobreak\
scalars}:\begin{cases}&m_{h_{1}}^{2}=\lambda_{1}v^{2}+\lambda_{4}v_{\sigma}^{2}-\sqrt{(\lambda_{1}v^{2}-\lambda_{4}v_{\sigma}^{2})^{2}+\lambda_{8}^{2}v^{2}v_{\sigma}^{2}}\,,\\\\[6.00006pt]
&m_{h_{2}}^{2}=\lambda_{1}v^{2}+\lambda_{4}v_{\sigma}^{2}+\sqrt{(\lambda_{1}v^{2}-\lambda_{4}v_{\sigma}^{2})^{2}+\lambda_{8}^{2}v^{2}v_{\sigma}^{2}}\,,\\\\[6.00006pt]
&m_{\eta_{R}}^{2}=\mu_{\eta}^{2}+\dfrac{1}{2}\left(v^{2}(\lambda_{5}+\lambda_{6}+\lambda_{7})+\lambda_{10}v_{\sigma}^{2}\right)\,,\\\\[6.00006pt]
&m_{\varphi_{R}}^{2}=\mu_{\varphi}^{2}+\dfrac{1}{2}(\lambda_{9}v^{2}+v_{\sigma}^{2}\left(\lambda_{12}+\lambda_{13}\right))\,.\end{cases}$
$\displaystyle\mathrm{CP-odd\leavevmode\nobreak\
scalars}:\begin{cases}&m_{\sigma_{I}}^{2}=4\mu_{sb}^{2}\,,\\\\[6.00006pt]
&m_{\eta_{I}}^{2}=\mu_{\eta}^{2}+\dfrac{1}{2}\left(v^{2}(\lambda_{5}+\lambda_{6}-\lambda_{7})+\lambda_{10}v_{\sigma}^{2}\right)\,,\\\\[6.00006pt]
&m_{\varphi_{I}}^{2}=\mu_{\varphi}^{2}+\dfrac{1}{2}(\lambda_{9}v^{2}+v_{\sigma}^{2}\left(\lambda_{12}-\lambda_{13}\right))\,.\end{cases}$

Here, the lightest scalar $m_{h_{1}}$ is taken to be the SM Higgs boson of
mass $\sim 125$ GeV. The above expressions can be inverted to express 7 out of
the 13 quartic couplings in terms of the physical scalar masses. Thus, the
physical masses can be taken as the input parameters in the numerical
analysis. The relations that are used in this work are,

$\displaystyle\lambda_{1}$
$\displaystyle=\frac{m^{2}_{h_{1}}-\sqrt{(m^{4}_{h_{1}}-2m^{2}_{h_{2}}m^{2}_{h_{2}}+m^{4}_{h_{2}}-4\lambda_{8}^{2}v_{1}^{2}v_{\sigma}^{2})}+m^{2}_{h_{2}}}{4v_{1}^{2}}$
(2.4) $\displaystyle\lambda_{4}$
$\displaystyle=\frac{m^{2}_{h_{1}}+\sqrt{\left(m^{4}_{h_{1}}-2m^{2}_{h_{1}}m^{2}_{h_{2}}+m^{4}_{h_{2}}-4\lambda_{8}^{2}v_{1}^{2}v_{\sigma}^{2}\right)}+m^{2}_{h_{2}}}{4v_{\sigma}^{2}}$
$\displaystyle\lambda_{5}$
$\displaystyle=\frac{-2\mu_{\eta}^{2}+2m_{\eta^{\pm}}^{2}-\lambda_{10}v_{\sigma}^{2}}{v_{1}^{2}}$
$\displaystyle\lambda_{6}$
$\displaystyle=\frac{m_{\eta_{I}}^{2}+m_{\eta_{R}}^{2}-2m_{\eta^{\pm}}^{2}}{v_{1}^{2}}$
$\displaystyle\lambda_{7}$
$\displaystyle=\frac{-m_{\eta_{I}}^{2}+m_{\eta_{R}}^{2}}{v_{1}^{2}}$
$\displaystyle\lambda_{9}$
$\displaystyle=\frac{-2\mu_{sb}^{2}+m_{\eta_{I}}^{2}+m_{\eta_{R}}^{2}-v_{\sigma}^{2}\lambda_{12}}{v_{1}^{2}}$
$\displaystyle\lambda_{13}$
$\displaystyle=\frac{m_{\left(\varphi\right)_{R}}^{2}-m_{\left(\varphi\right)_{I}}^{2}}{2v_{\sigma}^{2}}$

###  2.4 Implications for collider physics and Dark Matter phenomenology

Finally, to close this section we wish to provide a brief discussion of the
key phenomenological implications of the multi-scalar sector in the considered
model. It is expected that the presence of electrically charged scalars in the
inert doublet will provide an extra contribution to the Higgs diphoton decay
rate $h\to\gamma\gamma$, which will yield a deviation of the Higgs diphoton
signal strength from the SM expectation. However, that deviation is expected
to occur in the experimentally allowed range in a large region of the model
parameter space. Besides, given that the singlet scalar field acquires a VEV
at the TeV scale, its mixing with the CP-even part of the active doublet is
suppressed implying that the couplings of the $125$ GeV SM like Higgs boson
will be very close to the SM expectation realising the so-called Higgs
alignment limit. Thus, such a scenario appears to be consistent with the SM in
such a limit, at least, at the leading order (LO) level.

It is generally expected that the triple Higgs coupling can be affected
compared to its SM value in multi-scalar extensions of the SM. In this study,
the trilinear Higgs coupling, $\lambda_{hhh}$, deviates from the SM prediction
through a mixing with the $\sigma$ field. In particular, at LO, a simple
analytical expression is obtained and reads as

$\displaystyle\lambda^{\mathrm{LO}}_{hhh}=6\cos^{3}\alpha_{h}v\lambda_{1}+6\sin^{3}\alpha_{h}v_{\sigma}\lambda_{4}+3\cos\alpha_{h}\sin^{2}\alpha_{h}v\lambda_{8}+3\cos^{2}\alpha_{h}\sin\alpha_{h}v_{\sigma}\lambda_{8}\,,$
(2.5)

where $\alpha_{h}$ is the CP-even mixing angle between $h_{1}$ (Higgs) and
$h_{2}$. It is not hard to notice that in the limit of $\alpha_{h}\rightarrow
0$, the tree-level SM Higgs coupling
$\lambda_{hhh}^{\mathrm{SM}}=6v\lambda_{1}$ is successfully reproduced as we
effectively work very close to the alignment limit. Hence, the potential
corrections to $\lambda^{\mathrm{LO}}_{hhh}$ are very suppressed and we expect
no significant deviations from the SM expectation. On the other hand, previous
works [21, 22, 23] have suggested that one-loop corrections to the
$\lambda_{hhh}$ can actually be quite sizable in some BSM extensions,
resulting, in some cases, in an increase of order $\mathcal{O}(50\%)$ for
certain combinations coupling/mass parameters. In our case, next-to-leading
order (NLO) contributions will be dominant and represent the bulk of the new
physics contribution to $\lambda_{hhh}$. So, for a correct understanding of
the impact on $\lambda_{hhh}$, one must take these contributions into account.
In this regard, we have computed all one-loop contributions to the triple
Higgs coupling using the effective potential approach as described in Ref.
[24]. Note that here we are considering the limit in which the external
momentum vanishes, i.e. $p^{2}\to 0$, when the kinematical suppression from
$p^{2}\neq 0$ corrections is less or similar to the one produced by the top
quark [25]. Therefore, we have considered all possible contributions coming
from all physical scalars of the model ($\eta^{\pm}$, $h_{1}$, $h_{2}$,
$\eta_{R}$, $\varphi_{R}$, $\sigma_{I}$, $\eta_{I}$ and $\varphi_{I}$), as
well as the loop diagram involving the virtual exchange of the top quark. The
full one-loop formulas are listed in Appendix A. In what follows, in our
numerical analysis we will take into account the existing constraints on
$\lambda_{hhh}$ available from the LHC [26, 27]. Future collider measurements
[28, 29] can further probe possible deviations of $\lambda_{hhh}$ from the SM
value.

Another interesting implication the considered model is the presence of
potential DM candidates. First, the scalar one can be identified with the
lightest among the $\mathop{\mathrm{R}e}\varphi$,
$\mathop{\mathrm{I}m}\varphi$, $\mathop{\mathrm{R}e}\eta$ and
$\mathop{\mathrm{I}m}\eta$ fields. Assuming for simplicity that the
cosmological DM is dominated by such a scalar DM candidate, it would
annihilate mainly into $WW$, $ZZ$, $t\overline{t}$, $b\overline{b}$ and
$h_{1}h_{1}$ channels in the early Universe via a Higgs portal scalar
interaction. These annihilation channels will contribute to the DM relic
density, which can be accommodated for appropriate values of the coupling of
the Higgs portal scalar interaction and of the scalar DM mass, which in most
cases are at the TeV scale. In the low mass region, namely $m_{DM}<m_{h}/2$,
the constraints coming from the invisible Higgs decays should be taken into
account. In our numerical analysis, for simplicity, these decays are
kinematically forbidden. That is, we set the DM mass to be above $m_{h}/2$.
Let us note that in the DM resonant region, i.e. $m_{DM}\sim m_{h}/2$, the
relic density is successfully accounted for. In addition, DM direct detection
implies that the scalar DM candidate would scatter off a nuclear target in a
detector via Higgs boson exchange in the $t$-channel, giving rise to a
constraint on the Higgs portal scalar interaction coupling.

The model may also feature a possible fermionic DM candidate which can be
identified as the lightest among the right-handed Majorana neutrinos
$N_{k,\mathrm{R}}$ ($k=1,2$). Such particles are expected to annihilate into a
pair of active neutrinos via $t$-channel exchanges of the neutral CP-even and
CP-odd components of the inert scalar doublet $\eta$. The DM relic density
constraints can be successfully accommodated for an appropriate region of
parameter space. A detailed analysis of the implications of our model for DM
studies is beyond the scope of the present work and is left for future
investigations.

##  3 Neutrino mass and lepton flavor violation

Once $\Phi$ and $\sigma$ develop VEVs, a non-trivial contribution to the
neutrino mass matrix is generated at two-loop level, whose diagram is
illustrated in Fig. 1. The active light neutrino mass matrix in this case is
given as,

$\displaystyle\left(M_{\nu}\right)_{ij}=$
$\displaystyle\frac{\left(y_{N}\right)_{in}\left(y_{\Omega}^{\ast}\right)_{nk}\left(y_{\Omega}^{{\dagger}}\right)_{kr}\left(y_{N}^{T}\right)_{rj}m_{\Psi_{k}}}{4(4\pi)^{4}}$
(3.1)
$\displaystyle\int_{0}^{1}d\alpha\int_{0}^{1-\alpha}d\beta\frac{1}{\alpha(1-\alpha)}\Biggl{[}G\left(m_{\Psi_{k}}^{2},m_{RR}^{2},m_{RI}^{2}\right)-G\left(m_{\Psi_{k}}^{2},m_{IR}^{2},m_{II}^{2}\right)\Biggr{]}\,,$

where $n,k,r=1,2$. Here, the loop integral $I$ can be written as [30]

$\displaystyle G(x^{2},y^{2},z^{2})\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{x^{2}y^{2}\log\left(\displaystyle\frac{y^{2}}{x^{2}}\right)+y^{2}z^{2}\log\left(\displaystyle\frac{z^{2}}{y^{2}}\right)+z^{2}x^{2}\log\left(\displaystyle\frac{x^{2}}{z^{2}}\right)}{(x^{2}-y^{2})(x^{2}-z^{2})}\,,$
$\displaystyle m_{ab}^{2}\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\frac{\beta m_{\left(\eta\right)_{a}}^{2}+\alpha
m_{\left(\varphi\right)_{b}}^{2}}{\alpha(1-\alpha)}\quad(a,b=R:\mathrm{or}:I)\,.$
(3.2)

Note that the lightest active neutrino is massless which implies that there is
only one Majorana phase in the PMNS matrix. A convenient way to incorporate
the bounds from the neutrino oscillation data is to express the Yukawa
coupling matrix using the Casas-Ibarra parametrization. For the model
considered in this paper, one can express $y_{N}$ as,

$y_{N}=U_{\mathrm{PMNS}}^{*}\sqrt{M_{\nu}^{\mathrm{diag}}}R^{T}A^{-1}\,,$
(3.3)

where $U_{\mathrm{P}MNS}^{*}$ is the $3\times 3$ PMNS neutrino mixing matrix,
$M_{\nu}^{\mathrm{diag}}=\leavevmode\nobreak\ \text{diag}(m_{1},m_{2},m_{3})$
is the diagonal neutrino mass matrix, $R$ is a $2\times 3$ matrix that is
given as,

$R=\begin{pmatrix}0&\text{cos}z&-\text{sin}z\\\
0&\text{sin}z&\text{cos}z\end{pmatrix}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{for normal hierarchy (NH) of light neutrino
masses}\leavevmode\nobreak\ \leavevmode\nobreak\ (m_{1}=0)\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}$
$R=\begin{pmatrix}\text{cos}z&-\text{sin}z&0\\\
\text{sin}z&\text{cos}z&0\end{pmatrix}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \text{for inverted hierarchy (IH) of light neutrino
masses}\leavevmode\nobreak\ \leavevmode\nobreak\ (m_{3}=0).$

The matrix $A$ is given by

$A=y_{\Omega}^{*}\Lambda y_{\Omega}^{\dagger}\,,$

where $\Lambda=\begin{pmatrix}\Lambda_{1}&0\\\ 0&\Lambda_{2}\end{pmatrix}$
with,

$\Lambda_{k}=\frac{m_{\Psi_{k}}}{4(4\pi)^{4}}\int_{0}^{1}d\alpha\int_{0}^{1-\alpha}d\beta\frac{1}{\alpha(1-\alpha)}\Biggl{[}G\left(m_{\Psi_{k}}^{2},m_{RR}^{2},m_{RI}^{2}\right)-G\left(m_{\Psi_{k}}^{2},m_{IR}^{2},m_{II}^{2}\right)\Biggr{]}\,.$

Under this parameterization, the experimentally observed neutrino mass
differences and the PMNS mixing matrix can be given as input in the numerical
scan, with the Yukawa coupling being found as output. Since we do not have
control over the magnitude of the Yukawa couplings, we only allow solutions
which respect perturbativity, i.e.
$\absolutevalue{Y_{N}},\absolutevalue{Y_{\Omega}}<\sqrt{4\pi}$. For the
purpose of numerical scans, we only consider the normal ordering scenario for
the neutrino masses (that is, $m_{1}=0$).

The presence of the Majorana fermions can induce LFV decays, such as
$\mu\rightarrow e\gamma$, which are strongly constrained by experiment. In our
model, such decays are mediated at one-loop level via virtual exchanges of the
neutral fermions and the charged scalars. The branching fraction for the two-
body decay process $\ell_{i}\rightarrow\ell_{j}\gamma$, where $i=e,\mu,\tau$
is given as [31, 32, 33, 34]

$\mathrm{BR}\left(l_{i}\rightarrow
l_{j}\gamma\right)=\frac{3\left(4\pi\right)^{3}\alpha_{\mathrm{em}}}{4G_{F}^{2}}\left|\frac{x_{is}^{\left(\nu\right)}x_{js}^{\left(\nu\right)}}{2\left(4\pi\right)^{2}m_{\eta^{\pm}}^{2}}F\left(\frac{m_{N_{sR}}^{2}}{m_{\eta^{\pm}}^{2}}\right)\right|^{2}\mathrm{BR}\left(l_{i}\rightarrow
l_{j}\nu_{i}\overline{\nu_{j}}\right)\,,$ (3.4)

with $s=1,2$. Here, $\alpha_{\mathrm{em}}=1/137$ is the fine-structure
constant,
$x_{is}^{\left(\nu\right)}=\sum_{k=1}^{3}(y_{N})_{ks}^{\left(\nu\right)}\left(V_{lL}^{\dagger}\right)_{ik}$,
with $V$ being the left-handed charged lepton mixing matrix,
$G_{F}=1.166364\times 10^{-5}\leavevmode\nobreak\ \mathrm{GeV^{-2}}$ is the
Fermi constant, which in our case is the identity matrix, $m_{\varphi^{\pm}}$
are the masses of the charged scalar components of the $\mathrm{SU(2)_{L}}$
inert doublet $\varphi$, and $m_{N_{sR}}$ ($s=1,2$) correspond to the masses
of the right-handed Majorana neutrinos $N_{sR}$, which due to the charge
assignment, are generated at one-loop level. The loop function $F$ is written
as

$F\left(x\right)=\frac{1-6x+3x^{2}+2x^{3}-6x^{2}\ln x}{6\left(1-x\right)}\,.$
(3.5)

The most stringent bounds for LFV come from muon decay measurements, namely,
$\mu\rightarrow e\gamma$. Current experimental results put an upper bound on
the branching ratio, which reads as $\mathrm{BR}\left(\mu\rightarrow
e\gamma\right)<4.2\times 10^{-13}$ [35]. The model can also give rise to
$\mu\rightarrow e$ conversion in the atomic nucleus. But the existing bounds
are not as restrictive as the ones from $\mathrm{BR}\left(\mu\rightarrow
e\gamma\right)$. However, among various planned experiments, the
$\mu\rightarrow e$ conversion in the Al nuclei has the best projected
sensitivity of $\sim 10^{-17}$ [36]. The branching ratio for this process can
approximately be expressed as [34],

$\mathrm{CR}(\mu\,\mathrm{Al}\rightarrow
e\,\mathrm{Al})\approx\frac{1}{350}\mathrm{BR}\left(\mu\rightarrow
e\gamma\right)\,.$ (3.6)

Figure 2: The predictions for $\mathrm{BR}\left(\mu\rightarrow e\gamma\right)$
(green points) and $\mathrm{CR}(\mu\,\mathrm{Al}\rightarrow e\,\mathrm{Al})$
(red points) as functions of Tr$[y_{N}^{\dagger}y_{N}]$. All the points shown
are within the existing experimental bounds marked by the horizontal dashed
line. The cyan points correspond to the region within the projected
sensitivity for $\mu\,\mathrm{Al}\rightarrow e\,\mathrm{Al}$ conversion.

In this regards, in Fig. 2 shows the predictions for
$\mathrm{BR}\left(\mu\rightarrow e\gamma\right)$ (green points) and
$\mathrm{CR}(\mu\,Al\rightarrow e\,Al)$ (red points) as functions of
Tr$[y_{N}^{\dagger}y_{N}]$. Only the values that are within the existing
bounds are shown, whose limit is indicated by a dashed horizontal line in the
plot. The cyan points correspond to the region within the projected
sensitivity for $\mu\,\mathrm{Al}\rightarrow e\,\mathrm{Al}$ conversion.

For all the points shown in Fig. 2, neutrino oscillation data is respected,
with the relevant experimental parameters (namely, neutrino mass differences
and the PMNS mixing angles) being given as input and allowed to vary within
their respective $3\sigma$ uncertainties [37, 38, 39]. The Majorana phase is
also varied in the range $[0,\pi]$ and $y_{\Omega}$ is taken to be diagonal
where the entries are allowed to vary in the range $[0,1]$. Here, $z$ is
defined as $x_{1}+y_{1}i$ with $x_{1}$ and $y_{1}$ being varied in the range
$[-2,+2]$, the fermions $\Psi_{kR}$ and $N_{kR}$ are taken to be degenerate,
with their masses varied in the ranges $[0.5,5]$ TeV and $[0.1,60]$ GeV,
respectively. For the scalar potential parameters, we employ the inverted
equations in (2.4) to allow for the physical masses to be given as input. In
this regard, we have considered the following parameter ranges:

  * •

$m_{\eta_{I}},m_{\varphi_{I}}\rightarrow[0.5,3.2]$ TeV,
$m_{\eta_{R}}\rightarrow[m_{\eta_{I}},m_{\eta_{I}}+5\,\textrm{GeV}]$,
$m_{\varphi_{R}}\rightarrow[m_{\varphi_{I}},m_{\varphi_{I}}+5\,\textrm{GeV}]$,
$m_{\eta^{\pm}}\rightarrow[m_{\eta_{I}}-5\,\textrm{GeV},m_{\eta_{I}}]$

  * •

$v_{1}=246$ GeV, $m_{h_{2}}\rightarrow[0.5,5]$ TeV,
$v_{\sigma}\rightarrow[5,30]$ TeV,

  * •

$\lambda_{2,5,8,9,10,11,12}\rightarrow[0,1]$. The remaining $\lambda$’s are
treated as outputs in the inversion procedure. Notice that we have inverted
the expressions for $\lambda_{5}$ and $\lambda_{9}$, which can be solved such
that $\mu_{\eta}^{2}$ and $\mu_{\varphi}^{2}$ are instead given as output
parameters, where we have taken them to be both greater than $0$. Also, in the
parameter space that has been considered, we have ensured that all the quartic
couplings, $\lambda_{i}>0$, even though this is only a sufficient condition
but not necessary.

##  4 Phase transitions and primordial gravitational waves

One of the more interesting phenomenological consequences of extended scalar
sectors compared to the SM framework is the possibility for the existence of
FOPTs. Indeed, such a phenomenon is not present in the SM (both the EW and QCD
phase transitions are of second order), making it pure beyond-SM driven (BSM)
physics. Such transitions can lead to the generation of a SGWB. With the
detection of GWs at LIGO [40], a new era of multi-messenger astronomy has
arisen, where such events may give us unique perspectives into new physics
beyond the SM. It has been recently proposed that the cosmological phase
transitions between different vacua at finite temperatures (such as those
associated with symmetry breaking in BSM models) may give rise to
gravitational imprints that should be observable in future ground- and space-
based experiments [41, 42]. The SGWB signals can also provide hints into the
nature of the neutrino masses in the context of distinct seesaw mechanisms
[43].

As the universe expands and cools down, its temperature dependent vacuum
configuration evolves giving rise to transitions between different phases. If
such transitions are sufficiently strongly first order, they can lead to the
generation of fast expanding vacuum bubbles whose dynamics may have left a
visible imprint in the form of a SGWB [44, 45]. The characteristics of the
phase transition such as its strength and inverse duration depend on the
underlying particle physics model described at finite temperatures and are
highly sensitive, in particular, to the properties of the potential barrier
between the phases and to the critical temperature of the transition (see e.g.
[46, 41, 47, 48, 49, 50] and references therein). The phenomenon of strong
FOPTs is quite generic for multi-scalar extensions of the SM (in particular,
for those that undergo several symmetry breaking steps) such as the one we
discuss in this article. Therefore, a SGWB from cosmological FOPTs must be
regarded as an important source of phenomenological information about multi-
scalar extensions of the SM which may be complementary to constraints coming
from collider searches.

In order to fully classify and study the SGWB emergent in our model, there is
a total of five main ingredients. The first one is inherently dependent on the
considered model structure – the finite temperature effective potential [51,
52, 53]. Generically, it can be written as [54, 55]

$V_{\mathrm{eff}}=V^{(0)}+V_{\mathrm{CW}}^{(1)}+V_{T\neq
0}+V_{\mathrm{CT}}\,,$ (4.1)

where the first term, $V^{(0)}$, is the tree-level scalar potential given in
Eq. (2.2). The second term is the Coleman-Weinberg one-loop correction,

$V^{(1)}_{\mathrm{CW}}=\sum_{a}(-1)^{F}n_{a}\frac{m_{a}^{2}(\phi_{b})}{64\pi^{2}}\Bigg{[}\ln(\frac{m_{a}^{2}(\phi_{b})}{\mu^{2}})-C_{a}\Bigg{]}\,,$
(4.2)

where $F$ is 0 for bosons and 1 for fermions, $m_{a}(\phi_{b})$ is the
$\phi_{a}$-field dependent mass of the particle $a$, $n_{a}$ is the number of
degrees of freedom for each particle $a$ and $\mu$ is a renormalisation scale.
The particle degrees of freedom can be computed as $(-1)^{2s}QN(2s+1)$, where
$s$ is the spin of the particle, $Q=1(2)$ for neutral (charged) particles, and
$N$ is the number of colours. In this work, we consider the
$\overline{\mathrm{MS}}$ renormalisation scheme such that we have $C_{a}=3/2$
for scalars, fermions and longitudinally polarized gauge bosons, and
$C_{a}=1/2$ for transverse gauge bosons. Additionally, the renormalisation
scale is fixed as $\mu=\prod_{i=1}^{n}m_{\phi_{i}}^{1/n}$ where $i$ runs over
all BSM scalars in the model. Since we allow to vary the masses of the scalar
fields, $\mu$ will be different for each sampled point.

The third term encodes thermal corrections to the potential,

$V_{T\neq
0}=\frac{T^{4}}{2\pi^{2}}\sum_{a}(-1)^{F}n_{a}J_{a}\quantity[\frac{m_{a}^{2}(\phi_{a})}{T^{2}}]\,,$
(4.3)

where $T$ is the temperature, $J_{0}$ (for bosons) and $J_{1}$ (for fermions)
are the standard thermal integrals found e.g. in Ref. [54]. At the leading
$(m/T)^{2}$ order in thermal expansions of $J_{a}$, Eq. (4.3) can be
approximated as

$\Delta
V(T)\simeq\frac{T^{4}}{2\pi^{2}}\quantity(\Tr{M^{2}(\phi_{\alpha})}+\sum_{i=W,Z,\gamma}n_{i}m_{i}^{2}(\phi_{\alpha})+\sum_{i=f_{i}}\frac{n_{i}}{2}m_{i}^{2}(\phi_{\alpha}))\,,$
(4.4)

where in the last sum, all the fermions of the model are included. Here, $M$
is the field-dependent Hessian matrix in the classical field approximation
where all fields are fixed to their classical (background) components, that is

$\Phi\rightarrow\begin{pmatrix}0\\\
v_{\Phi}/\sqrt{2}\end{pmatrix}\,,\quad\sigma\rightarrow
v_{\sigma}\,,\quad\eta\rightarrow\begin{pmatrix}0\\\
0\end{pmatrix}\,,\quad\varphi\rightarrow 0\,,$ (4.5)

and the degrees of freedom, $n_{i}$, for different states are as follows [54]:
for quarks $n_{a}=12$, for charged leptons $n_{a}=4$, for neutral leptons
$n_{a}=2$, for the gauge bosons $n_{W}=6$, $n_{Z}=3$, $n_{\gamma}=2$ and for
the longitudinal component of the photon $n_{\gamma_{L}}=1$. Additional
higher-order corrections must also be included in the form of Daisy (ring)
diagrams [51, 56, 57, 58]. These can be parameterized as thermal corrections
to the quadratic mass terms of the bare potential and read as

$\mu^{2}(T)\rightarrow\mu^{2}+c_{\phi}(T^{2})\equiv\mu^{2}+\frac{\delta^{2}\Delta
V(T,\phi_{a})}{\delta\phi_{a}}\,.$ (4.6)

For the considered model, we have

$\displaystyle
c_{\phi}(T^{2})=\frac{T^{2}}{24}\Big{[}3g^{2}+\frac{3}{2}(g^{2}+g_{Y}^{2})+12(y_{b}^{2}+y_{t}^{2})+4y_{\tau}^{2}+12\lambda_{1}+4\lambda_{5}+2(\lambda_{6}+\lambda_{8}+\lambda_{9})\Big{]}$
(4.7) $\displaystyle
c_{\sigma}(T^{2})=\frac{T^{2}}{24}(4(\lambda_{10}+\lambda_{8})+2\lambda_{12}+8\lambda_{4})\,,$

where $g$ ($g_{Y}$) are the SU(2) (U(1)) gauge couplings and $y_{i}$ are the
Yukawa couplings for each of the SM-like third-generation fermions. The gauge
boson masses also gain thermal corrections and, since there are no new gauge
bosons in the model, the formula is identical to what is found in the previous
literature (see e.g. Ref. [54]).

The last relevant part of the effective potential in Eq. (4.1) is the counter-
term potential $V_{\mathrm{CT}}$. The counter-terms are defined such that the
minima of the tree-level potential remains fixed at zero temperature (i.e. the
tree-level masses remain the same at one-loop level), which is guaranteed by
imposing that the first and second derivatives of the counter-term potential
match the first and second derivatives of the Coleman-Weinberg one-loop
potential as

$\left\langle\frac{\partial
V_{\mathrm{CT}}}{\partial\phi_{i}}\right\rangle=\left\langle-\frac{\partial
V^{(1)}_{\mathrm{C}W}}{\partial\phi_{i}}\right\rangle\,,\qquad\left\langle\frac{\partial^{2}V_{\mathrm{CT}}}{\partial\phi_{i}\partial\phi_{j}}\right\rangle=\left\langle-\frac{\partial^{2}V^{(1)}_{\mathrm{C}W}}{\partial\phi_{i}\partial\phi_{j}}\right\rangle\,,\qquad\text{with}\qquad\phi_{i}=\left\\{h_{1},h_{2}\right\\}\,.$
(4.8)

In our case, this results in the following conditions

$\displaystyle\delta\mu_{\phi}=\frac{3}{2v_{\Phi}}\frac{\partial
V^{(1)}_{\mathrm{CW}}}{\partial
h_{1}}-\frac{1}{2}\frac{\partial^{2}V^{(1)}_{\mathrm{CW}}}{\partial
h^{2}_{1}}+\frac{v_{\sigma}}{2v_{\Phi}}\frac{\partial^{2}V^{(1)}_{\mathrm{CW}}}{\partial
h_{2}\partial h_{1}}\,,\quad\quad$
$\displaystyle\delta\mu_{sb}=\frac{3}{4v_{\sigma}}\frac{\partial
V^{(1)}_{\mathrm{CW}}}{\partial
h_{2}}+\frac{v_{\Phi}}{4v_{\sigma}}\frac{\partial^{2}V^{(1)}_{\mathrm{CW}}}{\partial
h_{2}\partial
h_{1}}-\frac{1}{4}\frac{\partial^{2}V^{(1)}_{\mathrm{CW}}}{\partial
h_{2}^{2}}\,,$ (4.9)
$\displaystyle\delta\lambda_{1}=\frac{1}{2v_{\Phi}^{3}}\quantity(v_{\Phi}\frac{\partial^{2}V^{(1)}_{\mathrm{CW}}}{\partial
h_{1}^{2}}-\frac{\partial V^{(1)}_{\mathrm{CW}}}{\partial h_{1}})\,,$
$\displaystyle\delta\lambda_{4}=\frac{1}{2v_{\sigma}^{3}}\quantity(v_{\sigma}\frac{\partial^{2}V^{(1)}_{\mathrm{CW}}}{\partial
h_{2}^{2}}-\frac{\partial V^{(1)}_{\mathrm{CW}}}{\partial h_{2}})\,,$
$\displaystyle\delta\lambda_{8}=\frac{1}{v_{\Phi}v_{\sigma}}\frac{\partial
V^{(1)}_{\mathrm{CW}}}{\partial h_{2}\partial h_{1}}\,,$

while for other parameters, we set $\delta\lambda_{i}=\delta\mu_{i}=0$.

The remaining four ingredients determine the FOPT dynamics, which include [45]

  * •

The inverse time-scale of the phase transition, $\beta$. This parameter
essentially determines the bubble nucleation rate. Normalising it
conventionally to the Hubble expansion rate, $H$, one writes

$\frac{\beta}{H}=T_{n}\frac{\partial}{\partial
T}\quantity(\frac{S(V_{\mathrm{eff}})}{T})\Bigr{|}_{\begin{subarray}{c}T=T_{n}\end{subarray}}\,,$
(4.10)

where $S(V_{\mathrm{eff}})$ is the Euclidean action which depends on the
structure of the effective potential in Eq. (4.1).

  * •

Nucleation temperature, $T_{n}$. It represents the temperature of the Universe
directly after the phase transition, and can be computed by requiring that the
probability for a single bubble nucleation per horizon volume is equal to one
[59].

  * •

Phase transition strength $\alpha$, which is related to the released latent
heat during the phase transition. It is found in terms of the difference
between the effective potential values before and after the transition via the
following relation [60, 61]

$\alpha=\frac{30}{g_{*}\pi^{2}T_{n}^{2}}\quantity[V_{\mathrm{eff},i}-V_{\mathrm{eff},f}-\frac{T}{4}\quantity(\frac{\partial
V_{\mathrm{eff},i}}{\partial T}-\frac{\partial V_{\mathrm{eff},f}}{\partial
T})]\,,$ (4.11)

being also normalised to the radiation energy density of the Universe at the
nucleation temperature. Here, $g_{*}$ is the effective number of relativistic
degrees of freedom in the cosmic plasma. The effective potential is evaluated
at initial (metastable) phase, $V_{\mathrm{eff},i}$, and at the final (stable)
phase, $V_{\mathrm{eff},f}$.

  * •

Bubble wall velocity, $v_{b}$. This quantity represents the speed of the new
phase at the interface of the nucleating bubble. In the current work we simply
assume the case of supersonic detonations with $v_{b}=0.95$. An improved
analysis using the recently developed methods for the determination of the
speed of sound [62] and wall velocity [63] is left for future work.

The authors are also aware of recently developed techniques accounting for a
great improvement in the calculation of thermal effects, for instance, using
dimensional reduction [64, 65, 66, 67, 68]. Improved studies incorporating
such methods are also left for future dedicated work.

The power spectrum associated to the production of a SGWB can arise from three
distinct sources,

$h^{2}\Omega_{\mathrm{GW}}=h^{2}\Omega_{\mathrm{BC}}+h^{2}\Omega_{\mathrm{SW}}+h^{2}\Omega_{\mathrm{TURB}}\,.$
(4.12)

Here, $h^{2}\Omega_{BC}$ is the GW power spectrum produced via collision of
vacuum bubbles. It offers a leading contribution when such bubbles “run-away”,
that is, the friction created by the thermal plasma at the boundaries is not
enough to stop the acceleration of the wall [45]. This effect is typically
subleading and can usually be neglected [69], unless one is considering
supercooled FOPTs [70, 71, 72]. The second term, $h^{2}\Omega_{\mathrm{SW}}$
is the power spectrum of GWs emerging from sound-waves triggered once the
bubble wall sweeps through the plasma at relativistic speeds. In this case, a
vast portion of the released energy is converted to the waves’ propagation in
the plasma that surrounds the walls as sound waves [45]. The last term is
denotes the magnetohydrodynamic turbulence contribution, arising from non-
linear effects of the sound waves the propagate on the plasma. The
determination of such a contribution is still uncertain. Nonetheless, it is
expected to be less relevant for the position of the SGWB peak amplitude and
frequency in comparison to the sound waves component for supersonic
detonations. Semi-analytical expressions describing each term of Eq. (4.12)
with the various parameters of the FOPTs introduced above, can be found in
[45, 73, 74] and references therein.

###  4.1 Numerical results

With a rich scalar potential, one of the immediate questions about the
considered model is if it can accommodate FOPTs and if they can be strong
enough to leave visible imprints in the form of a SGWB. To address this
question we have implemented the thermal effective potential of the scotogenic
two-loop neutrino mass model into the CosmoTransitions package [75] enabling
us to numerically compute the bounce action and hence the probability for
bubble nucleation at any given temperature. Among the key parameters that
control the dynamics of the phase transitions, $\beta/H$ depends on the
derivative of the Euclidean effective action, and hence can exhibit numerical
instabilities if the action is not numerically well-behaved. To minimize the
corresponding uncertainties, we perform an interpolation of the action on a
point-by-point basis. The technique for the smoothing of the Euclidean action
is described in detail in some of the author’s previous work [48].
Furthermore, for the scan, we have extended the ranges of various parameters
to be more inclusive. Namely, we now perform a logarithmic scan over the
quartic couplings in the ranges of $[10^{-8},\,4\pi]$ and the masses of the
BSM scalars to be above $200\leavevmode\nobreak\ \mathrm{GeV}$. Constraints
from neutrino physics and LFV observables are also imposed in the data points.

First, we show in Fig. 3 the GW peak signal,
$h^{2}\Omega_{\mathrm{GW}}^{\mathrm{peak}}$, as a function of the peak
frequency, $f_{\mathrm{peak}}$ in Hz; both axes are shown in logarithmic
scale. Focusing first on the panels (a), (b) and (c) where we display in the
colour axis the relevant GW observables for the dynamics of the signal,
namely, the phase transition strength (a), the inverse time-scale of the phase
transition (b), and the nucleation temperature (c). First, we begin by noting
that various points in the parameter space can give rise to potentially
detectable SGWB signals in different planned GW interferometers. The latter
include the Laser Interferometer Space Antenna (LISA) [76, 74, 77], the Big
Bang Observer (BBO) [78] and at the Deci-hertz Interferometer Gravitational
wave Observatory (DECIGO) [79], indicating that the gravitational channel may
offer a complementary approach to standard colliders channels in probing our
model.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Peak amplitude of the GW signal, $h^{2}\Omega_{\mathrm{GW}}$ as a
function of the peak frequency, $f_{\mathrm{peak}}$, in Hz. In colour, we have
the logarithm of the phase transition strength $\alpha$ (a), the logarithm of
the inverse time-scale of the phase transitions normalized to the Hubble
parameter (b), the nucleation temperature in GeV (c), the VEV of the singlet
$\sigma$ field in TeV (d), the mass of the second heaviest CP-even scalar (e)
and logarithm of the branching ratio of $\mu\rightarrow e\gamma$ (f). For all
panels, the $x-$axis is displayed in logarithmic scale. All points shown here
are also consistent with neutrino physics constraints as well as with
constraints from the triple Higgs coupling measurements [26, 27].

As expected, a strong correlation between the GW observables and the signal
strength is found. Indeed, higher values of $\alpha$ imply a stronger FOPT
which in turn implies lower values for both the nucleation temperature $T_{n}$
and $\beta/H$. In particular, points which can be observed at LISA fall within
the ranges $\alpha=[0.13,\,0.88]$, $\beta/H=[81.13,\,362.82]$ and
$T_{n}=[99.77,\,520.70]\leavevmode\nobreak\ \mathrm{GeV}$. For the points
within the LISA sensitivity reach, we have also computed the signal-to-noise
ratio (SNR), as showcased in Fig. 4, where we note that for the strongest
transitions, the corresponding SNR is greater than 200, indicating that most
of these points can indeed be probed by the LISA experiment [80, 74, 60]

Figure 4: Signal-to-noise ratio (SNR) plot for a LISA mission profile of 3
years. In the $x$-axis we show the strength of the phase transition ($\alpha$)
and in the $y$-axis the inverse time duration in Hubble units. In the colour
axis, we indicate the magnitude of the phase transition
$h^{2}\Omega^{\mathrm{peak}}_{\mathrm{GW}}$. The coloured isolines
(corresponding to the values 1, 5, 20, 50, 100 and 200) are representative of
the expected values for the SNR. As a general rule of thumb, points with an
SNR greater than 5 can be potentially observed. The grey area corresponds to
the region where the sound-wave component is dominant. All points shown in the
plot correspond to the points potentially detectable by the LISA experiment
(that is, above the red LISA curves of Fig. 3).

Now, focusing our attention on the bottom panels of Fig. 3, where we display
various models parameters in the colour axis, namely, the VEV of the singlet
field $\sigma$ (d), the mass of the second heaviest CP-even scalar $h_{2}$ (e)
and the branching ratio $\mu\rightarrow e\gamma$ (f). Here, we note that there
is no strong correlation between the physical parameters of the model and the
GW observables. Nevertheless, we find that for the larger values of the
singlet VEV and of the mass of $h_{2}$, while allowing for the existence of
FOPTs, they give rise to higher frequency SGWB outside the sensitivity of
future GW experiments. Indeed, for the parameter space regions which are
accessible at LISA, the $\sigma$ VEV is typically below $5\leavevmode\nobreak\
\mathrm{TeV}$ whereas the $h_{2}$ mass is below $3\leavevmode\nobreak\
\mathrm{TeV}$. As for the last panel (f), we do not find any correlation with
$\mu\rightarrow e\gamma$ process. Indeed, we observe that various distinct
possibilities for its branching ratio can be accommodated across different
values of the frequency and the signal strength.

In general, finding certain complementarity between the physical parameters of
the model and SGWB observables is relevant for exploring the model’s parameter
space. In the context of GWs studies, it has been noted in [81, 82, 83, 84]
that potentially large corrections to the trilinear Higgs coupling compared to
the SM expectation can significantly affect the type and the strength of the
cosmological phase transitions. Indeed, strong FOPTs can be driven by large
trilinear scalar couplings as they can play a sizeable effect on the potential
barrier between different vacua.

(a)

(b)

Figure 5: The BSM-to-SM ratio of trilinear Higgs couplings
$\kappa\equiv\lambda^{\mathrm{BSM}}_{hhh}/\lambda^{\mathrm{SM}}_{hhh}$ as a
function of the sine of the CP-even mixing angle $\alpha_{h}$. On the colour
axis we show the logarithm of the GW intensity. On the panel (a), the results
from the full collected dataset are showcased, whereas on (b) a cut on
$h^{2}\Omega^{\mathrm{peak}}_{\mathrm{peak}}$ was imposed, such that only
points which are visible at future GW experiments (LISA, BBO and DECIGO) are
displayed. We note that all constraints from neutrino physics and LFV are
applied here. Regions marked by gray bands are excluded at 95% confidence
level (CL) based on the latest results from the LHC, the blue band region
corresponds to the expected sensitivity at the high-luminosity phase of the
LHC [28] while the red band region corresponds to the expected sensitivity at
the high-energy phase of the LHC [29].

We show in Fig. 5 scatter plots for the BSM-to-SM ratio of the trilinear Higgs
coupling computed for our BSM scenario versus its SM value, i.e.
$\kappa\equiv\lambda^{\mathrm{BSM}}_{hhh}/\lambda^{\mathrm{SM}}_{hhh}$, as a
function of the sine of the CP-even mixing angle, with the logarithm of the GW
signal strength in colour. Here,
$\lambda_{hhh}^{\mathrm{BSM}}=\lambda_{hhh}^{\mathrm{LO}}+\lambda_{hhh}^{\mathrm{NLO}}$
is the full NLO result for the trilinear Higgs coupling in our model and
$\lambda_{hhh}^{\mathrm{SM}}$ corresponds to the SM prediction.

Most of the generated points satisfy $\sin\alpha_{h}<0.07$. Indeed, for such
low values of the scalar mixing angle, Eq. (2.5) asymptotically tends towards
the SM prediction. However, NLO contributions result in a sizeable
contribution to $\lambda_{hhh}$. We can also observe from panel (a) that
points with larger SGWB signal tend to accumulate as one approaches the
alignment limit. However, this does not represent a sufficient condition as we
also find scenarios where the SGWB peak amplitude is well below the
sensitivity ranges of LISA, BBO and DECIGO at low values of the scalar mixing
angle, as can be readily noted in panel (b). Besides, one observes that there
is a part of the parameter space that can already be excluded by the current
data (gray bands) and it can be further reduced by future measurements both at
the high-luminosity and high-energy phases of the LHC. However, the take home
message is that even within those regions there are still several points that
feature FOPTs potentially detectable at future GW facilities as a key feature
of the considered model.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: Scatter plots of various physical parameters of the considered model
(namely, quartic couplings, physical scalar masses and VEVs) with the
logarithm of the GW intensity shown in the colour axis. In particular, we show
$\lambda_{2}$ as a function of $\lambda_{5}$ (a), $\lambda_{2}$ as a function
of $\lambda_{8}$ (b), $\lambda_{4}$ as a function of $\lambda_{10}$ (c), the
mass of $\eta_{R}$ as a function of the $\varphi_{R}$ mass (d), the $h_{2}$
mass as a function of the singlet $\sigma$ VEV (e), and the mass of the
charged scalar $\eta^{\pm}$ as a function of the branching ratio of
$\mu\rightarrow e\gamma$ decay. The masses are given in units of GeV. Here, a
cut on $h^{2}\Omega^{\mathrm{peak}}_{\mathrm{peak}}$ was imposed, such that
only points which are visible at future GW experiments (either at LISA, BBO or
DECIGO) are displayed. In panel (f), the $x$-axis is shown in logarithmic
scale.

In Fig. 6 we show a selection of plots showing the impact of different quartic
couplings and physical masses in the energy density amplitude of the GW
spectra (colour axis). For panels (a), (b) and (c), we find the preferred
ranges of the model’s quartic couplings that lead to observable SGWB. In
particular we one obtains both $\mathcal{O}(1)$ couplings, such as e.g.
$\lambda_{2}$ and $\lambda_{10}$ and much smaller couplings such as
$\lambda_{8}$. Besides, we note that $\lambda_{8}$ is correlated with the size
of the mixing for the CP-even states, and its smallness is related to the fact
that we work very close to the alignment limit. In panels (d), (e) and (f) we
notice that a wide range of physical scalar masses can be successfully
accommodated with the observable GW signals. We highlight the fact that
neither we did find a correlation nor additional restrictions on the quartic
couplings $\lambda_{7}$ and $\lambda_{13}$ involved in the neutrino mass
generation. These parameters can take either small or large values without
modifying the conclusions of this section. It can also be noted that no
apparent correlation between the $\mu\rightarrow e\gamma$ decay rate and the
mass of the charged scalar field is present, exhibiting a distinct flat
distribution where most of the charged scalar masses tend to be below 1 TeV.

##  5 Conclusions

In this work, we have presented a new scotogenic two-loop neutrino mass model
where the SM gauge group is extended with an additional global
$\mathrm{U(1)}\times\mathcal{Z}_{2}$ symmetry. The scalar sector of this
framework contains the standard Higgs doublet, inert scalar doublet, singlet
scalar fields as well as an active scalar singlet. Additionally, the model
features the sector of right-handed Majorana neutrinos which in combination
with the new exotic scalars leads to the generation of neutrino mass and
mixing at two-loop level. The presence of Majorana neutrinos, together with a
physical charged scalar field, also generates highly suppressed contributions
to the tightly constrained LFV observables.

The presence of a rich scalar sector in this model opens up the possibility
for testing it at both collider and GW experiments. In particular, our model
features the generation of cosmological FOPTs in the early Universe that can
lead to the formation of potentially observable SGWB. Focusing primarily on
the phenomenologically favourable Higgs alignment limit, where the mixing
between the CP-even states is small, we have performed a parameter scan of the
model taking into account relevant constraints from LFV processes (such as
$\mu\rightarrow e\gamma$ decay) and those from neutrino physics. In this
article we have identified which regions of the parameter space can
successfully accommodate the existing constraints and giving rise to strong
FOPTs that lead to an observable SGWB at LISA. In addition, we have taken into
account the constraints on the trilinear Higgs coupling at the NLO, and found
a number of points with visible SGWB that can be probed in both collider and
GW experiments in a foreseeable future.

## Acknowledgments

This work was supported by the grants CERN/FIS-PAR/0021/2021, CERN/FIS-
PAR/0019/2021, CERN/FIS-PAR/0024/2021, CERN/FIS-PAR/0025/2021 and PTDC/FIS-
AST/3041/2020. A.E.C.H. acknowledges support by FONDECYT (Chile) under grant
No. 1210378, Milenio-ANID-ICN2019_044 and ANID PIA/APOYO AFB220004. The work
of C.B. was supported by FONDECYT grant No. 11201240. V.K.N. is supported by
ANID-Chile Fondecyt Postdoctoral grant No. 3220005. R.P. is supported in part
by the Swedish Research Council grant, contract number 2016-05996, as well as
by the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No 668679). A.P.M. and J.G.
are supported by the Center for Research and Development in Mathematics and
Applications (CIDMA) through the Portuguese Foundation for Science and
Technology (FCT - Fundação para a Ciência e a Tecnologia), references
UIDB/04106/2020 and UIDP/04106/2020. A.P.M. is also supported by national
funds (OE), through FCT, I.P., in the scope of the framework contract foreseen
in the numbers 4, 5 and 6 of the article 23, of the Decree-Law 57/2016, of
August 29, changed by Law 57/2017, of July 19. J.G. is also directly funded by
FCT through the doctoral program grant with the reference 2021.04527.BD. The
authors also would like to acknowledge the FCT Advanced Computing Project to
provide computational resources via the project CPCA/A00/7395/2020. This work
was partially produced with the support of INCD funded by FCT and FEDER under
the project 01/SAICT/2016 nº 022153. C.B. would like to acknowledge the
hospitality and support from the ICTP through the Associates Programme
(2023-2028).

##  Appendix A One-loop expression for the physical trilinear Higgs coupling

In the zero external momentum approximation, only fields with the mass near or
above the Higgs field provide with leading corrections to the Higgs triple
coupling, which in our model includes the top quark and all BSM scalar fields
at the NLO. The neutral heavy Majorana fields do not contribute to the
amplitude, since they do not couple directly to the Higgs boson before and
after the EW symmetry breaking. Following Ref. [24], the one-loop corrections
involving the top quark reads as

$\displaystyle\lambda_{hhh}^{t}=\frac{1}{32\pi^{2}}\quantity[-36(\mathcal{F}(m_{t},m_{t})-1)\lambda^{(0)}_{tth}\lambda^{(0)}_{tthh}-24\mathcal{F}(m_{t},m_{t},m_{t})(\lambda^{(0)}_{tth})^{3}]\,,$
(A.1)

where $m_{t}\leavevmode\nobreak\ 172.76\leavevmode\nobreak\ \mathrm{GeV}$ is
the top mass. The NLO contribution coming from the scalar fields reads

$\displaystyle\lambda_{hhh}^{\phi}=\frac{1}{32\pi^{2}}\Big{[}$ $\displaystyle
4\mathcal{F}(m_{\eta^{\pm}},m_{\eta^{\pm}},m_{\eta^{\pm}})(\lambda^{(0)}_{h_{1}\eta^{\pm}\eta^{\pm}})^{3}+2\sum_{k=2}^{8}\mathcal{F}(m_{k},m_{k},m_{k})(\lambda^{(0)}_{h_{1}kk})^{3}+$
(A.2) $\displaystyle
6\mathcal{F}(m_{h_{2}},m_{h_{2}},m_{h_{1}})(\lambda^{(0)}_{h_{1},h_{1},h_{2}})^{2}\sum_{k=2}^{3}\lambda_{h_{1}kk}^{(0)}+2\mathcal{F}(m_{\eta_{R}},m_{\eta_{R}},m_{\eta_{R}})(\lambda^{(0)}_{h_{1}\eta_{R}\eta_{R}})^{3}-$
$\displaystyle
3\sum_{k=2}^{8}\quantity(\mathcal{F}(m_{k},m_{k})-1)\lambda_{h_{1}kk}^{(0)}\lambda_{h_{1}h_{1}kk}^{(0)}+6\quantity(\mathcal{F}(m_{\eta^{\pm}},m_{\eta^{\pm}})-1)\lambda_{h_{1}\eta^{\pm}\eta^{\pm}}^{(0)}\lambda_{h_{1}h_{1}\eta^{\pm}\eta^{\pm}}^{(0)}+$
$\displaystyle
6\quantity(\mathcal{F}(m_{h_{1}},m_{h_{2}})-1)\lambda_{h_{1}h_{1}h_{2}}^{(0)}\lambda_{h_{1}h_{1}h_{2}h_{2}}^{(0)}\Big{]}\,,$

where for simplicity of presentation we have $k=1\dots
8\equiv\eta^{\pm},h_{1},h_{2},\eta_{R},\varphi_{R},\sigma_{I},\eta_{I},\varphi_{I}$.
The superscript “(0)” indicates tree-level accuracy. The loop function
$\mathcal{F}$ is defined as [24]

$\mathcal{F}(m_{a_{1}}\dots
m_{a_{N}})=\sum_{x=1}^{N}\frac{m_{a_{x}}^{2}\log\left(\frac{m_{a_{x}}^{2}}{\mu^{2}}\right)}{\Pi_{y\neq
x}\left(m_{a_{x}}^{2}-m_{a_{y}}^{2}\right)}\,,$ (A.3)

where $\mu$ is the renormalisation scale. Here, we have used the same
definition for the scale as in the calculation of the Coleman-Weinberg
potential, i.e. the renormalisation scale is defined as
$\mu=\prod_{i=1}^{n}m_{\phi_{i}}^{1/n}$, where $i$ runs over all BSM scalars
in the model. The various $\lambda$ couplings in the previous expressions are
written as

$\displaystyle\lambda_{thh}^{(0)}=4\cos(\alpha_{h})y_{t}^{2}v\,,$ (A.4)
$\displaystyle\lambda_{tthh}^{(0)}=4\cos^{2}(\alpha_{h})y_{t}^{2}\,,$ (A.5)
$\displaystyle\lambda_{h_{1}\eta^{\pm}\eta^{\pm}}^{(0)}=\sin(\alpha_{h})\lambda_{10}v_{\sigma}+\cos(\alpha_{h})\lambda_{5}v\,,$
(A.6)
$\displaystyle\lambda_{h_{1}h_{1}h_{1}}^{(0)}=6\cos^{3}(\alpha_{h})v\lambda_{1}+6\sin^{3}(\alpha_{h})v_{\sigma}\lambda_{4}+3\cos(\alpha_{h})\sin^{2}(\alpha_{h})v\lambda_{8}+3\cos^{2}(\alpha_{h})\sin(\alpha_{h})v_{\sigma}\lambda_{8}\,,$
(A.7)
$\displaystyle\begin{split}&\lambda_{h_{1}h_{2}h_{2}}^{(0)}=6\text{$\lambda_{1}$}v\sin^{2}(\alpha_{h})\cos(\alpha_{h})+\text{$\lambda_{8}$}v\cos^{3}(\alpha_{h})-2\text{$\lambda_{8}$}v\sin^{2}(\alpha_{h})\cos(\alpha_{h})-\\\
&\hskip
45.00006pt-6\text{$\lambda_{4}$}v_{\sigma}\sin(\alpha_{h})\cos^{2}(\alpha_{h})-\text{$\lambda_{8}$}v_{\sigma}\sin^{3}(\alpha_{h})+2\text{$\lambda_{8}$}\sin(\alpha_{h})\cos^{2}(\alpha_{h})\,,\end{split}$
(A.8)
$\displaystyle\lambda_{h_{1}\eta_{R}\eta_{R}}^{(0)}=v\cos(\alpha_{h})(\text{$\lambda_{5}$}+\text{$\lambda_{7}$})-\text{$\lambda_{10}$}v_{\sigma}\sin(\alpha_{h})\,,$
(A.9)
$\displaystyle\lambda_{h_{1}\varphi_{R}\varphi_{R}}^{(0)}=\text{$\lambda_{9}$}v\cos(\alpha_{h})-v_{\sigma}\sin(\alpha_{h})(\text{$\lambda_{12}$}+\text{$\lambda_{13}$})\,,$
(A.10)
$\displaystyle\lambda_{h_{1}\sigma_{I}\sigma_{I}}^{(0)}=\text{$\lambda_{8}$}v_{1}\cos(\alpha_{h})-2\text{$\lambda_{4}$}v_{\sigma}\sin(\alpha_{h})\,,$
(A.11)
$\displaystyle\lambda_{h_{1}\eta_{I}\eta_{I}}^{(0)}=v_{1}\cos(\alpha_{h})(\text{$\lambda_{5}$}-\text{$\lambda_{7}$})-\text{$\lambda_{10}$}v_{\sigma}\sin(\alpha_{h})\,,$
(A.12)
$\displaystyle\lambda_{h_{1}\varphi_{I}\varphi_{I}}^{(0)}=\text{$\lambda_{9}$}v_{1}\cos(\alpha_{h})-v_{\sigma}\sin(\alpha_{h})(\text{$\lambda_{12}$}-\text{$\lambda_{13}$})\,,$
(A.13)
$\displaystyle\begin{split}&\lambda_{h_{1}h_{1}h_{2}}^{(0)}=6\text{$\lambda_{1}$}v\sin(\alpha_{h})\cos^{2}(\alpha_{h})+\text{$\lambda_{8}$}v\sin^{3}(\alpha_{h})-2\text{$\lambda_{8}$}v\sin(\alpha_{h})\cos^{2}(\alpha_{h})+\\\
&\hskip
45.00006pt6\text{$\lambda_{4}$}v_{\sigma}\sin^{2}(\alpha_{h})\cos(\alpha_{h})+\text{$\lambda_{8}$}v_{\sigma}\cos^{3}(\alpha_{h})-2\text{$\lambda_{8}$}v_{\sigma}\sin^{2}(\alpha_{h})\cos(\alpha_{h})\,,\end{split}$
(A.14)
$\displaystyle\lambda_{h_{1}h_{1}\eta^{\pm}\eta^{\pm}}^{(0)}=\text{$\lambda_{10}$}\sin^{2}(\alpha_{h})+\text{$\lambda_{5}$}\cos^{2}(\alpha_{h})\,,$
(A.15)
$\displaystyle\lambda_{h_{1}h_{1}h_{1}h_{1}}^{(0)}=6\text{$\lambda_{1}$}\cos^{4}(\alpha_{h})+6\text{$\lambda_{4}$}\sin^{4}(\alpha_{h})+6\text{$\lambda_{8}$}\sin^{2}(\alpha_{h})\cos^{2}(\alpha_{h})\,,$
(A.16)
$\displaystyle\lambda_{h_{1}h_{1}h_{2}h_{2}}^{(0)}=6\text{$\lambda_{1}$}\sin^{2}(\alpha_{h})\cos^{2}(\alpha_{h})+6\text{$\lambda_{4}$}\sin^{2}(\alpha_{h})\cos^{2}(\alpha_{h})+\text{$\lambda_{8}$}\sin^{4}(\alpha_{h})+\text{$\lambda_{8}$}\cos^{4}(\alpha_{h})-$
(A.17) $\displaystyle\hskip
45.00006pt4\text{$\lambda_{8}$}\sin^{2}(\alpha_{h})\cos^{2}(\alpha_{h})\,,$
(A.18)
$\displaystyle\lambda_{h_{1}h_{1}\eta_{R}\eta_{R}}^{(0)}=\text{$\lambda_{10}$}\sin^{2}(\alpha_{h})+\cos^{2}(\alpha_{h})(\text{$\lambda_{5}$}+\text{$\lambda_{7}$})\,,$
(A.19)
$\displaystyle\lambda_{h_{1}h_{1}\varphi_{R}\varphi_{R}}^{(0)}=\sin^{2}(\alpha_{h})(\text{$\lambda_{12}$}+\text{$\lambda_{13}$})+\text{$\lambda_{9}$}\cos^{2}(\alpha_{h})\,,$
(A.20)
$\displaystyle\lambda_{h_{1}h_{1}\sigma_{I}\sigma_{I}}^{(0)}=2\text{$\lambda_{4}$}\sin^{2}(\alpha_{h})+\text{$\lambda_{8}$}\cos^{2}(\alpha_{h})\,,$
(A.21)
$\displaystyle\lambda_{h_{1}h_{1}\eta_{I}\eta_{I}}^{(0)}=\text{$\lambda_{10}$}\sin^{2}(\alpha_{h})+\cos^{2}(\alpha_{h})(\text{$\lambda_{5}$}-\text{$\lambda_{7}$})\,,$
(A.22)
$\displaystyle\lambda_{h_{1}h_{1}\varphi_{I}\varphi_{I}}^{(0)}=\sin^{2}(\alpha_{h})(\text{$\lambda_{12}$}-\text{$\lambda_{13}$})+\text{$\lambda_{9}$}\cos^{2}(\alpha_{h})\,,$
(A.23)

such that the trilinear Higgs coupling at one-loop accuracy is determined as
$\lambda_{hhh}^{\mathrm{NLO}}=\lambda_{h_{1}h_{1}h_{1}}^{(0)}+\lambda_{hhh}^{t}+\lambda_{hhh}^{\phi}$.

## References

  * [1] Y. Cai, J. Herrero-García, M. A. Schmidt, A. Vicente, and R. R. Volkas, From the trees to the forest: a review of radiative neutrino mass models, Front. in Phys. 5 (2017) 63, [1706.08524].
  * [2] P. Minkowski, $\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?, Phys. Lett. B 67 (1977) 421–428.
  * [3] T. Yanagida, Horizontal gauge symmetry and masses of neutrinos, Conf. Proc. C 7902131 (1979) 95–99.
  * [4] S. L. Glashow, The Future of Elementary Particle Physics, NATO Sci. Ser. B 61 (1980) 687.
  * [5] R. N. Mohapatra and G. Senjanovic, Neutrino Mass and Spontaneous Parity Nonconservation, Phys. Rev. Lett. 44 (1980) 912.
  * [6] M. Gell-Mann, P. Ramond, and R. Slansky, Complex Spinors and Unified Theories, Conf. Proc. C 790927 (1979) 315–321, [1306.4669].
  * [7] J. Schechter and J. W. F. Valle, Neutrino Masses in $SU(2)\times U(1)$ Theories, Phys. Rev. D 22 (1980) 2227.
  * [8] J. Schechter and J. W. F. Valle, Neutrino Decay and Spontaneous Violation of Lepton Number, Phys. Rev. D 25 (1982) 774.
  * [9] C. Bonilla, E. Ma, E. Peinado, and J. W. F. Valle, Two-loop Dirac neutrino mass and WIMP dark matter, Phys. Lett. B 762 (2016) 214–218, [1607.03931].
  * [10] S. Baek, H. Okada, and Y. Orikasa, A Two Loop Radiative Neutrino Model, Nucl. Phys. B 941 (2019) 744–754, [1703.00685].
  * [11] S. Saad, Origin of a two-loop neutrino mass from SU(5) grand unification, Phys. Rev. D 99 (2019), no. 11 115016, [1902.11254].
  * [12] T. Nomura and H. Okada, A two loop induced neutrino mass model with modular $A_{4}$ symmetry, Nucl. Phys. B 966 (2021) 115372, [1906.03927].
  * [13] C. Arbeláez, A. E. Cárcamo Hernández, R. Cepedello, M. Hirsch, and S. Kovalenko, Radiative type-I seesaw neutrino masses, Phys. Rev. D 100 (2019), no. 11 115021, [1910.04178].
  * [14] S. Saad, Combined explanations of $(g-2)_{\mu}$, $R_{D^{(*)}}$, $R_{K^{(*)}}$ anomalies in a two-loop radiative neutrino mass model, Phys. Rev. D 102 (2020), no. 1 015019, [2005.04352].
  * [15] Z.-z. Xing and D. Zhang, On the two-loop radiative origin of the smallest neutrino mass and the associated Majorana CP phase, Phys. Lett. B 807 (2020) 135598, [2005.05171].
  * [16] C.-H. Chen and T. Nomura, Two-loop radiative seesaw, muon g $-$ 2, and $\tau$-lepton-flavor violation with DM constraints, JHEP 09 (2021) 090, [2001.07515].
  * [17] T. Nomura, H. Okada, and Y. Uesaka, A two-loop induced neutrino mass model, dark matter, and LFV processes $\ell_{i}\to\ell_{j}\gamma$, and $\mu e\to ee$ in a hidden local $U(1)$ symmetry, Nucl. Phys. B 962 (2021) 115236, [2008.02673].
  * [18] S. Jana, P. K. Vishnu, and S. Saad, Minimal realizations of Dirac neutrino mass from generic one-loop and two-loop topologies at $d=5$, JCAP 04 (2020) 018, [1910.09537].
  * [19] C. Arbeláez, R. Cepedello, J. C. Helo, M. Hirsch, and S. Kovalenko, How many 1-loop neutrino mass models are there?, JHEP 08 (2022) 023, [2205.13063].
  * [20] C. Bonilla, S. Centelles-Chuliá, R. Cepedello, E. Peinado, and R. Srivastava, Dark matter stability and Dirac neutrinos using only Standard Model symmetries, Phys. Rev. D 101 (2020), no. 3 033011, [1812.01599].
  * [21] P. Basler, M. Mühlleitner, and J. Müller, Electroweak Phase Transition in Non-Minimal Higgs Sectors, JHEP 05 (2020) 016, [1912.10477].
  * [22] H. Bahl, J. Braathen, and G. Weiglein, New Constraints on Extended Higgs Sectors from the Trilinear Higgs Coupling, Phys. Rev. Lett. 129 (2022), no. 23 231802, [2202.03453].
  * [23] P. Osland, P. N. Pandita, and L. Selbuz, Trilinear Higgs couplings in the two Higgs doublet model with CP violation, Phys. Rev. D 78 (2008) 015003, [0802.0060].
  * [24] J. E. Camargo-Molina, A. P. Morais, R. Pasechnik, M. O. P. Sampaio, and J. Wessén, All one-loop scalar vertices in the effective potential approach, JHEP 08 (2016) 073, [1606.07069].
  * [25] J. R. Ellis, G. Ridolfi, and F. Zwirner, Radiative corrections to the masses of supersymmetric Higgs bosons, Phys. Lett. B 257 (1991) 83–91.
  * [26] ATLAS Collaboration, Constraining the Higgs boson self-coupling from single- and double-Higgs production with the ATLAS detector using $pp$ collisions at $\sqrt{s}=13$ TeV, tech. rep., CERN, Geneva, 2022.  All figures including auxiliary figures are available at https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/CONFNOTES/ATLAS-CONF-2022-050.
  * [27] ATLAS Collaboration, Combination of searches for non-resonant and resonant Higgs boson pair production in the $b\bar{b}\gamma\gamma$, $b\bar{b}\tau^{+}\tau^{-}$ and $b\bar{b}b\bar{b}$ decay channels using $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector, tech. rep., CERN, Geneva, 2021.  All figures including auxiliary figures are available at https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/CONFNOTES/ATLAS-CONF-2021-052.
  * [28] ATLAS Collaboration, Snowmass White Paper Contribution: Physics with the Phase-2 ATLAS and CMS Detectors, tech. rep., CERN, Geneva, 2022.  All figures including auxiliary figures are available at https://atlas.web.cern.ch/Atlas/GROUPS/PHYSICS/PUBNOTES/ATL-PHYS-PUB-2022-018.
  * [29] S. Homiller and P. Meade, Measurement of the Triple Higgs Coupling at a HE-LHC, JHEP 03 (2019) 055, [1811.02572].
  * [30] Y. Kajiyama, H. Okada, and T. Toma, Multicomponent dark matter particles in a two-loop neutrino model, Phys. Rev. D 88 (2013), no. 1 015029, [1303.7356].
  * [31] E. Ma and M. Raidal, Neutrino mass, muon anomalous magnetic moment, and lepton flavor nonconservation, Phys. Rev. Lett. 87 (2001) 011802, [hep-ph/0102255]. [Erratum: Phys.Rev.Lett. 87, 159901 (2001)].
  * [32] T. Toma and A. Vicente, Lepton Flavor Violation in the Scotogenic Model, JHEP 01 (2014) 160, [1312.2840].
  * [33] A. Vicente and C. E. Yaguna, Probing the scotogenic model with lepton flavor violating processes, JHEP 02 (2015) 144, [1412.2545].
  * [34] M. Lindner, M. Platscher, and F. S. Queiroz, A Call for New Physics : The Muon Anomalous Magnetic Moment and Lepton Flavor Violation, Phys. Rept. 731 (2018) 1–82, [1610.06587].
  * [35] MEG Collaboration, A. M. Baldini et al., Search for the lepton flavour violating decay $\mu^{+}\rightarrow\mathrm{e}^{+}\gamma$ with the full dataset of the MEG experiment, Eur. Phys. J. C 76 (2016), no. 8 434, [1605.05081].
  * [36] R. H. Bernstein and P. S. Cooper, Charged Lepton Flavor Violation: An Experimenter’s Guide, Phys. Rept. 532 (2013) 27–64, [1307.5787].
  * [37] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo, Global constraints on absolute neutrino masses and their ordering, Phys. Rev. D 95 (2017), no. 9 096014, [2003.08511]. [Addendum: Phys.Rev.D 101, 116013 (2020)].
  * [38] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, The fate of hints: updated global analysis of three-flavor neutrino oscillations, JHEP 09 (2020) 178, [2007.14792].
  * [39] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle, 2020 global reassessment of the neutrino oscillation picture, JHEP 02 (2021) 071, [2006.11237].
  * [40] LIGO Scientific, Virgo Collaboration, B. P. Abbott et al., Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116 (2016), no. 6 061102, [1602.03837].
  * [41] A. Greljo, T. Opferkuch, and B. A. Stefanek, Gravitational Imprints of Flavor Hierarchies, Phys. Rev. Lett. 124 (2020), no. 17 171802, [1910.02014].
  * [42] S. F. King, S. Pascoli, J. Turner, and Y.-L. Zhou, Gravitational Waves and Proton Decay: Complementary Windows into Grand Unified Theories, Phys. Rev. Lett. 126 (2021), no. 2 021802, [2005.13549].
  * [43] A. Addazi, A. Marcianò, A. P. Morais, R. Pasechnik, R. Srivastava, and J. W. F. Valle, Gravitational footprints of massive neutrinos and lepton number breaking, Phys. Lett. B 807 (2020) 135577, [1909.09740].
  * [44] A. Kosowsky and M. S. Turner, Gravitational radiation from colliding vacuum bubbles: envelope approximation to many bubble collisions, Phys. Rev. D 47 (1993) 4372–4391, [astro-ph/9211004].
  * [45] K. Hashino, R. Jinno, M. Kakizaki, S. Kanemura, T. Takahashi, and M. Takimoto, Selecting models of first-order phase transitions using the synergy between collider and gravitational-wave experiments, Phys. Rev. D 99 (2019), no. 7 075011, [1809.04994].
  * [46] A. P. Morais and R. Pasechnik, Probing multi-step electroweak phase transition with multi-peaked primordial gravitational waves spectra, JCAP 04 (2020) 036, [1910.00717].
  * [47] J. M. Cline, A. Friedlander, D.-M. He, K. Kainulainen, B. Laurent, and D. Tucker-Smith, Baryogenesis and gravity waves from a UV-completed electroweak phase transition, Phys. Rev. D 103 (2021), no. 12 123529, [2102.12490].
  * [48] F. F. Freitas, G. Lourenço, A. P. Morais, A. Nunes, J. a. Olívia, R. Pasechnik, R. Santos, and J. a. Viana, Impact of SM parameters and of the vacua of the Higgs potential in gravitational waves detection, JCAP 03 (2022), no. 03 046, [2108.12810].
  * [49] Z. Zhang, C. Cai, X.-M. Jiang, Y.-L. Tang, Z.-H. Yu, and H.-H. Zhang, Phase transition gravitational waves from pseudo-Nambu-Goldstone dark matter and two Higgs doublets, JHEP 05 (2021) 160, [2102.01588].
  * [50] P. Di Bari, D. Marfatia, and Y.-L. Zhou, Gravitational waves from first-order phase transitions in Majoron models of neutrino mass, JHEP 10 (2021) 193, [2106.00025].
  * [51] L. Dolan and R. Jackiw, Symmetry Behavior at Finite Temperature, Phys. Rev. D 9 (1974) 3320–3341.
  * [52] S. Weinberg, Gauge and Global Symmetries at High Temperature, Phys. Rev. D 9 (1974) 3357–3378.
  * [53] R. H. Brandenberger, Quantum Field Theory Methods and Inflationary Universe Models, Rev. Mod. Phys. 57 (1985) 1.
  * [54] M. Quiros, Finite temperature field theory and phase transitions, in ICTP Summer School in High-Energy Physics and Cosmology, pp. 187–259, 1, 1999.  hep-ph/9901312.
  * [55] D. Curtin, P. Meade, and H. Ramani, Thermal Resummation and Phase Transitions, Eur. Phys. J. C 78 (2018), no. 9 787, [1612.00466].
  * [56] R. R. Parwani, Resummation in a hot scalar field theory, Phys. Rev. D 45 (1992) 4695, [hep-ph/9204216]. [Erratum: Phys.Rev.D 48, 5965 (1993)].
  * [57] P. B. Arnold and O. Espinosa, The Effective potential and first order phase transitions: Beyond leading-order, Phys. Rev. D 47 (1993) 3546, [hep-ph/9212235]. [Erratum: Phys.Rev.D 50, 6662 (1994)].
  * [58] J. R. Espinosa and M. Quiros, Improved metastability bounds on the standard model Higgs mass, Phys. Lett. B 353 (1995) 257–266, [hep-ph/9504241].
  * [59] J. M. Moreno, M. Quiros, and M. Seco, Bubbles in the supersymmetric standard model, Nucl. Phys. B 526 (1998) 489–500, [hep-ph/9801272].
  * [60] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Shape of the acoustic gravitational wave power spectrum from a first order phase transition, Phys. Rev. D 96 (2017), no. 10 103520, [1704.05871]. [Erratum: Phys.Rev.D 101, 089902 (2020)].
  * [61] M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Numerical simulations of acoustically generated gravitational waves at a first order phase transition, Phys. Rev. D 92 (2015), no. 12 123009, [1504.03291].
  * [62] T. V. I. Tenkanen and J. van de Vis, Speed of sound in cosmological phase transitions and effect on gravitational waves, JHEP 08 (2022) 302, [2206.01130].
  * [63] W.-Y. Ai, B. Laurent, and J. van de Vis, Model-independent bubble wall velocities in local thermal equilibrium, 2303.10171.
  * [64] A. Ekstedt, P. Schicho, and T. V. I. Tenkanen, DRalgo: a package for effective field theory approach for thermal phase transitions, 2205.08815.
  * [65] P. M. Schicho, T. V. I. Tenkanen, and J. Österman, Robust approach to thermal resummation: Standard Model meets a singlet, JHEP 06 (2021) 130, [2102.11145].
  * [66] L. Niemi, P. Schicho, and T. V. I. Tenkanen, Singlet-assisted electroweak phase transition at two loops, Phys. Rev. D 103 (2021), no. 11 115035, [2103.07467].
  * [67] P. Schicho, T. V. I. Tenkanen, and G. White, Combining thermal resummation and gauge invariance for electroweak phase transition, JHEP 11 (2022) 047, [2203.04284].
  * [68] D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, and G. White, Theoretical uncertainties for cosmological first-order phase transitions, JHEP 04 (2021) 055, [2009.10080].
  * [69] D. Bodeker and G. D. Moore, Electroweak Bubble Wall Speed Limit, JCAP 05 (2017) 025, [1703.08215].
  * [70] M. Kierkla, A. Karam, and B. Swiezewska, Conformal model for gravitational waves and dark matter: a status update, JHEP 03 (2023) 007, [2210.07075].
  * [71] T. Prokopec, J. Rezacek, and B. Świeżewska, Gravitational waves from conformal symmetry breaking, JCAP 02 (2019) 009, [1809.11129].
  * [72] J. Ellis, M. Lewicki, and V. Vaskonen, Updated predictions for gravitational waves produced in a strongly supercooled phase transition, JCAP 11 (2020) 020, [2007.15586].
  * [73] LISA Cosmology Working Group Collaboration, P. Auclair et al., Cosmology with the Laser Interferometer Space Antenna, 2204.05434.
  * [74] C. Caprini et al., Detecting gravitational waves from cosmological phase transitions with LISA: an update, JCAP 03 (2020) 024, [1910.13125].
  * [75] C. L. Wainwright, CosmoTransitions: Computing Cosmological Phase Transition Temperatures and Bubble Profiles with Multiple Fields, Comput. Phys. Commun. 183 (2012) 2006–2013, [1109.4189].
  * [76] P. Amaro-Seoane et al., Laser Interferometer Space Antenna, 1702.00786.
  * [77] C. Caprini et al., Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions, JCAP 04 (2016) 001, [1512.06239].
  * [78] V. Corbin and N. J. Cornish, Detecting the cosmic gravitational wave background with the big bang observer, Class. Quant. Grav. 23 (2006) 2435–2446, [gr-qc/0512039].
  * [79] H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Detecting a gravitational-wave background with next-generation space interferometers, Phys. Rev. D 73 (2006) 064006, [gr-qc/0511145].
  * [80] D. J. Weir, PTPlot: a tool for exploring the gravitational wave power spectrum from first-order phase transitions, Aug., 2022.
  * [81] R. Zhou, W. Cheng, X. Deng, L. Bian, and Y. Wu, Electroweak phase transition and Higgs phenomenology in the Georgi-Machacek model, JHEP 01 (2019) 216, [1812.06217].
  * [82] A. Ahriche, K. Hashino, S. Kanemura, and S. Nasri, Gravitational Waves from Phase Transitions in Models with Charged Singlets, Phys. Lett. B 789 (2019) 119–126, [1809.09883].
  * [83] J. Bernon, L. Bian, and Y. Jiang, A new insight into the phase transition in the early Universe with two Higgs doublets, JHEP 05 (2018) 151, [1712.08430].
  * [84] T. Biekötter, S. Heinemeyer, J. M. No, M. O. Olea-Romacho, and G. Weiglein, The trap in the early Universe: impact on the interplay between gravitational waves and LHC physics in the 2HDM, JCAP 03 (2023) 031, [2208.14466].