# Phenomenology of a flavoured multiscalar BGL-like model with three
generations of massive neutrinos

P.M. Ferreira1,2<EMAIL_ADDRESS>Felipe F. Freitas3
<EMAIL_ADDRESS>João Gonçalves3<EMAIL_ADDRESS>Antonio P. Morais3
<EMAIL_ADDRESS>Roman Pasechnik4<EMAIL_ADDRESS>Vasileios
Vatellis3<EMAIL_ADDRESS>
1Instituto Superior de Engenharia de Lisboa — ISEL, 1959-007 Lisboa, Portugal
2Centro de Física Teórica e Computacional, Faculdade de Ciências, Universidade
de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal
3Departamento de Física, Universidade de Aveiro and CIDMA, Campus de Santiago,
3810-183 Aveiro, Portugal
4Department of Astronomy and Theoretical Physics, Lund University, SE-223 62
Lund, Sweden

###### Abstract

In this paper, we present several possible anomaly-free implementations of the
Branco-Grimus-Lavoura (BGL) model with two Higgs doublets and one singlet
scalar. The model also includes three generations of massive neutrinos that
get their mass via a type-I seesaw mechanism. A particular anomaly-free
realization, which we dub $\nu$BGL-1 scenario, is subjected to an extensive
phenomenological analysis, from the perspective of flavour physics and
collider phenomenology.

##  I Introduction

With the discovery of the Higgs boson at the LHC in 2012 Aad  _et al._ (2012);
Chatrchyan  _et al._ (2012), the entire particle spectrum of the Standard
Model (SM) has been observed and experimentally validated. However, despite
its remarkable success, there is experimental evidence for the existence of
Dark Matter (DM) and neutrino masses, while the observed hierarchical
structure that characterizes the fermion sector finds no explanation in the
SM. With a plethora of New Physics (NP) extensions that aim at offering
solutions to the aforementioned phenomena, the vast majority of them come in
the form of extended scalar sectors. One of the most popular and thoroughly
explored multi-scalar scenarios is the Two-Higgs doublet model (THDM), as
first proposed by T.D. Lee Lee (1973). The introduction of a new complex
scalar doublet, in addition to the one that is already present in the SM,
results in an extended particle spectrum with a new electrically charged
scalar and two new neutral ones, one of which is even while the other is odd
under discrete CP (charge times parity) transformations.

This type of models provide a rich breeding ground for phenomenological
studies, with the possibility of having a suitable DM candidate in its
spectrum Cárcamo Hernández  _et al._ (2021); Arcadi  _et al._ (2022); Bell
_et al._ (2017); Dorsch  _et al._ (2017); Camargo  _et al._ (2019); Ivanov
(2017), spontaneous CP-violation Branco  _et al._ (2012); Grzadkowski  _et
al._ (2014); Boto  _et al._ (2020); Grzadkowski  _et al._ (2016); Sokolowska
_et al._ (2008) and generation of flavour-changing neutral currents (FCNCs)
mediated by the new neutral scalars Crivellin  _et al._ (2013); Kim  _et al._
(2015). However, the presence of FCNCs can be troublesome since they are very
constrained experimentally, and in the models we are discussing they can
manifest themselves already at tree-level, whereas equivalent processes in the
SM are of pure radiative origin, thus highly suppressed. As such, one either
requires FCNCs couplings to be fine-tuned or rather there should exist an
underlying symmetry that naturally mitigates the effect of these interactions.
A concrete mechanism for symmetry-suppressed FCNCs has been proposed by
Branco, Grimus and Lavoura Branco  _et al._ (1996) where, due to the presence
of an abelian flavour symmetry such as $\mathbb{Z}_{2}$ or $\mathrm{U(1)}$,
FCNC couplings are kept under control by making such interactions proportional
to off-diagonal elements of the Cabibbo–Kobayashi–Maskawa (CKM) matrix. Such a
class of models is commonly dubbed as BGL models.

A natural extension of the THDM-BGL framework can emerge in the form of a
gauge $\mathrm{U(1)}$ symmetry, which necessarily requires a further
enlargement of the scalar sector with a new EW singlet, thus falling in the
family of Next-to-minimal THDM (or NTHDM) for short. Recently, a complete list
of all possible anomaly-free implementations of the gauged NTHDM with two
families of right-handed neutrinos and a type-I seesaw mechanism was
introduced in Ref. Ordell  _et al._ (2020), while a generalization to three
neutrino generations was discussed in Ref. Nottensteiner  _et al._ (2019). For
previous phenomenological considerations in a set of THDMs featuring the
absence of tree-level FCNCs and accommodating a seesaw mechanism for neutrino
mass generation, see Ref. Campos  _et al._ (2017). Among all allowed textures
in the Yukawa sector one of them replicates the BGL structure. The model
explored in this article is inspired by such a scenario but for the case of a
global $\mathrm{U(1)^{\prime}}$ flavour symmetry. The latter is explicitly
broken via soft terms in the scalar potential as well as by the vacuum
expectation value (VEV) in the real component of the complex singlet scalar
field. As a result, in addition to a SM-like Higgs boson candidate, five new
physical scalars emerge, including an electrically charged one, $H^{\pm}$, two
CP-even neutral ones, $H_{2}$ and $H_{3}$, and two pseudoscalars, $A_{2}$ and
$A_{3}$. Note that in the context of the gauge $\mathrm{U(1)^{\prime}}$
flavour symmetry Ordell  _et al._ (2020); Nottensteiner  _et al._ (2019) one
of those pseudoscalars, as a pure Goldstone particle, becomes a longitudinal
polarization mode of a $Z^{\prime}$ vector boson. In turn, tight constraints
imposed by direct searches at the LHC Aad  _et al._ (2020a); Aaboud  _et al._
(2018a); Aad  _et al._ (2021a); Sirunyan  _et al._ (2020a, 2021, 2018a) on new
gauge bosons imply a rather large $\mathrm{U(1)}^{\prime}$ breaking scale, of
at least a few TeV, naturally inducing masses of the same order to the new
scalars, unless a certain degree of fine-tuning is in place. Our current goal
is to study the impact of the extended scalar sectors not too far from the EW
scale. For this purpose, a global version of the BGL scenario instead of the
gauged one proposed in Ref. Ordell  _et al._ (2020); Nottensteiner  _et al._
(2019) is considered in this work.

We test the viability of the scalar sector against phenomenological
observables, in particular, by making sure that FCNCs obey current
experimental constraints and that there exists one CP-even neutral scalar that
replicates the properties of the SM Higgs boson, defined as $H_{1}$ in this
article. For a rigorous analysis we employ widely used open-source software,
namely, SARAH Staub (2014), for model implementation, and SPheno Porod and
Staub (2012), to numerically calculate masses and interaction vertices, which
are then used to interface with HiggsBounds and HiggsSignals Bechtle  _et al._
(2015, 2014), where the viability of the Higgs sector is verified, and flavio
Straub (2018), where flavour observables are calculated.

The paper is organised as follows. In Sec. II we introduce the model, with the
focus on the Yukawa and scalar sectors of the theory that obey the BGL
structure. In Sec. III we detail the method by which one can obtain anomaly-
free models, as first discussed in Ref. Nottensteiner  _et al._ (2019),
highlighting the possible textures in the lepton sector. In Sec. III we show
all possible anomaly-free implementations of the considered class of THDSM
scenarios. In Sec. IV we discuss the particular version of the model
considered in this article, with a discussion on the neutrino masses and
mixings. In Sec. V we describe the methodology behind the phenomenological
studies and present the numerical results that follow. In Sec. VI we confront
our results with direct scalar searches at the LHC. In Sec. VII we conclude
with a summary of our findings.

##  II A generic NTHDM with a BGL structure

In this work, we extend the SM with a flavour non-universal
$\mathrm{U}(1)_{\mathrm{}}^{\prime}$ global symmetry alongside three
generations of right-handed neutrinos $\nu_{R}^{1,2,3}$, a scalar singlet $S$
and a second Higgs doublet111The BGL model versions discussed in this article
were first introduced in the Appendix of Ordell (2020). $\Phi_{2}$. In
addition, we demand that the singlet and at least one of the doublets
transform non-trivially under $\mathrm{U(1)}^{\prime}$ such that at least two
VEVs contribute to the breaking of the flavor symmetry.

The Yukawa interactions that define the model to be discussed read as

$\displaystyle\begin{split}-\mathcal{L}_{\mathrm{Yukawa}}&=\overline{q_{L}^{0}}\Gamma_{a}\Phi^{a}d_{R}^{0}+\overline{q_{L}^{0}}\Delta_{a}\tilde{\Phi}^{a}u_{R}^{0}+\overline{\ell_{L}^{0}}\Pi_{a}\Phi^{a}e_{R}^{0}+\overline{\ell_{L}^{0}}\Sigma_{a}\tilde{\Phi}^{a}\nu_{R}^{0}\\\
&+\frac{1}{2}\overline{\nu_{R}^{c\,0}}\left(\mathrm{A}+\mathrm{B}S+\mathrm{C}S^{\ast}\right)\nu_{R}^{0}+\mathrm{h.c.}\,,\end{split}$
(1)

where $\Gamma$, $\Delta$, $\Pi$ and $\Sigma$ are the $3\times 3$ Dirac-like
down quark, up quark, charged lepton and neutrino Yukawa coupling matrices,
respectively. The index $a$ runs over the two scalar doublets (with implicit
summation) and, as usual, $\tilde{\Phi}\equiv i\sigma_{2}{\Phi}^{*}$. While B
and C are the Majorana-like Yukawa matrices involving right-handed neutrinos
and the scalar singlet, A is a Majorana mass term written in the flavour
basis, defined by the fields expressed with the $0$ superscript. To comply
with a BGL structure the transformation laws of the different fields under the
$\mathrm{U}(1)_{\mathrm{}}^{\prime}$ global symmetry are chosen so that the
quark Yukawa textures are given as

$\displaystyle\begin{split}\Gamma_{1}&:\begin{pmatrix}\times&\times&\times\\\
\times&\times&\times\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\hskip
5.69054pt\Gamma_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\times&\times&\times\end{pmatrix},\hskip 5.69054pt\\\
\Delta_{1}&:\begin{pmatrix}\times&\times&\makebox[0.0pt][c]{0}\\\
\times&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\hskip
4.2679pt\Delta_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix}.\end{split}$
(2)

Note that this choice of textures implies that tree-level FCNCs will appear
only in the down quark sector, since the up-quark Yukawa matrices can be
diagonalized simultaneously. Indeed, this represents the more conservative
scenario, since down-type FCNCs are more tightly constrained by experiment
when compared to those in the up-sector. Of course, one can also extend the
BGL mechanism to the lepton sector, as first proposed in Ref. Botella  _et
al._ (2011), which is beyond the scope of the present work.

In Eq. (1), the scalar fields can be expanded as

$\displaystyle\begin{split}\Phi_{a}&\equiv\frac{1}{\sqrt{2}}\begin{pmatrix}\sqrt{2}\phi^{+}_{a}\\\
v_{a}\,e^{i\varphi_{a}}+R_{a}+iI_{a}\end{pmatrix}\,,\\\\[2.27621pt]
S&\equiv\frac{1}{\sqrt{2}}\left(v_{S}\,e^{i\varphi_{S}}+\rho+i\eta\right)\,,\end{split}$
(3)

where $v_{1,2}$ and $v_{S}$ denote the doublet and singlet VEVs, respectively.
As the scalar fields have, in general, different $\mathrm{U(1)_{Y}\times
U(1)^{\prime}}$ charge assignment, we can choose two of the VEVs real, for
example $v_{1}$ and $v_{2}$, i.e. $\varphi_{1,2}=0$. In multi-Higgs SM
extensions there is a possibility to trigger spontaneous CP breaking via
complex VEVs. For the purposes of this work, however, we are only focused on a
CP conserving scalar phase. CP violation will emerge via the CKM matrix as in
the SM, due to the fact that the entries of the Yukawa matrices of equation
(2) are complex.

Once the two doublets acquire VEVs, both the quarks and the charged leptons
obtain tree-level masses which are given by

$\displaystyle
M_{u}^{0}\equiv\frac{1}{\sqrt{2}}\left(v_{1}\Delta_{1}+v_{2}\Delta_{2}\right),$
(4) $\displaystyle
M_{d}^{0}=\frac{1}{\sqrt{2}}\left(v_{1}\Gamma_{1}+v_{2}\Gamma_{2}\right),$
$\displaystyle
M_{e}^{0}=\frac{1}{\sqrt{2}}\left(v_{1}\Pi_{1}+v_{2}\Pi_{2}\right).$

One can rotate the spectrum to the mass basis via bi-unitary transformations
of the form

$D_{f}=U_{f\mathrm{L}}^{\dagger}M_{f}^{0}\,U_{f\mathrm{R}},$ (5)

with $D_{f}$ a diagonal mass form with ordered fermion masses and where
$U_{\mathrm{f\,L,R}}$ are unitary matrices with the subscript $f$ indicating
the fermion flavour, i.e. $f=e,u,d$. Following the standard notation, one can
define the following matrices

$\displaystyle
N_{u}^{0}=\frac{1}{\sqrt{2}}\big{(}v_{2}{\Delta_{1}}-v_{1}{\Delta_{2}}\big{)}\,\quad\quad
N_{d}^{0}=\frac{1}{\sqrt{2}}\big{(}v_{2}{\Gamma_{1}}-v_{1}{\Gamma_{2}}\big{)}\,,$
(6)

whose off-diagonal elements are responsible for inducing tree-level FCNC
interactions. The key feature of the BGL model is that such matrices can be
re-expressed solely in terms of quark masses, CKM mixing elements
$V_{ij}\equiv\sum_{k}(U_{u\mathrm{L}})_{ik}(U_{d\mathrm{L}}^{\dagger})_{kj}$,
with $i,j,k=1,2,3$, and the ratio $t_{\beta}\equiv\tan\beta\equiv v_{1}/v_{2}$
as follows

$\displaystyle\begin{split}\left(N_{u}\right)_{ij}&=\left(t_{\beta}\delta_{ij}-\left(t_{\beta}+t_{\beta}^{-1}\right)\delta_{ij}\delta_{j3}\right)m_{u_{j}},\\\
\left(N_{d}\right)_{ij}&=\left(t_{\beta}\delta_{ij}-\left(t_{\beta}+t_{\beta}^{-1}\right)V_{3i}^{*}V_{3j}\right)m_{d_{j}}\,,\end{split}$
(7)

such that the down-specific FCNCs are controlled by the smallness of the off-
diagonal CKM mixing elements $V_{3i}$.

For the charged lepton sector, the generic coupling combinations can be found
in complete analogy with the quarks, i.e.

$\displaystyle
M_{e}^{0}=\frac{1}{\sqrt{2}}\left(v_{1}\Pi_{1}+v_{2}\Pi_{2}\right)\,,\quad\quad
N_{e}^{0}=\frac{1}{\sqrt{2}}\left(v_{2}\Pi_{1}-v_{1}\Pi_{2}\right)\,.$ (8)

Last but not least, the neutrino sector is equipped with a seesaw mechanism of
type-I, as soon as the singlet $S$ develops a VEV. In particular, defining

$\displaystyle n_{L}^{0}\equiv\begin{pmatrix}\nu_{L}^{0}\\\
\nu_{R}^{c}\end{pmatrix},$ (9)

one can recast the neutrino mass Lagrangian as

$\displaystyle\begin{split}-\mathcal{L}_{\nu}^{\mathrm{mass}}=&\frac{1}{2}\overline{n_{L}^{0}}\mathcal{M}n_{L}^{0,c}+\text{h.c.}\,,\end{split}$
(10)

where

$\displaystyle\mathcal{M}\equiv\begin{pmatrix}\mathbb{0}&m_{D}\\\
m_{D}^{T}&M_{R}\end{pmatrix},$ (11)

and with the Dirac and Majorana neutrino masses given by

$\displaystyle\begin{split}m_{D}&\equiv\frac{1}{\sqrt{2}}\left(v_{1}\Sigma_{1}+v_{2}\Sigma_{2}\right),\\\
M_{R}&\equiv\mathrm{A}+\frac{v_{S}}{\sqrt{2}}\left(\mathrm{B}+\mathrm{C}\right)\,,\end{split}$
(12)

respectively. Here, we have also considered that the singlet VEV is real. At
low energies, in the limit of $M_{R}\gg m_{D}/m_{D}^{T}$, the effective
Lagrangian for the light active neutrinos takes the form

$\displaystyle\begin{split}-\mathcal{L}_{\nu}^{\mathrm{eff}}\equiv\frac{1}{2}\overline{\nu_{L}^{0}}m_{\nu}\nu_{L}^{0,c}+\text{h.c.},\end{split}$
(13)

with

$\displaystyle m_{\nu}\equiv-m_{D}M_{R}^{-1}m_{D}^{T}\,.$ (14)

The numerical analysis performed in this study relies on an inversion
procedure analogous to that developed in Das  _et al._ (2021). In particular,
for the quark sector, one uses their physical masses as well as the CKM mixing
elements as input parameters. However, due to a rich neutrino content, one
cannot analytically invert the entire lepton sector. Instead, while electron,
muon and tau masses are given as input parameters, one performs a numerical
fit to neutrino masses alongside an extended version of the Pontecorvo-Maki-
Nakagawa-Sakata (PMNS) matrix, as it is discussed in Sec. IV.

The most generic form of the scalar potential can be written as

$\displaystyle V=V_{0}+V_{1}$ (15)

with

$\displaystyle\begin{split}V_{0}=\quad&\mu_{i}^{2}|\Phi^{i}|^{2}+\lambda_{i}|\Phi^{i}|^{4}+\lambda_{3}|\Phi_{1}|^{2}|\Phi_{2}|^{2}+\lambda_{4}|\Phi_{1}^{\dagger}\Phi_{2}|^{2}+{\mu_{S}}^{2}|S|^{2}+\lambda_{1}^{\prime}|S|^{4}\\\
&+\lambda_{2}^{\prime}|\Phi_{1}|^{2}|S|^{2}+\lambda_{3}^{\prime}|\Phi_{2}|^{2}|S|^{2}\qquad\textrm{and}\\\
V_{1}=\quad&\mu_{3}^{2}\Phi_{2}^{\dagger}\Phi_{1}+\frac{1}{2}\mu_{b}^{2}S^{2}+a_{1}\Phi_{1}^{\dagger}\Phi_{2}S+a_{2}\Phi_{1}^{\dagger}\Phi_{2}S^{\dagger}+a_{3}\Phi_{1}^{\dagger}\Phi_{2}S^{2}+a_{4}\Phi_{1}^{\dagger}\Phi_{2}{S^{\dagger}}^{2}+\mathrm{h.c.}\,.\end{split}$
(16)

While in $V_{0}$ we adopt the usual notation where $\mu_{i,S}$ represent
quadratic mass parameters and $\lambda_{i}$, $\lambda_{i}^{\prime}$ are the
quartic couplings, in $V_{1}$ $a_{1,2}$ denote cubic couplings between the
Higgs doublets and the singlet while $a_{3,4}$ are the quartic couplings
between the two Higgs doublets and the scalar singlet field. Given that the
singlet $S$ carries a non-trivial $\mathrm{U}(1)_{\mathrm{}}^{\prime}$ charge
$X_{S}$, then, out of the four $a_{1,2,3,4}$ terms, only one is allowed in the
limit of an exact $\mathrm{U}(1)_{\mathrm{}}^{\prime}$. However, both $a_{1}$
and $a_{2}$, as well as $\mu_{b}^{2}$, can be introduced to softly break the
flavour symmetry and are allowed to coexist with either $a_{3}$ or $a_{4}$
Nottensteiner  _et al._ (2019); Ordell  _et al._ (2020).

##  III Anomaly-free BGL implementation

The transformation properties of the quark and scalar fields under the
$\mathrm{U}(1)_{\mathrm{}}^{\prime}$ flavor symmetry are given, as usual, by

$\Psi_{i}\rightarrow e^{i\theta X_{\Psi_{i}}}\Psi_{i}\,,$ (17)

where $\Psi_{i}\subset(q_{L}^{0},u_{R}^{0},d_{R}^{0},\Phi_{1},\Phi_{2},S)$,
and where $X_{\Psi_{i}}$ are the corresponding $\mathrm{U(1)}^{\prime}$
charges with $\theta$ a global phase parameter. It then follows that the
Yukawa matrices transform as

$\displaystyle\begin{split}(\Gamma_{a})_{ij}&=e^{i\theta(X_{q_{i}}-X_{d_{j}}-X_{\Phi_{a}})}(\Gamma_{a})_{ij}\,,\\\
(\Delta_{a})_{ij}&=e^{i\theta(X_{q_{i}}-X_{u_{j}}+X_{\Phi_{a}})}(\Delta_{a})_{ij}\,,\end{split}$
(18)

such that textures on $(\Gamma_{a})_{ij}$ result from a set of linear
constraints that read as

$\displaystyle\begin{split}(\Gamma_{a})_{ij}&=\mathrm{any\;\;if}\;\;X_{q_{i}}-X_{d_{j}}=X_{\Phi_{a}}\,,\\\
(\Gamma_{a})_{ij}&=0\;\;\;\;\;\;\mathrm{if}\;\;X_{q_{i}}-X_{d_{j}}\neq
X_{\Phi_{a}}\,.\end{split}$ (19)

The same can be done in the up-quark sector by trading
$(\Gamma_{a})_{ij}\to(\Delta_{a})_{ij}$, $X_{d_{j}}\to X_{u_{j}}$ and
$X_{\Phi_{a}}\to-X_{\Phi_{a}}$. Among many possible solutions found in Refs.
Nottensteiner  _et al._ (2019) and Ordell  _et al._ (2020), the well-known
quark BGL textures (see Eq. (2)), that are used in this work, satisfy the
following 36 constraints

$\displaystyle\begin{split}&X_{q_{1,2}}-X_{d_{1,2,3}}=X_{\Phi_{1}},\;\;X_{q_{3}}-X_{d_{1,2,3}}\neq
X_{\Phi_{1}}\,,\\\
&X_{q_{3}}-X_{d_{1,2,3}}=X_{\Phi_{2}},\;\;X_{q_{1,2}}-X_{d_{1,2,3}}\neq
X_{\Phi_{2}}\,,\\\
&X_{q_{1,2}}-X_{u_{1,2}}=-X_{\Phi_{1}},\;\;X_{q_{3}}-X_{u_{1,2,3}}\neq-
X_{\Phi_{1}}\,,\\\ &X_{q_{1,2}}-X_{u_{3}}\neq-
X_{\Phi_{1}},\;\;X_{q_{3}}-X_{u_{3}}=-X_{\Phi_{2}}\,,\\\
&X_{q_{1,2}}-X_{u_{1,2,3}}\neq-X_{\Phi_{2}},\;\;X_{q_{3}}-X_{u_{1,2}}\neq-
X_{\Phi_{2}}\,.\end{split}$ (20)

Focusing now on the lepton and neutrino sectors, one must consider all
possible textures of the Dirac-like $\Pi_{a}$ and $\Sigma_{a}$ and Majorana-
like $B$ and $C$ Yukawa matrices, as well as on the $A$ mass matrix. To do so
one considers that

  * (i)

there are no massless charged leptons yielding $\text{det}\hskip
1.42262ptM_{e}\neq 0$;

  * (ii)

there are three generations of massive neutrinos such that222Note that there
exists no anomaly-free implementation of the BGL textures with two generations
of massive neutrinos, as shown in Ref. Ordell  _et al._ (2020).
$\text{det}\hskip 1.42262ptM_{\nu}\neq 0$;

  * (iii)

there is a non-zero complex phase in the PMNS matrix which implies that
$\text{det}\hskip 1.42262pt[M_{e}M_{e}^{\dagger}]\neq 0$ and $\text{det}\hskip
1.42262pt[M_{\nu}M_{\nu}^{\dagger}]\neq 0$ .

Notice that the second condition is the only one differing from those first
introduced in Ref. Nottensteiner  _et al._ (2019) and Ordell  _et al._ (2020).
For a type-I seesaw mechanism such a condition also translates into
$M_{\mathrm{R}}$ and $M_{\mathrm{D}}$ being $3\times 3$ matrices with non-zero
determinant. Note also that any two models that can be mapped from one another
via a permutation are equivalent as it would simply correspond to a
relabelling of flavor indices. Defining the new primed Yukawa and mass
matrices as

$\displaystyle\begin{split}&\Gamma^{\prime}_{1,2}=\mathcal{P}_{i}^{\mathrm{T}}\Gamma_{1,2}\mathcal{P}_{j}\,,\;\;\Delta^{\prime}_{1,2}=\mathcal{P}_{i}^{\mathrm{T}}\Delta_{1,2}\mathcal{P}_{k}\,,\\\
&\Pi^{\prime}_{1,2}=\mathcal{P}_{l}^{\mathrm{T}}\Pi_{1,2}\mathcal{P}_{m}\,,\;\;\Sigma^{\prime}_{1,2}=\mathcal{P}_{l}^{\mathrm{T}}\Sigma_{1,2}{\mathcal{P}}_{n}\,,\\\
&A^{\prime}={\mathcal{P}}_{n}^{\mathrm{T}}A\hskip
1.13809pt{\mathcal{P}}_{n}\,,\;\;B^{\prime}={\mathcal{P}}_{n}^{\mathrm{T}}B\hskip
1.13809pt{\mathcal{P}}_{n}\,,\;\;C^{\prime}={\mathcal{P}}_{n}^{\mathrm{T}}C\hskip
1.13809pt{\mathcal{P}}_{n}\,,\end{split}$ (21)

with $\mathcal{P}$ being a three-dimensional representation of the permutation
group $S_{3}$, with all indices running from one to six, one excludes all such
additional textures that can be obtained by permutations in the flavour space.

### Majorana neutrino sector

In total, there are 11 minimal textures for $A$, $B$ and $C$ that fulfil the
constraint of $M_{\mathrm{R}}$ being a $3\times 3$ symmetric matrix with a
non-zero determinant. These are

$\displaystyle\begin{split}(1)&\;\;\;\;A:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;B:\mathbb{0},\;\;C:\mathbb{0},\\\
(2)&\;\;\;\;A:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;B:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;C:\mathbb{0},\\\\[2.84526pt]
(3)&\;\;\;\;A:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;B:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;C:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix},\\\\[2.84526pt]
(4)&\;\;\;\;A:\begin{pmatrix}\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;B:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;C:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\end{split}$
(22)

where textures (1) and (2) come in three and six versions, respectively,
taking into account all possible permutations of $A$, $B$ and $C$. For
textures (3) and (4), on the other hand, we only need to consider those shown
above as permutations of rows and columns solely correspond to a relabelling
of flavor indices. In the presence of the $U(1)^{\prime}$ flavour symmetry one
must also consider the transformation laws

$\displaystyle A_{ij}$
$\displaystyle=e^{i\alpha(X_{\nu_{i}}+X_{\nu_{j}})}A_{ij},$ (23)
$\displaystyle B_{ij}$
$\displaystyle=e^{i\alpha(X_{\nu_{i}}+X_{\nu_{j}}+X_{S})}B_{ij},$
$\displaystyle C_{ij}$
$\displaystyle=e^{i\alpha(X_{\nu_{i}}+X_{\nu_{j}}-X_{S})}C_{ij}.$

Therefore, texture (1) implies that

$2X_{\nu_{1}}=0,\;\;X_{\nu_{2}}+X_{\nu_{3}}=0,$ (24)

while texture (4) results in

$\displaystyle X_{\nu_{1}}+X_{\nu_{2}}$
$\displaystyle=0,\;\;X_{\nu_{1}}+X_{\nu_{3}}=X_{S},$ (25) $\displaystyle
X_{\nu_{2}}+X_{\nu_{3}}$ $\displaystyle=-X_{S},$

and so on. In addition to this, one of the following four conditions (see Ref.
Ordell  _et al._ (2020))

$\displaystyle
X_{S}=\pm\left(X_{\Phi_{1}}-X_{\Phi_{2}}\right),\;\;X_{S}=\pm\frac{1}{2}\left(X_{\Phi_{1}}-X_{\Phi_{2}}\right),$
(26)

can be extracted from the scalar potential $V_{1}$ depending on which of the
$a_{1,2,3,4}$ terms is invariant under the flavour
$\mathrm{U}(1)_{\mathrm{}}^{\prime}$ symmetry. However, if only $A$ has a non-
zero texture then $S$ is neutral under the flavour symmetry and the four
conditions above are all automatically valid. In that case the cubic $a_{1,2}$
and quartic $a_{3,4}$ can simultaneously coexist.

### Charged lepton sector

For the charged lepton sector, the same minimal textures already discussed in
Ref. Nottensteiner  _et al._ (2019) and Ordell  _et al._ (2020) apply, i.e.

$\displaystyle\begin{split}(1)\;\;\;\;\;\;\Pi_{1}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix},\;\;\Pi_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\\\
(2)\;\;\;\;\;\;\Pi_{1}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Pi_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix},\\\
(3)\;\;\;\;\;\;\Pi_{1}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Pi_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix},\\\
(4)\;\;\;\;\;\;\Pi_{1}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Pi_{2}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix}.\end{split}$
(27)

### Dirac neutrino sector

For the Dirac neutrino textures, permutation of rows and columns are no longer
independent from those in the charged lepton and Majorana neutrino sectors. As
a result, besides fulfilling the constraint of $M_{\mathrm{D}}$ having a non-
zero determinant, we must also consider textures that are equivalent up to
permutations. In total, this amounts to six possible minimal combined
textures,

$\displaystyle\begin{split}1:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix},\;\;2:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;3:\begin{pmatrix}\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix},\\\
4:\begin{pmatrix}\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;5:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;6:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\end{split}$

which in turn corresponds to 48 possible textures for $\Sigma_{1,2}$ – eight
for each of the textures displayed above; 111, 112, 121, 211, 122, 212, 221
and 222, where the numbers indicate the corresponding non-zero entry in
$\Sigma_{1}$ or $\Sigma_{2}$. As an example, the eight possibilities for
texture number 2 are given by

$\displaystyle\begin{split}(111)\;\;\;\;\;\;\Sigma_{1}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Sigma_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\\\
(112)\;\;\;\;\;\;\Sigma_{1}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Sigma_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\\\
(121)\;\;\;\;\;\;\Sigma_{1}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Sigma_{2}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\\\
(211)\;\;\;\;\;\;\Sigma_{1}:\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix},\;\;\Sigma_{2}:\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix},\end{split}$

and so on for 122, 212, 221, 222.

### Anomaly cancellation conditions

While global anomalies are typically regarded as benign, one must recall that
the work previously developed in Nottensteiner  _et al._ (2019); Ordell
(2020), where this article is inspired on, considers a local
$\mathrm{U}(1)_{\mathrm{}}^{\prime}$ symmetry where gauge anomalies must
necessarily be forbidden. With this in mind, and with the purpose of making
the considered model consistent with a gauged version (to be studied
elsewhere), one must also include a set of restrictions that incorporate the
$\mathrm{U}(1)_{\mathrm{}}^{\prime}$ charges in the following triangle
anomalies

$\displaystyle\begin{split}\left[\mathrm{U(1)^{\prime}}\right]^{3},\;\;\mathrm{U(1)^{\prime}}\left[\mathrm{Gravity}\right]^{2},\\\
\mathrm{U(1)^{\prime}}\left[\mathrm{U(1)}_{\mathrm{Y}}\right]^{2},\;\;\mathrm{U(1)^{\prime}}\left[\mathrm{SU(2)_{L}}\right]^{2},\\\
\mathrm{U(1)^{\prime}}\left[\mathrm{SU(3)_{C}}\right]^{2},\;\;\left[\mathrm{U(1)^{\prime}}\right]^{2}\mathrm{U(1)_{Y}}\,.\end{split}$
(28)

As the right-handed neutrinos are solely charged under
$\mathrm{U(1)^{\prime}}$, confer Tab. 2, the only anomalies altered in
comparison to Nottensteiner  _et al._ (2019) are those coming from
$\left[\mathrm{U(1)^{\prime}}\right]^{3}$ and
$\mathrm{U(1)^{\prime}}\left[\mathrm{Gravity}\right]^{2}$, i.e.

$\displaystyle\begin{split}A_{{}_{\mathrm{U(1)^{\prime}}\mathrm{U(1)^{\prime}}\mathrm{U(1)^{\prime}}}}&\equiv\sum_{i=1}^{3}\left(6X_{q_{i}}^{3}+2X_{l_{i}}^{3}-3X_{u_{i}}^{3}-3X_{d_{i}}^{3}-X_{e_{i}}^{3}-X_{\nu_{i}}^{3}\right)=0,\\\
A_{{}_{\mathrm{g}\mathrm{g}\mathrm{U(1)^{\prime}}}}&\equiv\sum_{i=1}^{3}\left(6X_{q_{i}}+2X_{l_{i}}-3X_{u_{i}}-3X_{d_{i}}-X_{e_{i}}-X_{\nu_{i}}\right)=0\,,\end{split}$
(29)

where $g$ denotes gravitational interactions while
$A_{{}_{X,Y,Z}}\propto\mathrm{Tr}[\mathrm{T}_{X}\\{\mathrm{T}_{Y},\mathrm{T}_{Z}\\}]$
with $\mathrm{T}$ being the generators of the corresponding Lie Algebra. The
remaining anomaly cancellation conditions are given in Eq. 55. Overall, this
results in a series of equations for anomaly cancellation, restricting the
possible values of $X_{\Psi_{i}}$ from where one can extract all viable,
anomaly-free model versions.

### Anomaly-free solutions

Focusing on the solutions compatible with the BGL model, there are in total
three anomaly-free possibilities equipped with a type-I seesaw mechanism.
These are denoted as $\nu$BGL-I, $\nu$BGL-IIa and $\nu$BGL-IIb, being
characterized by their lepton sector textures as follows

  1. 1.

$\nu$BGL-I Scenario

$\displaystyle\begin{split}\Pi_{1},\Sigma_{1},B=\begin{pmatrix}\times&\times&\makebox[0.0pt][c]{0}\\\
\times&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix}\,,\quad\Pi_{2},\Sigma_{2}=\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix}\,,\\\
A=\mathbb{0}\,,\quad
C=\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\times&\times&\makebox[0.0pt][c]{0}\end{pmatrix}\,.\end{split}$ (30)

  2. 2.

$\nu$BGL-IIa Scenario

$\displaystyle\begin{split}\Pi_{1}&,\Sigma_{1}=\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix}\,,\\\
\Pi_{2}&=\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix}\,,\quad\Sigma_{2}=\begin{pmatrix}\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\end{pmatrix}\,\\\
A&=\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix}\,,\quad
B=\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix}\,,\quad
C=\mathbb{0}\,.\end{split}$ (31)

  3. 3.

$\nu$BGL-IIb Scenario

$\displaystyle\begin{split}A=\begin{pmatrix}\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\times\\\
\makebox[0.0pt][c]{0}&\times&\makebox[0.0pt][c]{0}\end{pmatrix}\,,\quad
B=\begin{pmatrix}\times&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\\\
\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}&\makebox[0.0pt][c]{0}\end{pmatrix}\,,\quad
C=\mathbb{0}\,.\end{split}$ (32)

where what distinguishes the $\nu$BGL-IIa and $\nu$BGL-IIb scenarios are the
textures in the Majorana neutrino sector. Last but not least we show in Tab. 1
the generic $\mathrm{U}(1)_{\mathrm{}}^{\prime}$ charges for each of the BGL
scenarios.

Charges |  $\nu$BGL-I |  $\nu$BGL-IIa |  $\nu$BGL-IIb |
---|---|---|---|---
$q_{L}$ | $\left[\begin{array}[]{c}x\\\ x\\\ x_{tL}\end{array}\right]$ | $--$ | $--$ |
$u_{R}$ | $\left[\begin{array}[]{c}y\\\ y\\\ x_{tR}\end{array}\right]$ | $--$ | $--$ |
$d_{R}$ | $\left[\begin{array}[]{c}2x-y\\\ 2x-y\\\ 2x-y\end{array}\right]$ | $--$ | $--$ |
$\ell_{L}$ | $\left[\begin{array}[]{c}-3x\\\ -3x\\\ 21x-6y\end{array}\right]$ | $\left[\begin{array}[]{c}x-y\\\ -7x+y\\\ 21x-6y\end{array}\right]$ | $\dfrac{1}{3}\left[\begin{array}[]{c}-x-2y\\\ -17x+2y\\\ 39x-12y\end{array}\right]$ |
$e_{R}$ | $\left[\begin{array}[]{c}-2x-y\\\ -2x-y\\\ 30x-9y\end{array}\right]$ | $\left[\begin{array}[]{c}2x-2y\\\ -6x\\\ 30x-9y\end{array}\right]$ | $\dfrac{1}{3}\left[\begin{array}[]{c}2x-5y\\\ -14x-y\\\ 58x-19y\end{array}\right]$ |
$\nu_{R}$ | $\left[\begin{array}[]{c}-4x+y\\\ -4x+y\\\ 12x-3y\end{array}\right]$ | $\left[\begin{array}[]{c}0\\\ -8x+2y\\\ 12x-3y\end{array}\right]$ | $\dfrac{1}{3}\left[\begin{array}[]{c}-4x+y\\\ -20x+5y\\\ 20x-5y\end{array}\right]$ |
$\Phi$ | $\left[\begin{array}[]{c}-x+y\\\ -9x+3y\end{array}\right]$ | $\left[\begin{array}[]{c}-x+y\\\ -9x+3y\end{array}\right]$ | $\dfrac{1}{3}\left[\begin{array}[]{c}3(-x+y)\\\ -19x+7y\end{array}\right]$ |
$S$ | $8x-2y$ | $-4x+y$ | $\dfrac{8x-2y}{3}$ |
Table 1: Allowed charges for the various models. For model $\nu$BGL-I and IIa
we have $x_{tL}=-7x+2y$ and $x_{tR}=-16x+5y$. Model $\nu$BGL-IIb has
$x_{tL}=(-13x+4y)/3$ and $x_{tR}=(-32x+11y)/3$. In order for the BGL textures
to be preserved, we additionally require that $y\neq 4x$.

##  IV Chosen scenario

In what follows we will focus on the $\nu$BGL-I scenario, whose allowed
charges are constrained by the first column in Tab. 1. We make this choice
after having investigated the other possibilities available within this model
and verified that the $\nu$BGL-I scenario is the one for which our numerical
analysis and parameter space scans were more efficient. This, of course, does
not exclude the possible interest of the other scenarios, but our main goal is
to understand the typical phenomenology of this class of models. Thus, for
$\nu$BGL-I, choosing $x=1$ and $y=1/3$ the model considered in this article is
then defined by the quantum numbers in Tab. 2.

| $\Phi_{1}$ | $\Phi_{2}$ | $S$ | $q_{1}$ | $q_{2}$ | $q_{3}$ | $u_{R_{1}}$ | $u_{R_{2}}$ | ${u}_{R_{3}}$ | $d_{R_{1}}$ | $d_{R_{2}}$ | $d_{R_{3}}$
---|---|---|---|---|---|---|---|---|---|---|---|---
$\mathrm{U(1)_{Y}}$ | 1/2 | 1/2 | 0 | 1/6 | 1/6 | 1/6 | 2/3 | 2/3 | 2/3 |  $-$1/3 |  $-$1/3 |  $-$1/3
$\mathrm{SU(2)_{L}}$ | 2 | 2 | 1 | 2 | 2 | 2 | 1 | 1 | 1 | 1 | 1 | 1
$\mathrm{SU(3)_{C}}$ | 1 | 1 | 1 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3
$\mathrm{U(1)^{\prime}}$ |  $-$2/3 |  $-$8 | 22/3 | 1 | 1 |  $-$19/3 | 1/3 | 1/3 |  $-$43/3 | 5/3 | 5/3 | 5/3

| $\ell_{1}$ | $\ell_{2}$ | $\ell_{3}$ | $e_{R_{1}}$ | $e_{R_{2}}$ | $e_{R_{3}}$ | $\nu_{R_{1}}$ | $\nu_{R_{2}}$ | $\nu_{R_{3}}$
---|---|---|---|---|---|---|---|---|---
$\mathrm{U(1)_{Y}}$ |  $-$1/2 |  $-$1/2 |  $-$1/2 |  $-$1 |  $-$1 |  $-$1 | 0 | 0 | 0
$\mathrm{SU(2)_{L}}$ | 2 | 2 | 2 | 1 | 1 | 1 | 1 | 1 | 1
$\mathrm{SU(3)_{C}}$ | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1
$\mathrm{U(1)^{\prime}}$ |  $-$3 |  $-$3 | 19 |  $-$7/3 |  $-$7/3 | 27 |  $-$11/3 |  $-$11/3 | 11

Table 2: Charge assignments of the SU(2) scalar doublets and SM fermions under
the flavour symmetry $\mathrm{U(1)^{\prime}}$ and the SM gauge group,
$\mathrm{SU(3)_{C}\times SU(2)_{L}\times U(1)_{Y}}$. The corresponding
$\mathrm{U(1)^{\prime}}$ charges are computed using the formulas from the
first column of Tab. 1 taking $x=1$ and $y=1/3$.

Notice that the assigned quantum numbers for the up-type quarks can also be
flipped in such a way that there are in total three variations of the
$\nu$BGL-I scenario. This depends, in practice, on the choice of up-type quark
associated with the $\Delta_{2}$ matrix in Eq. (2) and the designation of the
model is typically associated with that, i.e. variation $u$, $c$ or $t$. In
the reminder of this article we focus on variation $t$.

As already stated, one of the novel features of the model considered is the
presence of three generations of right-handed neutrinos. In order to comply
with experimental data a realistic lepton mixing and neutrino mass hierarchy
must be assured. Both the charged neutral sectors can be rotated to the mass
basis via bi-unitary matrices, that is,

$\displaystyle
M_{\mathrm{L}}^{\mathrm{diag}}=U^{\dagger}_{\mathrm{L}}M_{\mathrm{L}}U_{\mathrm{R}},$
(33) $\displaystyle
M_{\mathrm{\nu}}^{\mathrm{diag}}=U^{\dagger}_{\nu}M_{\nu}U_{\nu},$

where $M_{\mathrm{L}}$ is a $3\times 3$ matrix expressed in the basis
$\\{e_{\mathrm{L}}^{1},e_{\mathrm{L}}^{2},e_{\mathrm{L}}^{3}\\}\otimes\\{e_{\mathrm{R}}^{1},e_{\mathrm{R}}^{2},e_{\mathrm{R}}^{3}\\}$
as usual, while $M_{\mathrm{\nu}}$ is a $6\times 6$ matrix written in the
basis
$\\{\nu_{\mathrm{L}}^{1},\nu_{\mathrm{L}}^{2},\nu_{\mathrm{L}}^{3},\nu_{\mathrm{R}}^{1},\nu_{\mathrm{R}}^{2},\nu_{\mathrm{R}}^{3}\\}$.
On the other hand, $U_{\mathrm{L}}$, $U_{\mathrm{R}}$ and $U_{\nu}$ are
unitary matrices. We can then define the PMNS matrix as

$V_{P}=U^{\dagger}_{\mathrm{L}}PU_{\nu}=\begin{pmatrix}V_{\mathrm{PMNS}}^{\mathrm{SM}}&\mathcal{O}_{P}\end{pmatrix},\quad
P=\begin{pmatrix}\mathbb{1}_{3\times 3}\hskip 5.0pt|\hskip 5.0pt0_{3\times
3}\end{pmatrix},$ (34)

where $V_{P}$ is a $3\times 6$ matrix with the first $3\times 3$ entries
encoding the SM charged-lepton and neutrinos interactions and
$\mathcal{O}_{P}$ contains the interactions between leptons and the heavy
neutrinos. In our analysis we require $\mathcal{O}_{P}<10^{-5}$ such that
charged currents involving active and Beyond Standard Model (BSM) neutrinos
are suppressed. Also, we require the masses of such active neutrinos to be
between [0,$10^{-9}$] GeV while their heavy counterparts have masses around
the TeV scale. While an in-depth study of neutrino physics is beyond the scope
of this work, we have performed a fit to the light neutrinos masses in order
to respect the normal hierarchy Zyla  _et al._ (2020a),

$\begin{split}M_{\nu_{2}}^{2}-M_{\nu_{1}}^{2}\approx
10^{-5}~{}\mathrm{eV^{2}},\\\ M_{\nu_{3}}^{2}-M_{\nu_{2}}^{2}\approx
10^{-3}~{}\mathrm{eV^{2}}.\end{split}$ (35)

This was numerically required for all points that pass the theoretical and
experimental constraints, as it will be discussed in upcoming sections.
Although a full fit was not done across the entirety of the parameter space,
and given that our studies of the scalar and quark sectors are independent of
neutrino physics, we have confirmed and can safely argue that the discussed
model is compatible with neutrino data.

### Inversion procedure for the scalar sector

Following the same methodology used previously in Das  _et al._ (2021), our
goal is to provide the physical masses and mixing angles as inputs and then
search for valid regions in the parameter space. In the considered $\nu$BGL-I
scenario the scalar potential is given in Eq. 16 with $a_{3},a_{4}=0$.
Selecting a vacuum configuration with real VEVs $(v_{1},v_{2},v_{S})$, one
obtains the following minimization conditions

$\displaystyle\begin{split}&\mu_{1}^{2}=\quad-\left(v_{1}^{2}\lambda_{1}+\frac{v_{S}^{2}\lambda^{\prime}_{2}}{2}+\frac{v_{2}\hat{\mu}}{2v_{1}}\right)\>,\\\
&\mu_{2}^{2}=\quad-\left(v_{2}^{2}\lambda_{2}+\frac{v_{S}^{2}\lambda^{\prime}_{3}}{2}+\frac{v_{1}\hat{\mu}}{2v_{2}}\right)\>,\\\
&\mu_{S}^{2}=\quad-\frac{\sqrt{2}v_{1}v_{2}\left(a_{1}+a_{2}\right)+v_{S}\left(v_{1}^{2}\lambda_{2}^{\prime}+v_{2}^{2}\lambda^{\prime}_{3}+2v_{S}^{2}\lambda^{\prime}_{1}+2\mu_{b}^{2}\right)}{2v_{S}}\,,\end{split}$
(36)

where one defines

$\displaystyle\begin{split}\hat{\mu}=v_{S}\sqrt{2}\left(a_{1}+a_{2}\right)+v_{1}v_{2}\left(\lambda_{3}+\lambda_{4}\right)+2\mu_{3}^{2}\,.\end{split}$
(37)

The physical scalar spectrum contains three CP-even scalars $H_{1,2,3}$, two
pseudo-scalars $A_{2,3}$, and one charged Higgs boson $H^{\pm}$. In a
configuration close to the vacuum alignment in one of the Higgs doublets, the
corresponding CP-even scalar corresponds to the SM Higgs boson with mass
$125.09~{}\mathrm{GeV}$ and with SM-like interactions to the gauge bosons and
fermions. Notice that, if such interactions coincide with the SM ones, one
often refers to it as the exact alignment limit. However, direct searches for
new scalars Aaboud  _et al._ (2018b); Khachatryan  _et al._ (2016a, b); Aad
_et al._ (2020b, 2021b); Sirunyan  _et al._ (2020b) still keep a window opened
for possible off-alignment deviations that we consider in our analysis below.
In order to keep the physical parameters of the theory under control one
performs an inversion procedure analogous to that developed in Das  _et al._
(2021), expressing the parameters of the scalar potential in Eq. 36 in terms
of physical masses, mixing angles (including the Higgs alignment parameter)
and VEVs as we discuss in what follows.

First, the matrix $\mathcal{O}_{\beta}$, which rotates the gauge eigenstates
into the Higgs basis333The Higgs basis corresponds to a redefinition of the
physically-equivalent doublets so that only one of them contains a VEV., is
introduced

$\displaystyle\begin{split}\mathcal{O}_{\beta}=\begin{pmatrix}\cos\beta&\sin\beta&0\\\
-\sin\beta&\cos\beta&0\\\
0&0&1\end{pmatrix}=\begin{pmatrix}v_{1}/v&v_{2}/v&0\\\ -v_{2}/v&v_{1}/v&0\\\
0&0&1\end{pmatrix}\end{split}$ (38)

where one defines

$\displaystyle\begin{split}&v_{1}=v\cos\beta\,,\\\ &v_{2}=v\sin\beta\,,\\\
&v=\sqrt{{v_{1}}^{2}+{v_{2}}^{2}}\approx 246\ \mathrm{GeV}\,.\end{split}$ (39)

Starting with the scalar potential in Eq. 15 and expanding it using the field
definitions in Eq. 3, the CP-odd sector mass matrix can be obtained as

$\displaystyle\frac{1}{2}\begin{pmatrix}I_{1}&I_{2}&\eta\end{pmatrix}{\mathcal{M}_{S}}^{2}\begin{pmatrix}I_{1}\\\
I_{2}\\\ \eta\end{pmatrix}\,.$ (40)

Using $\mathcal{O}_{\beta}$ one can rotate $\mathcal{M}_{S}^{2}$ to the Higgs
basis, where it acquires a block-diagonal form

$\displaystyle{\mathfrak{M}_{S}}^{2}=\mathcal{O}_{\beta}\mathcal{M}_{S}^{2}\mathcal{O}_{\beta}^{T}\,.$
(41)

We can further diagonalize ${\mathfrak{M}_{S}}^{2}$ via an additional rotation
matrix $\mathcal{O}_{\gamma}$ as

$\displaystyle\begin{pmatrix}0&0&0\\\ 0&{M_{A_{2}}}^{2}&0\\\
0&0&{M_{A_{3}}}^{2}\end{pmatrix}=\mathcal{O}_{\gamma}{\mathfrak{M}_{S}}^{2}{\mathcal{O}_{\gamma}}^{T}=\mathcal{O}_{\gamma}\mathcal{O}_{\beta}\mathcal{M}_{S}^{2}\mathcal{O}_{\beta}^{T}{\mathcal{O}_{\gamma}}^{T}$
(42)

where one defines

$\displaystyle\mathcal{O}_{\gamma}=\begin{pmatrix}1&0&0\\\
0&\cos\gamma_{1}&-\sin\gamma_{1}\\\
0&\sin\gamma_{1}&\cos\gamma_{1}\end{pmatrix}\,.$ (43)

Notice that while the neutral CP-odd Goldstone boson can be gauged away from
the particle spectrum and gets absorbed as the longitudinal polarisation of
the $\mathrm{Z^{0}}$ boson, one of the two physical pseudo-scalars is in fact
a pseudo-Goldstone state that emerges due to both spontaneous and explicit
breaking of the global $\mathrm{U}(1)_{\mathrm{}}^{\prime}$ symmetry.

Similarly, the CP-even states can be rotated to the mass basis as

$\displaystyle\begin{pmatrix}{M_{H_{1}}}^{2}&0&0\\\ 0&{M_{H_{2}}}^{2}&0\\\
0&0&{M_{H_{3}}}^{2}\end{pmatrix}=\mathcal{O}_{\alpha}\mathcal{O}_{\beta}\mathcal{M}_{H}^{2}\mathcal{O}_{\beta}^{T}{\mathcal{O}_{\alpha}}^{T}$
(44)

with

$\displaystyle\mathcal{O}_{\alpha}=R_{3}\dotproduct R_{2}\dotproduct R_{1}$
(45)

and

$\displaystyle\begin{split}R_{1}=\begin{pmatrix}\cos{\alpha_{1}}&-\sin{\alpha_{1}}&0\\\
\sin{\alpha_{1}}&\cos{\alpha_{1}}&0\\\ 0&0&1\\\
\end{pmatrix},R_{2}=\begin{pmatrix}\cos{\alpha_{2}}&0&\sin{\alpha_{2}}\\\
0&1&0\\\ -\sin{\alpha_{2}}&0&\cos{\alpha_{2}}\\\
\end{pmatrix},R_{3}=\begin{pmatrix}1&0&0\\\
0&\cos{\alpha_{3}}&\sin{\alpha_{3}}\\\
0&-\sin{\alpha_{3}}&\cos{\alpha_{3}}\end{pmatrix}\,.\end{split}$ (46)

The presence of a CP-even state behaving as a SM-like Higgs boson is ensured
not only requiring that

$\displaystyle M_{H_{1}}=125.09\ \textrm{GeV}\ ,$ (47)

but also that the couplings to fermions and gauge bosons of this particle
closely reproduce the SM predictions, which are controlled by the coupling
modifier

$\left[\mathcal{O}_{\alpha}\mathcal{O}_{\beta}\right]_{12}=-\cos\alpha_{2}\sin\delta\qquad\textrm{with}\qquad\delta\equiv\alpha_{1}-\beta$
(48)

such that the coupling of the 125 GeV Higgs boson to the electroweak gauge
bosons can be expressed as

$\displaystyle\lambda^{\prime}_{h\,V\,V}=-\lambda^{SM}_{h\,V\,V}\cos\alpha_{2}\sin\delta$
(49)

with $\lambda^{SM}_{h\,V\,V}$ denoting the SM Higgs boson coupling to gauge
bosons, $V\,=\,Z,\,W$ in the SM. This means that the alignment condition
concerning Higgs gauge interactions is governed not by a single angle but
instead two angular parameters, $\delta$ and $\alpha_{2}$. Furthermore, as the
LHC has been putting stringent limits on deviations from the SM, one expects
that $\cos\alpha_{2}\sin\delta\sim\mathcal{O}(1)$, i.e. $\delta$ must not be
far from $\pi/2$ while $\alpha_{2}$ must be small. Note also that Eq. (48)
ensures that the interactions of the SM-like Higgs boson are dominantly
inherited from the second doublet which is the one that couples to the top-
quark. Further conditions for alignment will be obtained from the Higgs’
interactions with fermions. For instance, given that the top mass in this
model comes only from the matrix $\Delta_{2}$, it will be proportional to
$v_{2}$ and therefore the coupling modifier for the Higgs-top interations will
be

$\displaystyle\lambda^{\prime}_{h\,t\,\bar{t}}=-\lambda^{SM}_{h\,t\,\bar{t}}\frac{\cos\alpha_{2}\sin\delta}{\sin\beta}\,.$
(50)

With $\sin\delta=\sin(\alpha_{1}-\beta)\sim\mathcal{O}(1)$ as mentioned
earlier, we see that a further condition for alignment are
$\cos\alpha_{2}\simeq-\sin\beta\simeq-\cos\alpha_{1}$. Additional conditions
would be obtained from analysing the remaining Higgs-fermion interactions.
Notice that the minus signs present in Eqs. (48) and (50) are not necessarily
meaningful. We know, from the measured Higgs diphoton branching ratio at the
LHC, that the relative sign of $\lambda^{\prime}_{h\,t\,\bar{t}}$ and
$\lambda^{\prime}_{h\,W\,W}$ ought to be the same as that of
$\lambda^{SM}_{h\,t\,\bar{t}}$ and $\lambda^{SM}_{h\,W\,W}$. An overall minus
sign affecting both couplings can therefore be physically meaningless. In our
numerical analysis of the next section, however, we took into account all
possible relative signs, and compared the predicted Higgs branching ratios and
production cross sections with the known experimental values, thus ensuring
that, whatever relative signs there may be between Higgs couplings in this
model, the LHC constraints are satisfied.

Last but not least, the charged Higgs sector contains only one physical scalar
as the charged Goldstone boson is gauged away from the particle spectrum and
plays a role of the longitudinal polarization of the $\mathrm{W}^{\pm}$ boson.
In particular, the mass form reads as

$\displaystyle\begin{pmatrix}0&0\\\ 0&{M_{H^{\pm}}}^{2}\\\
\end{pmatrix}={\mathcal{O}_{\gamma}}^{\prime}{\mathcal{O}_{\beta}^{\prime}}\mathcal{M}^{\pm}{\mathcal{O}_{\beta}^{\prime}}^{T}{{\mathcal{O}_{\gamma}}^{\prime}}^{T}$
(51)

with

$\displaystyle\mathcal{O}_{\beta}^{\prime}=\begin{pmatrix}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{pmatrix}\ ,\
\mathcal{O}_{\gamma}^{\prime}=\begin{pmatrix}\cos\gamma_{2}&\sin\gamma_{2}\\\
-\sin\gamma_{2}&\cos\gamma_{2}\end{pmatrix}\,.$ (52)

The invertion procedure described above is applied to the scalar sector with
the purpose of trading the quartic couplings for physical masses and mixing
angles. The corresponding analytical formulas are lengthy and can be found in
Appendix B. In our numerical analysis, the physical masses
$M_{H_{1}},M_{H_{2}},M_{H_{3}},M_{A_{2}},M_{A_{3}},M_{H^{\pm}}$, the angles
$\beta,\alpha_{1},\alpha_{2},\alpha_{3},\gamma_{1},\gamma_{2}$ and the cubic
couplings $a_{1},a_{2},a_{3}$ are then used as input parameters limited to be
randomly sampled within physically reasonable ranges as detailed below in Sec.
V.

##  V Parameter space and phenomenological constraints

The parameter space scan discussed in this section relies on the inversion
procedure described previously, where both scalar masses and mixing angles are
used as input, allowing one to determine the gauge eigenbasis couplings of the
Higgs potential in Eq. 16. The ranges used in the scan can be found in Tab. 3.

Parameter |  $\alpha_{2}$,$\alpha_{3}$,$\gamma_{1}$ | $\tan\beta$ | $\delta$ | $a_{3},a_{4}$
---|---|---|---|---
range | $[-\pi,\pi]$ | $[0.5,30]$ | $[\tfrac{\pi}{2}-1,\tfrac{\pi}{2}+1]$ | $0$

Parameter |  $M_{A_{2}}$, $M_{H^{\pm}}$ | $M_{A_{3}}$ |  $M_{H_{2}}$,$M_{H_{3}}$ | $a_{2}$
---|---|---|---|---
range [GeV] | $[0.5,1600]$ | $[30,2000]$ | $[126,1800]$ | $[-1,1]$

Table 3: Parameter intervals used in numerical scans.

According to the discussion below Eq. 49, the chosen scanning ranges for
$\delta$ are centered around $\pi/2$ where $\sin\delta=1$. While $\alpha_{2}$
must be small, we let it range between $-\pi$ and $\pi$ in order to let
$\cos\alpha_{2}$ capture both positive and negative values. Furthermore, the
somewhat large ranges set for both $\delta$ and $\alpha_{2}$ were deliberately
picked in order to ensure that the allowed off-alignment regions are entirely
covered in the numerical analysis, as we discuss further ahead. It is also
worth mentioning that the SM-like Higgs boson was chosen to be the lightest of
the CP-even scalars and the absolute value of the cubic couplings is such that
$\absolutevalue{a_{2}}\leq 1$.

Using the inversion procedure we then extract the Lagrangian parameters and
feed them to SPheno Porod and Staub (2012) in order to calculate both the STU
electroweak precision observables (or oblique parameters) as well as Higgs
decay widths and branching fractions. The most up to date electroweak fit for
the STU parameter Zyla  _et al._ (2020b) can be summarized as

$\displaystyle\begin{split}\begin{matrix}S&=-0.01\pm 0.10\\\ T&=0.03\pm
0.12\\\ U&=0.02\pm 0.11\end{matrix}\ ,\
\quad\rho_{ij}=\begin{pmatrix}1&0.92&-0.80\\\ 0.92&1&-0.93\\\ -0.80&-0.93&1\\\
\end{pmatrix}\end{split}$ (53)

with $\rho_{ij}$ denoting a $92\%$ correlation between S and T, while U-S and
U-T are $80\%$ and $-93\%$ anti-correlated, respectively. For the purpose of
this analysis, we required $\Delta\chi^{2}<7.815$, which translates into a 95%
confidence level (CL) agreement with the electroweak fit, where

$\displaystyle\Delta\chi^{2}=\sum_{ij}\left(\Delta\mathcal{O}_{i}-\Delta\mathcal{O}_{i}^{(0)}\right)\left[(\sigma^{2})^{-1}\right]_{ij}\left(\Delta\mathcal{O}_{j}-\Delta\mathcal{O}_{j}^{(0)}\right)$
(54)

with $\left[\sigma^{2}\right]_{ij}\equiv\sigma_{i}\rho_{ij}\sigma_{j}$ the
covariance matrix and $\sigma_{i}$ the standard deviations in Eq. 53. The
notation employed here is such that $\Delta\mathcal{O}_{i}^{(0)}$ indicates
the $S,T,U$ experimental fit values, while $\Delta\mathcal{O}_{i}$ denotes the
$S,T,U$ calculated values in SPheno. In Fig. 1 we show the points that pass
the STU restrictions with a CL of, at least, 95% (blue shades) in comparison
with the points that are excluded (grey). The colour bar represents the
$\Delta\chi^{2}$ values. The majority of points that were simulated are
rejected due to large $T$ values, as anticipated. However, one of the benefits
of 2HDM-like models is that there are parameter space regions for which all
constraints on the oblique parameters are satisfied, as is the case here.

Figure 1: Electroweak precision observables for all simulated points. The
colored points are those who pass the STU analysis with a confidence level
(CL) of at least 95%. Grey points are excluded by precision EW fit data.

While the STU analysis offers a strong razor of the parameter space, it is
still necessary to match these points to both Higgs and flavour observables.
First, let us focus on the impact of Higgs physics leaving the flavour
analysis for later. In Fig. 2, we show the allowed regions in the
$\delta$-$\alpha_{2}$ plane that survive not only the EW precision tests but
also that comply with the properties of a SM-like Higgs boson and are not
excluded by direct searches for new scalars at the LHC. As usual, this is done
by using the SPheno output SLHA files as input cards to HiggsSignals (HS) and
HiggsBounds (HB) Bechtle  _et al._ (2015, 2014) for every single point that
passes the STU test.

Figure 2: Angles entering in the Higgs coupling modifier of Eq. 49. Grey
points are excluded by EW precision observables while orange ones pass the STU
$\Delta\chi^{2}$ test but fail Higgs physics constraints. Red points pass all
restrictions and represent the regions where the Higgs boson is closely
aligned with a SM-like one.

The points that survive both HS and HB are shown in red and were found to pile
in a region around that of the exact alignment limit where $\delta=\pi/2$ and
$\alpha_{2}=-\pi\,,0\,,\pi$. In particular, we found that deviations from an
exact alignment are primarily driven by the $\delta$ angle which can vary in a
wider range, whereas $\alpha_{2}$ is tightly constrained in the aforementioned
regions, with $0.60<-\cos\alpha_{2}\sin\delta<1.0$ for the upper red island,
$0.60<\cos\alpha_{2}\sin\delta<0.99$ in the middle region, and
$0.60<-\cos\alpha_{2}\sin\delta<0.98$ at the $\alpha_{2}=-\pi$ line. In Fig. 3
we show the dependence of the charged Higgs boson mass $M_{H^{\pm}}$ in terms
of the masses of the $A_{2}$ and $A_{3}$ pseudoscalars (left panels) as well
as of the masses of the non-SM CP-even states $H_{2}$ and $H_{3}$ (right
panels). The colour scheme is the same as in Fig. 2 and, as the red region
demonstrates, we have obtained valid points that cover a wide range of masses
for new scalars. In particular, there are viable solutions with charged Higgs
bosons that can be as light as approximately $100~{}\mathrm{GeV}$ extending
all the way up to about $1.6~{}\mathrm{TeV}$, the light pseudo-scalar is
limited to the range that goes from $100~{}\mathrm{GeV}$ up to
$1.5~{}\mathrm{TeV}$ while its heavy counterpart is heavier than
$300~{}\mathrm{GeV}$ but lighter than $1.8~{}\mathrm{TeV}$. For the CP-even
sector, besides a SM-like state with mass $125.09~{}\mathrm{GeV}$, one finds
that $H_{2}$ is sharply limited from below at the Higgs boson mass, i.e.
$125.09~{}\mathrm{GeV}$, allowing valid solutions up to approximately
$1.4~{}\mathrm{TeV}$, while its heavy partner mass ranges between
$200~{}\mathrm{GeV}$ and approximately $1.7~{}\mathrm{TeV}$. Note that the
limits above must be understood in the context of a search of the parameter
space dedicated to find light scalars. In fact, the low mass region is
particularly interesting from the collider physics perspective and can
potentially be tested in the near future.

Figure 3: Mass of the charged Higgs versus the masses of the pseudoscalars
$A_{2}$ and $A_{3}$ (left panels) and of the CP-even scalars $H_{2}$ and
$H_{3}$ (right panels). Grey points are excluded by the oblique STU
parameters, orange points are excluded either by HB or HS, while
simultaneously passing the STU restrictions. Red points pass all constraints.

Besides electroweak and Higgs physics constraints, one also has to confront
the generated data with quark flavour violation (QFV) observables. According
to the textures of the chosen $\nu$BGL-I scenario, the most stringent
constraints are on the down-type quark sector rather than in the up-quark or
charged lepton ones. Notice also that the most stringent QFV restrictions are
indeed coming from processes involving down-type quarks, in particular those
shown in Tab. 4.

Channel | $\mathcal{O}_{SM}$ | $\sigma_{SM}$ | $\mathcal{O}_{Exp}$ | $\sigma_{Exp}$ | $\sigma$
---|---|---|---|---|---
$\mathrm{BR(B\rightarrow\chi_{s}\gamma)}$ | $3.29\times 10^{-4}$ | $1.87\times 10^{-5}$ | $3.32\times 10^{-4}$ | $0.16\times 10^{-4}$ | 0.075
$\mathrm{BR(B_{s}\rightarrow\mu\mu)}$ | $3.66\times 10^{-9}$ | $1.66\times 10^{-10}$ | $2.80\times 10^{-9}$ | $0.06\times 10^{-9}$ | 0.038
$\Delta M_{d}$ (GeV) | $3.97\times 10^{-13}$ | $5.07\times 10^{-14}$ | $3.33\times 10^{-13}$ | $0.013\times 10^{-13}$ | 0.11
$\Delta M_{s}$ (GeV) | $1.24\times 10^{-11}$ | $7.08\times 10^{-13}$ | $1.17\times 10^{-11}$ | $0.0014\times 10^{-11}$ | 0.054
$\epsilon_{K}$ (GeV) | $1.81\times 10^{-3}$ | $2.00\times 10^{-4}$ | $2.23\times 10^{-3}$ | $0.011\times 10^{-3}$ | 0.14

Table 4: Most relevant Quark Flavour Violation (QFV) observables. The first
column represents the prediction of the SM for each observable, the second
column indicates the 1$\sigma$ error of the SM theoretical calculation, the
third column is the measured experimental value, the fourth is the 1$\sigma$
experimental error and the last column is this measured deviation between the
SM prediction and the observed value from experiment.

We show in Fig. 4 the individual impact that each of the considered QFV
observable poses to the points that have survived the electroweak and Higgs
physics constraints.

Figure 4: Mass of the charged Higgs versus the mass of the lightest
pseudoscalar $A_{2}$. Grey points are excluded by STU observables, orange
points are excluded by HS or HB while still passing STU, and red points pass
STU, HS and HB constraints. In green, we showcase the points that pass a given
QFV observable, namely, we have a) $B\rightarrow\chi_{S}\gamma$, b) $\Delta
M_{S}$, c) $\epsilon_{K}$, d) $B_{S}\rightarrow\mu\mu$ and e) $\Delta M_{d}$,
while passing STU, HB and HS constraints.

The colour scheme is identical to that in Figs. 2 and 3 with the inclusion of
green points that represent points simultaneously allowed by EW, Higgs and QFV
observables. It is promptly noticeable that both $\epsilon_{K}$, panel c), and
$\mathrm{BR}(B\rightarrow\chi_{s}\gamma)$, panel a), represent the weaker
flavour physics constraints on the parameter space in such a way that all of
the points that have survived the EW and Higgs physics constraints are also
compatible with these two QFV observables. On the other hand, $\Delta M_{s}$,
panel b), $\Delta M_{d}$, panel e), and $\mathrm{BR}(B_{s}\rightarrow\mu\mu)$,
panel d), place a noticeable veto on the allowed parameter space. $\Delta
M_{s}$ is the main culprit, with the lowest number of allowed points that
survive such a constraint. We also note that, both $\Delta M_{d}$ and $\Delta
M_{s}$ limit the mass to be below $1~{}\mathrm{TeV}$ for the $A_{2}$
pseudoscalar. Indeed, the low mass regime is where the vast majority of the
allowed points tend to accumulate. For a cleaner analysis we shown in the
histograms of Fig. 5 the number of points that are allowed by $\Delta M_{s}$,
$\Delta M_{d}$, $\mathrm{BR}(B_{s}\rightarrow\mu\mu)$ and combinations of
them, in bins of the lightest pseudoscalar mass.

Figure 5: Histograms containing points that survive STU, HS, HB and a given
QFV (or pair of) in bins of the $A_{2}$ mass. The most restrictive is coloured
in  blue.

We also show in Tab. 5 the acceptance ratio, defined as the percentage of the
initial sample consistent with EW and Higgs physics that also survives a
certain QFV cut, for each of the considered QFV observables as well as the
three combinations in Fig. 5.

Set of QFV observables | Acceptance ratio
---|---
$\mathrm{BR}\left(B\rightarrow\chi_{s}\gamma\right)$ | 100.0%
$\mathrm{BR}\left(B_{s}\rightarrow\mu\mu\right)$ | 35.0%
$\Delta M_{d}$ (GeV) | 48.0%
$\Delta M_{s}$ (GeV) | 26.0%
$\epsilon_{K}$ (GeV) | 100.0%
$\mathrm{BR}\left(B_{s}\rightarrow\mu\mu\right)$ & $\Delta M_{s}$ | 9.39%
$\mathrm{BR}\left(B_{s}\rightarrow\mu\mu\right)$ & $\Delta M_{d}$ | 22.21%
$\Delta M_{s}$ & $\Delta M_{d}$ | 25.57%
Table 5: Acceptance ratio for each QFV (first five rows) and for pairs of QFV
variables (last three rows). The most restrictive set is marked in  blue.

Our results reveal that the sharpest razor of the parameter space is the
combination of $\mathrm{BR}(B_{s}\rightarrow\mu\mu)$ and $\Delta M_{s}$ with a
very low acceptance ratio of only 9.39%, while their individual contributions
are $35\%$ and $26\%$, respectively. This combination of observables also
prefers pseudoscalar masses to be in a range between $125~{}\mathrm{GeV}$ and
$1300~{}\mathrm{GeV}$. However, such a limited range must not be regarded as a
particular feature of the $\nu$BGL-I model but merely as a result of the
numerical scan, which is only efficient if all scalar masses are sufficiently
light444When using an invertion procedure that receives as input parameters
only physical masses and mixing angles, large hierarchies between the Higgs
boson mass and the remaining scalars often result in unacceptably large values
for the quartic couplings, which can be as large as several tenths of
thousands.. Indeed, it is well-known that multi-Higgs models typically observe
a decoupling limit when the masses of the new scalars tend to infinity.
Therefore, the study of parametric regions with large scalar masses requires a
different strategy to increase the acceptance ratio as, e.g. it was done in
Das  _et al._ (2021) by forcing the exact alignment limit and providing only
the Higgs boson mass as input.

Figure 6: Plots of the ratio between the BR for the decays
$B_{S}\rightarrow\mu^{+}\mu^{-}$ (on the left) and $B\rightarrow X_{s}\gamma$
(on the right) computed in the $\nu$BGL-I model and in the SM versus the mass
of the charged scalar. In the colour bar, we showcase the mass of the heaviest
pseudo-scalar $A_{3}$. Grey points are excluded from at least one QFV
observable from 4 at a 2$\sigma$ uncertainty. Dashed blue lines indicate the
$2\sigma$ bound and full blue lines indicate the 1$\sigma$ limit.

Figure 7: Summary plots for the QFV variables. Dark grey points are excluded
for at least one QFV observable in Table 4 while light grey points do not pass
the EW precision tests. Orange points are excluded either by HB or HS and
green (cyan) points are compatible with EW and Higgs physics as well as all
QFV constraints at a 2$\sigma$ (1$\sigma$) level.

For a closer inspection we show in Fig. 6 two of the QFV observables of Tab. 4
in terms of the charged Higgs and heavier pseudoscalar masses. For a moment we
ignore the EW and Higgs physics constraints and show in grey the points that
are excluded by at least one QFV observable while the colored points are
simultaneously allowed by all the QFV restrictions. Note that, even
considering the entire universe of sampled points there are only two points
outside the $\pm 2\sigma$ uncertainty band (dashed lines) that would fail
$\mathrm{BR}\left(B\rightarrow\chi_{s}\gamma\right)$ (right panel) while
$\mathrm{BR}\left(B_{s}\rightarrow\mu\mu\right)$ is visibly rather more
restrictive (left panel). It is also clear that larger scalar masses tend to
relax the QFV constraints as expected.

Last but not least, we show in Fig. 7 the combination of all EW, Higgs and
flavour physics constraints. We have considered the four mostly stringent QFV
observables in Tab. 4, leaving out only $\epsilon_{K}$ as it adds the least
amount of constraints to our model. It is clear that the majority of sampled
points fall in the dark grey area where at least one of the considered QFV
observables lies beyond the $2\sigma$ uncertainty range limited by the blue
dashed lines. As one can also see, flavour physics plays the major role in
strongly reducing the allowed parameter space. In particular, for points
complying with a $2\sigma$ uncertainty of all QFV observables (within the two
blue dashed curves), we note that the majority of allowed cases also prefer
the charged Higgs boson mass to be in a range between $400$ GeV and
$1.4~{}\mathrm{TeV}$. If we further demand that all QFV observables are
simultaneously reproduced within their 1$\sigma$ bound, then the viable
region, identified by the cyan points, tends to be further compressed. Once
again, it is important to mention that our scan was optimized to the low mass
region, becoming increasingly inefficient with increasing scalar masses.
Therefore one cannot extract bold conclusions about the high mass regime. On
the contrary, what is important to retain as a key result of this article is
that the $\nu$BGL-I model is compatible with the current EW, Higgs and QFV
observables in the limit where new scalar masses are near the EW-scale. This
opens up the possibility for testing the model in a foreseeable future,
potentially already at the LHC run-III, as discussed below.

##  VI Phenomenology at the LHC

In this section we compare the points that have survived all constraints with
current direct searches for neutral scalars at the LHC, focusing on $\tau\tau$
Aad  _et al._ (2020c), $ZZ$ Aad  _et al._ (2021b) and $WW$ Aaboud  _et al._
(2018c) final states. An interface between SPheno and MadGraph Alwall  _et
al._ (2014) allows one to determine the neutral scalar production cross
section times its decay branching ratio to a given final state.

In Fig. 8 we show the latest ATLAS limits for scalar and pseudoscalar searches
in the di-tau final state and confront them with our prediction for the
lighter BSM Higgs bosons $H_{2}$ and $A_{2}$. The green (cyan) points have the
same meaning as in Fig. 7, i.e. they are points for which the measured values
of the considered QFV observables are reproduced in our calculation within a
2$\sigma$ (1$\sigma$) uncertainty.

Figure 8: The signal strength for production of a CP-even (top panels) and a
CP-odd (bottom pannels) scalar via the gluon-gluon fusion mechanism (a), (c)
times its branching ratio to $\tau\bar{\tau}$, and (b), (d) with associated
production of a $b\bar{b}$-pair times its branching fraction to
$\tau\bar{\tau}$, as a function of the lightest CP-even mass, $M_{H_{2}}$ (top
panels) and lightest pseudoscalar mass $M_{A_{2}}$ (bottom panels). Green
points satisfy QFV observables within 2 standard deviations, whereas blue ones
pass QFV constraints at $1\sigma$. Green and yellow error bands correspond to
the limits set by the ATLAS experiment Aad  _et al._ (2020c) including run-II
data collected at $\sqrt{s}=13~{}\mathrm{TeV}$.

As one can observe, all points comfortably sit within the bounds set after the
ATLAS run-II search. One may also take note that the calculated cross-sections
in the $\nu$BGL-I model can become fairly small, even for masses close to the
EW scale, although various points with sizable cross-sections are present. For
the $H_{2}$ gluon-fusion production mechanism (Fig. 8a), a point with the
largest signal strength was found to be $0.083~{}\mathrm{pb}$, which can be
probed in future LHC runs. An estimation of the expected number of events at
the high-luminosity (HL) phase of the LHC, with
$\mathcal{L}=3000~{}\mathrm{fb^{-1}}$, leads to $N=\sigma\mathcal{L}\sim
249078$, whereas for run-III luminosities one expects $N\sim 24907$. This
suggests that the low $H_{2}$ mass regime can potentially be tested already
with run-III data in the di-tau final state, as well as in the HL phase. For
the case of gluon-fusion production in association with a $b\bar{b}$ pair
(Fig. 8b), we obtain slightly smaller values for $\sigma\cdot\mathrm{BR}$,
where, the point with the highest cross-section is the same as in the previous
process, with a mass of $132.27~{}\mathrm{GeV}$ and with the signal rate of
$\sigma\cdot\mathrm{BR}=0.0097~{}\mathrm{pb}$, which, at the (run III) HL
phase, corresponds to an expected number of events of around
$N\sim(2910)~{}29100$. As such, one can test this channel at both run-III and
HL-LHC. A better picture is expected for the pseudoscalar production as one
can see in the bottom panels of Fig. 8, where one can find cross sections 1 to
2 orders of magnitude larger when compared to the analogous $H_{2}$
productions and decay channels.

The small signal strengths found for the di-tau final states are, to a large
extent, a result of the small branching fractions generically observed for
$H_{2}\rightarrow\tau^{+}\tau^{-}$ decays. Indeed, as one can see in Tab. 6,
for a selection of five benchmark points, the values of
$\mathrm{BR}(H_{2}\to\tau^{+}\tau^{-})$ are always sub-dominant, with the
exception of benchmarks with sizable cross sections.

Similar conclusions are obtained for the pseudoscalar case as one can see in
the bottom plots of Fig. 8 and in the lower half of Tab. 6. However, as
mentioned above, $A_{2}$ production cross sections with subsequent decay into
taus tend to be greater when compared to $H_{2}$, with the largest value
encountered reaching $0.087~{}\mathrm{pb}$ for a pseudoscalar with mass
$174~{}\mathrm{GeV}$.

The benchmark points BP4 and BP2 in Tabs. 6 and 7 summarize the discussion of
the previous two paragraphs, which were selected as the ones offering the
larger $H_{2}$ and $A_{2}$ production and decay cross sections, respectively.
For various displayed benchmarks, while for the $H_{2}$ field one finds
multiple final states worth exploring, for the $A_{2}$ we note that quark
decays are typically preferred. For the pseudoscalar case, we highlight
benchmarks BP4 and BP5, which result in a final state with two tops and two
charm quarks, respectively. In the later case, it would result in a final
state with two light jets in the detector, whereas for the former case, one
would expect at least two $b$-jets plus two $W$ bosons in the final state,
following the decay of the top quarks.

Besides some points in the low-mass region, the small branching ratios are
generically obtained for the scalar and pseudoscalar decays into a pair of
taus, and as such one might wonder what then the dominant channels are. The
answer is in part given in the last column of Tab. 6, where we have found that
$H_{2}$ mostly decays to $b\bar{b}$ and $WW$ pairs, as well as into
$A_{2}\mathrm{Z^{0}}$ and two Higgs bosons, $H_{1}H_{1}$. The first would
result in a signature with two $b-$jets at the LHC while the second either in
a pair of jets or a pair of leptons or even a jet and a lepton, together with
missing energy. Note that the $H_{2}\rightarrow A_{2}\mathrm{Z^{0}}$ decay
channel results in a interesting signature with at least four jets or two jets
+ two leptons, where we assume that the $A_{2}$ scalar mostly decays into
c-jets, according to the discussion above. A more detailed analysis of this
channel deserves special attention and its collider phenomenology is being the
subject of a separate study currently in preparation. Notice that the
$H_{2}\to A_{2}Z^{0}$ decay channel becomes increasingly important for heavier
CP-even scalars as opposed to lighter $H_{2}$ scenarios where the pure SM
decays dominate. One must also mention that $H_{2}$ can decay in a pair of
Higgs bosons. Such a possibility is rather interesting as it can lead to a
final state with $4b$-jets (via the chain $pp\rightarrow H_{2}\rightarrow
H_{1}H_{1}\rightarrow b\bar{b}b\bar{b}$). Indeed, both ATLAS and CMS have
searched for Higgs boson pairs in the four b-jets final state Aaboud  _et al._
(2019); Sirunyan  _et al._ (2019, 2018b).

$\phi$ | ID | Mass (GeV) | $\mathrm{BR}(\phi\rightarrow\tau^{+}\tau^{-})$ |  $\sigma(gg\rightarrow\phi)\dotproduct\mathrm{BR}$ (pb) |  $\sigma(gg\rightarrow b\bar{b}\phi)\dotproduct\mathrm{BR}$ (pb) | Maximum BR
---|---|---|---|---|---|---
$H_{2}$ |  BP1 BP2 BP3 BP4 BP5 |  $160.21$ $347.99$ $129.26$ $132.27$ $668.49$ |  $4.99\times 10^{-3}$ $6.65\times 10^{-7}$ $0.0127$ $0.0357$ $6.14\times 10^{-6}$ |  $0.007354$ $3.25\times 10^{-8}$ $3.59\times 10^{-4}$ $0.0830$ $7.42\times 10^{-8}$ |  $8.01\times 10^{-4}$ $3.67\times 10^{-9}$ $4.3\times 10^{-5}$ $0.00967$ $1.01\times 10^{-8}$ |  $\mathrm{BR}(W^{+}W^{*})=0.881$ $\mathrm{BR}(H_{1}H_{1})=0.611$ $\mathrm{BR}(W^{+}W^{*})=0.377$ $\mathrm{BR}(b\bar{b})=0.590$ $\mathrm{BR}(A_{2}\mathrm{Z^{0}})=0.75$
$A_{2}$ |  BP1 BP2 BP3 BP4 BP5 |  $194.99$ $173.55$ $1077.21$ $937.61$ $126.74$ |  $0.0336$ $0.0264$ $2.65\times 10^{-5}$ $3.49\times 10^{-5}$ $2.48\times 10^{-3}$ |  $0.0546$ $0.0874$ $1.0\times 10^{-6}$ $6.36\times 10^{-9}$ $3.2\times 10^{-5}$ |  $0.00532$ $0.008249$ $1.62\times 10^{-7}$ $8.66\times 10^{-10}$ $3.0\times 10^{-6}$ |  $\mathrm{BR}(b\bar{b})=0.553$ $\mathrm{BR}(b\bar{b})=0.432$ $\mathrm{BR}(t\bar{t})=0.711$ $\mathrm{BR}(t\bar{t})=0.918$ $\mathrm{BR}(c\bar{c})=0.922$

Table 6: A selection of five benchmark points (BP) that respect all QFV,
electroweak, Higgs and theoretical constraints. All masses are given in GeV
and cross sections in pb. These correspond to $H_{2}$ early discovery
scenarios that maximize the signal strength to di-bosons (BP1,BP2), and to
later $H_{2}$ and $A_{2}$ discovery cases that maximize the di-tau channel
(BP4) and again (BP2), respectively. The (BP3) represents the second lightest
$M_{H_{2}}$ scenario while (BP5) was chosen to represent scenarios with
heavier $H_{2}$ which preferably decay to $A_{2}Z^{0}$. Di-tau final states
are showcased while in the last column we provide the dominant decay mode with
the respective branching ratio.

$\phi$ | ID | Mass (GeV) | $\mathrm{BR}(\phi\rightarrow W^{+}W^{-})$ | $\mathrm{BR}(\phi\rightarrow\mathrm{Z^{0}Z^{0}})$ |  $\sigma(gg\rightarrow\phi\rightarrow W^{+}W^{-})$ (pb) |  $\sigma(gg\rightarrow\phi\rightarrow\mathrm{Z^{0}Z^{0}})$ (pb) | Maximum BR
---|---|---|---|---|---|---|---
$H_{2}$ |  BP1 BP2 BP3 BP4 BP5 |  $160.21$ $347.99$ $129.26$ $132.27$ $668.49$ |  $0.881$ $0.129$ $0.377$ $0.239$ $0.416$ |  $0.0210$ $0.0589$ $0.0449$ $0.0299$ $0.203$ |  $1.0032$ $0.00722$ $3.0\times 10^{-6}$ $0.00366$ $0.00506$ |  $0.00054$ $0.003319$ $1.0\times 10^{-6}$ $0.001426$ $0.00244$ |  $\mathrm{BR}(W^{+}W^{*})=0.881$ $\mathrm{BR}(H_{1}H_{1})=0.611$ $\mathrm{BR}(W^{+}W^{*})=0.377$ $\mathrm{BR}(b\bar{b})=0.590$ $\mathrm{BR}(A_{2}\mathrm{Z^{0}})=0.75$
$H_{3}$ |  BP1 BP2 BP3 BP4 BP5 |  $929.20$ $823.89$ $1093.01$ $1156.49$ $754.64$ |  $0.124$ $0.0922$ $0.166$ $0.0598$ $0.0517$ |  $0.0616$ $0.0455$ $0.0824$ $0.0297$ $0.0254$ |  $0.00341$ $1.36\times 10^{-4}$ $0.002964$ $0.000192$ $8.4\times 10^{-4}$ |  $0.00167$ $6.7\times 10^{-5}$ $0.001468$ $9.4\times 10^{-5}$ $4.1\times 10^{-4}$ |  $\mathrm{BR}(H_{1}H_{1})=0.456$ $\mathrm{BR}(A_{2}A_{2})=0.228$ $\mathrm{BR}(H_{1}H_{2})=0.531$ $\mathrm{BR}(t\bar{t})=0.617$ $\mathrm{BR}(t\bar{t})=0.591$

Table 7: A selection of five benchmark points (BP) that respect all QFV,
electroweak, Higgs and theoretical constraints. All masses are given in GeV
and cross sections in pb. The benchmark points are the same as in Tab. 6. Di-
boson final states are showcased while in the last column we provide the
dominant decay mode with the respective branching ratio.

In Tab. 7 one can also see that, even when not the dominant contributions, the
decay of CP-even scalars into di-boson final states offer larger branching
ratios than those of the tau channel. Recall that, in this article, we do not
assume an exact Higgs alignment limit. Instead, we allow for small deviations
from it, while still guaranteeing a SM-like Higgs candidate to have its mass
and interactions within experimental bounds. In turn, the new CP-even scalars
can couple to both $\mathrm{Z^{0}}$ and $W^{\pm}$ bosons and constraints
arising from LHC direct searches, considering the $\phi\rightarrow W^{+}W^{-}$
and $\phi\rightarrow\mathrm{Z^{0}}\mathrm{Z^{0}}$ channels, with
$\phi=H_{2},H_{3}$, become relevant to our analysis. The results are
summarized in the plots of Fig. 9, for both the $H_{2}$ (left panels) and
$H_{3}$ (right panels). As expected, all points remain bellow the experimental
bounds but with larger signal strengths in comparison with the di-tau channel.
In fact, most of points shown on panels a) and c) can be at the reach of the
LHC run-III or its HL phase. In particular, we have found that for the diboson
channel with $\mathrm{Z^{0}Z^{0}}$, a big portion of the allowed points comes
close to the lower bound set by experiment, with the largest allowed cross-
sections being of the order of $\mathcal{O}(10^{-1})$ and
$\mathcal{O}(10^{-3})~{}\mathrm{pb}$. Therefore, it is possible the upcoming
measurements would be capable of probing the allowed parameter space of the
model. Note that for these two cases, the decay into a pair of W bosons is
indeed the dominant channel. These points are very close to the ATLAS limit
marking them as interesting benchmarks for an early probing, already with run-
III data. It is also worth mentioning that, by virtue of a sizeable branching
ratio to a W boson pair, the signal strength for the heavy CP-even state
$H_{3}$ is comparable to those obtained for the lighter $H_{2}$ state.
Finally, let us comment that, for the case of pseudoscalars, only the tau
channel is considered given that tree-level couplings of $A_{2}$ to vector
bosons are forbidden. Radiative effects where such processes can be induced
are beyond the scope of our discussion in this article.

Figure 9: The signal strength for production of a CP-even scalars via the
gluon-gluon fusion mechanism (a), (b) times its branching ratio to $ZZ$, and
(c), (d) times its branching fraction to $WW$, as a function of the lightest
CP-even mass, $M_{H_{2}}$ (left panels) and heaviest CP-even mass $M_{H_{3}}$
(right panels). The colour code is the same as in Fig. 8 and the exclusion
bounds were taken from Aad  _et al._ (2021b); Aaboud  _et al._ (2018c).

##  VII Conclusions

We have discussed anomaly-free implementations of a NTHDM-BGL model containing
three generations of right-handed neutrinos, imbuing masses to the active
left-handed neutrinos via a Type-I seesaw mechanism. We have specialized our
numerical analysis on a version denoted as $\nu$BGL-I where flavour violation
processes are expected only on the down-type quarks sector. We have
constrained the allowed parameter space of the model upon imposing electroweak
precision, Higgs and flavour physics constraints and scouted the parameter
space in a search for phenomenologically valid regions. We have successfully
assessed the viability of the low mass region and found that even for a number
of scenarios with new scalars at, or at least very close to the EW scale, the
$\nu$BGL-I model remains unconstrained. On the other hand, the majority of the
excluded scenarios met their fate due to a combined effect of the $\Delta
M_{s}$ and $\mathrm{BR}(B_{s}\rightarrow\mu\mu)$ QFV observables, which have
eliminated approximately $90.61\%$ of the sampled points.

For the points that have survived all constraints we have confronted our
results with existing direct searches at the LHC. We have found that new CP-
even scalars are largely favouring final states with a pair of W bosons for
masses not far from the EW scale. However, as their mass grows, the decay
channel $H_{2}\to A_{2}Z^{0}$ becomes dominant, thus a preferable option for
further searches at the LHC run-III or HL phases. Last but not least, while
our results confirm that the di-tau channel is well suited for pseudoscalar
searches, their decay branching ratios to a pair of quarks is in general
largely dominant.

Acknowledgements. Useful discussions and helpful correspondence with Mark
Goodsell are gratefully acknowledged. J.G., F.F.F., A.P.M. and V.V. 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., F.F.F., J.G. and V.V. are
supported by the project PTDC/FIS-PAR/31000/2017. 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 a doctoral program grant with the reference
2021.04527.BD. 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). P.F. is supported by CFTC-UL under FCT
contracts UIDB/00618/2020, UIDP/00618/2020, and by the projects CERN/FIS-
PAR/0002/2017 and CERN/FIS-PAR/0014/2019.

##  Appendix A Anomaly conditions and charge assignment

The anomaly cancellation conditions independent of the neutrino flavour
charges read as

$\displaystyle\begin{split}A_{{}_{\mathrm{U(1)_{Y}}\mathrm{U(1)_{Y}}\mathrm{U(1)^{\prime}}}}&\equiv\sum_{i=1}^{3}\left(X_{q_{i}}+3X_{l_{i}}-8X_{u_{i}}-2X_{d_{i}}-6X_{e_{i}}\right)=0,\\\
A_{{}_{\mathrm{U(1)_{Y}}\mathrm{U(1)^{\prime}}\mathrm{U(1)^{\prime}}}}&\equiv\sum_{i=1}^{3}\left(X_{q_{i}}^{2}-X_{l_{i}}^{2}-2X_{u_{i}}^{2}+X_{d_{i}}^{2}+X_{e_{i}}^{2}\right)=0,\\\
A_{{}_{\mathrm{SU(2)_{L}}\mathrm{SU(2)_{L}}\mathrm{U(1)^{\prime}}}}&\equiv\sum_{i=1}^{3}\left(3X_{q_{i}}+X_{l_{i}}\right)=0,\\\
A_{{}_{\mathrm{SU(3)_{C}}\mathrm{SU(3)_{C}}\mathrm{U(1)^{\prime}}}}&\equiv\sum_{i=1}^{3}\left(2X_{q_{i}}-X_{u_{i}}-X_{d_{i}}\right)=0\,,\end{split}$
(55)

with the $A_{XYZ}$ factor already defined in the main text below Eq. 29.

##  Appendix B Analytical equations for the Lagrangian parameters

$\displaystyle a_{1}=$ $\displaystyle
a_{2}+\frac{\left(M_{A_{2}}^{2}-M_{A_{3}}^{2}\right)\sin\left(2\gamma_{1}\right)}{\sqrt{2}v},$
$\displaystyle\mu_{3}^{2}=$
$\displaystyle\frac{M_{A_{3}}^{2}\sin^{2}\left(\gamma_{1}\right)\left(\cot\left(\gamma_{1}\right)v_{S}-v\sin(\beta)\cos(\beta)\right)-M_{A_{2}}^{2}\cos^{2}\left(\gamma_{1}\right)\left(\tan\left(\gamma_{1}\right)v_{S}+v\sin(\beta)\cos(\beta)\right)}{v}-\sqrt{2}a_{2}v_{S},$
$\displaystyle\mu_{b}^{2}=$
$\displaystyle-\frac{1}{2v_{S}}\big{(}\sqrt{2}a_{2}v^{2}\sin(\beta)\cos(\beta)+M_{A_{3}}^{2}\cos\left(\gamma_{1}\right)\left(\cos\left(\gamma_{1}\right)v_{S}-v\sin(\beta)\cos(\beta)\sin\left(\gamma_{1}\right)\right)$
$\displaystyle+M_{A_{2}}^{2}\sin\left(\gamma_{1}\right)\left(\sin\left(\gamma_{1}\right)v_{S}+v\sin(\beta)\cos(\beta)\cos\left(\gamma_{1}\right)\right)\big{)},$
$\displaystyle\lambda_{1}=$
$\displaystyle\frac{1}{4v^{2}}\big{(}\sec^{2}(\beta)\big{(}-\sin^{2}(\beta)\big{(}\big{(}M_{A_{2}}^{2}-M_{A_{3}}^{2}\big{)}\cos\big(2\gamma_{1}\big{missing})+M_{A_{2}}^{2}+M_{A_{3}}^{2}\big{)}$
$\displaystyle+\sin\big(\alpha_{2}\big{missing})\sin\big(2\alpha_{3}\big{missing})\big{(}M_{H_{2}}^{2}-M_{H_{3}}^{2}\big{)}\sin\big(2\big{(}\alpha_{1}-\beta\big{)}\big{missing})$
$\displaystyle+2\cos^{2}\big{(}\beta-\alpha_{1}\big{)}\big{(}\cos^{2}\big{(}\alpha_{2}\big{)}M_{H_{1}}^{2}+\sin^{2}\big{(}\alpha_{2}\big{)}\big{(}\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}+\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}\big{)}\big{)}$
$\displaystyle+2\sin^{2}\big{(}\beta-\alpha_{1}\big{)}\big{(}\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}+\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}\big{)}\big{)}\big{)},$
$\displaystyle\lambda_{2}=$
$\displaystyle\frac{1}{4v^{2}}\big{(}\csc^{2}(\beta)\big{(}-2\cos^{2}(\beta)\big{(}M_{A_{3}}^{2}\sin^{2}\big{(}\gamma_{1}\big{)}+M_{A_{2}}^{2}\cos^{2}\big{(}\gamma_{1}\big{)}\big{)}$
$\displaystyle+\sin\big(\alpha_{2}\big{missing})\sin\big(2\alpha_{3}\big{missing})\big{(}M_{H_{3}}^{2}-M_{H_{2}}^{2}\big{)}\sin\big(2\big{(}\alpha_{1}-\beta\big{)}\big{missing})+2\cos^{2}\big{(}\alpha_{2}\big{)}M_{H_{1}}^{2}\sin^{2}\big{(}\beta-\alpha_{1}\big{)}$
$\displaystyle+2\sin^{2}\big{(}\alpha_{2}\big{)}\sin^{2}\big{(}\beta-\alpha_{1}\big{)}\big{(}\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}+\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}\big{)}$
$\displaystyle+2\cos^{2}\big{(}\beta-\alpha_{1}\big{)}\big{(}\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}+\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}\big{)}\big{)}\big{)},$
$\displaystyle\lambda_{3}=$
$\displaystyle\frac{1}{2v^{2}}\big{(}\csc(\beta)\sec(\beta)\big{(}2\sqrt{2}\big{(}a_{1}+a_{2}\big{)}v_{S}-\sin\big(2\big{(}\alpha_{1}-\beta\big{)}\big{missing})\big{(}\cos^{2}\big{(}\alpha_{2}\big{)}M_{H_{1}}^{2}$
$\displaystyle+\cos^{2}\big{(}\alpha_{3}\big{)}\big{(}\sin^{2}\big{(}\alpha_{2}\big{)}M_{H_{3}}^{2}-M_{H_{2}}^{2}\big{)}\big{)}+\sin\big(\alpha_{2}\big{missing})\sin\big(2\alpha_{3}\big{missing})\big{(}M_{H_{2}}^{2}-M_{H_{3}}^{2}\big{)}\cos\big(2\big{(}\alpha_{1}-\beta\big{)}\big{missing})+4\mu_{3}^{2}\big{)}$
$\displaystyle+2M_{A_{3}}^{2}\sin^{2}\big{(}\gamma_{1}\big{)}+2M_{A_{2}}^{2}\cos^{2}\big{(}\gamma_{1}\big{)}+4M_{H^{\text{pm}}}^{2}+\sin^{2}\big{(}\alpha_{3}\big{)}\csc(2\beta)\sin\big(2\big{(}\alpha_{1}-\beta\big{)}\big{missing})\big{(}\cos\big(2\alpha_{2}\big{missing})M_{H_{2}}^{2}$
$\displaystyle-M_{H_{2}}^{2}+2M_{H_{3}}^{2}\big{)}\big{)},$
$\displaystyle\lambda_{4}=$
$\displaystyle-\frac{\csc(\beta)\sec(\beta)\left(\sqrt{2}\left(a_{1}+a_{2}\right)v_{S}+\sin(2\beta)M_{H^{\pm}}^{2}+2\mu_{3}^{2}\right)}{v^{2}},$
$\displaystyle\lambda_{1}^{\prime}=$
$\displaystyle\frac{1}{2v_{S}^{3}}\big{(}v\sin(\beta)\cos(\beta)\big{(}\sqrt{2}a_{2}v+\big{(}M_{A_{2}}^{2}-M_{A_{3}}^{2}\big{)}\sin\big(\gamma_{1}\big{missing})\cos\big(\gamma_{1}\big{missing})\big{)}+\sin^{2}\big{(}\alpha_{2}\big{)}M_{H_{1}}^{2}v_{S}$
$\displaystyle+\cos^{2}\big{(}\alpha_{2}\big{)}v_{S}\big{(}\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}+\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}\big{)}\big{)},$
$\displaystyle\lambda_{2}^{\prime}=$
$\displaystyle-\frac{1}{vv_{S}}\big{(}\sec(\beta)\big{(}\sin(\beta)\big{(}\sqrt{2}a_{2}v+\big{(}M_{A_{2}}^{2}-M_{A_{3}}^{2}\big{)}\sin\big(\gamma_{1}\big{missing})\cos\big(\gamma_{1}\big{missing})\big{)}$
$\displaystyle+\cos\big(\alpha_{2}\big{missing})\big{(}\sin\big(\alpha_{2}\big{missing})\cos\big(\beta-\alpha_{1}\big{missing})\big{(}\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}+\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}-M_{H_{1}}^{2}\big{)}$
$\displaystyle-\sin\big(\alpha_{3}\big{missing})\cos\big(\alpha_{3}\big{missing})\big{(}M_{H_{2}}^{2}-M_{H_{3}}^{2}\big{)}\sin\big(\beta-\alpha_{1}\big{missing})\big{)}\big{)}\big{)},$
$\displaystyle\lambda_{3}^{\prime}=$
$\displaystyle-\frac{1}{vv_{S}}\big{(}\csc(\beta)\big{(}\cos(\beta)\big{(}\sqrt{2}a_{2}v+\big{(}M_{A_{2}}^{2}-M_{A_{3}}^{2}\big{)}\sin\big(\gamma_{1}\big{missing})\cos\big(\gamma_{1}\big{missing})\big{)}$
$\displaystyle+\cos\big(\alpha_{2}\big{missing})\big{(}\sin\big(\alpha_{3}\big{missing})\cos\big(\alpha_{3}\big{missing})\big{(}M_{H_{2}}^{2}-M_{H_{3}}^{2}\big{)}\cos\big(\beta-\alpha_{1}\big{missing})$
$\displaystyle-\sin\big(\alpha_{2}\big{missing})\sin\big(\beta-\alpha_{1}\big{missing})\big{(}-\sin^{2}\big{(}\alpha_{3}\big{)}M_{H_{2}}^{2}-\cos^{2}\big{(}\alpha_{3}\big{)}M_{H_{3}}^{2}+M_{H_{1}}^{2}\big{)}\big{)}\big{)}\big{)}$

For completeness, in the case of a model with a neutral singlet (that is, with
a $\mathrm{U(1)}^{\prime}$ charge of zero), all $a_{1}$, $a_{2}$, $a_{3}$ and
$a_{4}$ would be simultaneously allowed, resulting in modified expressions for
the parameters of the scalar potential apart from $\lambda^{\prime}_{3}$.
These are

$\displaystyle a_{4}=$
$\displaystyle\frac{1}{2vv_{S}}\left(\left(-M_{A_{2}}^{2}+M_{A_{3}}^{2}\right)\sin{2\gamma_{1}}+\sqrt{2}v\left(a_{1}-a_{2}\right)+2vv_{S}a_{3}\right)\,,$
$\displaystyle\mu_{3}^{2}=$
$\displaystyle\frac{1}{4v}\Big{(}-v\left(M_{A_{2}}^{2}+M_{A_{3}}^{2}+\left(M_{A_{2}}^{2}-M_{A_{3}}^{2}\right)\cos{2\gamma_{1}}\sin{2\beta}\right)+\left(M_{A_{2}}^{2}-M_{A_{3}}^{2}\right)v_{s}\sin{2\gamma_{1}}$
$\displaystyle+vv_{S}\left(-\sqrt{2}(3a_{1}+a_{2})-4v_{S}a_{3}\right)\Big{)}\,,$
$\displaystyle\mu_{b}^{2}=$
$\displaystyle\frac{1}{4v_{S}}\Big{(}-\left(M_{A_{2}}^{2}+M_{A_{3}}^{2}\right)v_{s}+\left(M_{A_{2}}^{2}-M_{A_{3}}^{2}\right)\left(v_{S}\cos{2\gamma_{1}}+v\sin{2\beta}\sin{2\gamma_{1}}\right)$
$\displaystyle+v^{2}\cos{\beta}\sin{\beta}(-3\sqrt{2}a_{1}+\sqrt{2}a_{2}-8v_{S}a_{3})\Big{)}\,,$
$\displaystyle\lambda_{1}=$
$\displaystyle-\frac{1}{2v^{2}}\Big{(}\sec{\beta}^{2}\left(-M_{H_{1}}^{2}\cos^{2}{\alpha_{2}}\cos^{2}{\delta}+\frac{1}{2}\left(M_{A_{2}}^{2}+M_{A_{3}}^{2}+(M_{A_{2}}^{2}-M_{A_{3}}^{2})\cos{2\gamma_{1}}\right)\right)\sin^{2}{\beta}$
$\displaystyle-
M_{H_{2}}^{2}\left(\cos{\delta}\sin{\alpha_{2}}\sin{\alpha_{3}}+\cos{\alpha_{3}}\sin{\delta}\right)^{2}-M_{H_{3}}^{2}\left(\cos{\alpha_{3}}\cos{\delta}\sin{\alpha_{2}}-\sin{\alpha_{3}}\sin{\delta}\right)^{2}\Big{)}\,,$
$\displaystyle\lambda_{2}=$
$\displaystyle\frac{1}{2v^{2}}\csc^{2}{\beta}\Big{(}M_{H_{3}}^{2}\cos^{2}{\delta}\sin^{2}{\alpha_{3}}-\cos^{2}{\beta}\left(M_{A_{2}}^{2}\cos^{2}{\gamma_{1}}+M_{A_{3}}^{2}\sin^{2}{\gamma_{1}}\right)$
$\displaystyle+\left(M_{H_{1}}^{2}\cos^{2}{\alpha_{2}}+M_{H_{2}}^{2}\sin^{2}{\alpha_{2}}\sin^{2}{\alpha_{3}}\right)\sin^{2}{\delta}+\cos^{2}{\alpha}\left(M_{H_{2}}^{2}\cos^{2}{\delta}+M_{H_{3}}\sin^{2}{a2}\sin^{2}{\delta}\right)$
$\displaystyle+\left(-M_{H_{2}}+M_{H_{3}}^{2}\right)\cos{\alpha_{3}}\sin{\alpha_{2}}\sin{\alpha_{3}}\sin{2\delta}\Big{)}\,,$
$\displaystyle\lambda_{3}=$
$\displaystyle\frac{1}{8v^{2}}\Big{(}-4\left(M_{A_{2}}^{2}+M_{A_{3}}^{2}-4M_{H^{\pm}}^{2}\right)+4\left(-M_{A_{2}}^{2}+M_{A_{3}}^{2}\right)\cos{2\gamma_{1}}$
$\displaystyle+\csc{\beta}\sec{\beta}\big{(}4\left(M_{H_{2}}^{2}-M_{H_{3}}^{2}\right)\cos{2\delta}\sin{\alpha_{2}}\sin{2\alpha_{3}}+(2\left(-2M_{H_{1}}^{2}+M_{H_{2}}^{2}+M_{H_{3}}^{2}\right)\cos^{2}{\alpha_{2}}$
$\displaystyle-\left(M_{H_{2}}^{2}-M_{H_{3}}^{2}\right)\left(-3+\cos{2\alpha_{2}}\right)\cos{2\alpha_{3}})\sin{2\delta}\big{)}\Big{)}\,,$
$\displaystyle\lambda_{4}=$
$\displaystyle\frac{1}{v^{2}}\left(M_{A_{2}}^{2}+M_{A_{3}}^{2}-2M_{H^{\pm}}^{2}+\left(M_{A_{2}}^{2}-M_{A_{3}}^{2}\right)\cos{2\gamma_{1}}\right)\,,$
$\displaystyle\lambda_{1}^{\prime}=$
$\displaystyle\frac{1}{4v_{S}^{3}}\left(2v_{s}\left(M_{H_{1}}^{2}\sin^{2}{\alpha_{2}}+\cos^{2}{\alpha_{2}}\left(M_{H_{3}}^{2}\cos^{2}{\alpha_{3}}+M_{H_{2}}^{2}\sin^{2}{\alpha_{3}}\right)\right)+\sqrt{2}v^{2}\cos{\beta}Sin{\beta}\left(a_{1}+a_{2}\right)\right)\,,$
$\displaystyle\lambda_{2}^{\prime}=$
$\displaystyle\frac{1}{vv_{S}}\Big{(}\cos{\alpha_{2}}\sec{\beta}(\cos{\delta}\sin{\alpha_{2}}\left(M_{H_{1}}^{2}-M_{H_{3}}^{2}\cos^{2}{\alpha_{3}}-M_{H_{2}}^{2}\sin^{2}{\alpha_{3}}\right)$
$\displaystyle+\left(-M_{H_{2}}^{2}+M_{H_{3}}^{2}\right)\cos{\alpha_{3}}\sin{\alpha_{3}}\sin{\delta})+\left(M_{A_{2}}^{2}-M_{A_{3}}^{2}\right)\cos{\gamma
1}\sin{\gamma 1}\tan{\beta}$
$\displaystyle-v\left(\sqrt{2}a_{1}+2v_{S}a_{3}\right)\tan{\beta}\Big{)}\,,$

and where the minimization conditions only modify $\mu_{S}^{2}$ and
$\hat{\mu}$ as

$\mu_{S}^{2}=\quad-\frac{\sqrt{2}v_{1}v_{2}\left(a_{1}+a_{2}\right)+v_{S}\left(v_{1}\left(2v_{2}\left(a_{3}+a_{4}\right)+v_{1}\lambda_{2}^{\prime}\right)+v_{2}^{2}\lambda^{\prime}_{3}+2\left(v_{S}^{2}\lambda^{\prime}_{1}+\mu_{b}^{2}\right)\right)}{2v_{S}}$

and

$\hat{\mu}=v_{S}\left[\sqrt{2}\left(a_{1}+a_{2}\right)+v_{S}\left(a_{3}+a_{4}\right)\right]+v_{1}v_{2}\left(\lambda_{3}+\lambda_{4}\right)+2\mu_{3}^{2}\,.$

## References

  * Aad  _et al._ (2012) G. Aad _et al._ (ATLAS), Phys. Lett. B 716, 1 (2012), arXiv:1207.7214 [hep-ex] .
  * Chatrchyan  _et al._ (2012) S. Chatrchyan _et al._ (CMS), Phys. Lett. B 716, 30 (2012), arXiv:1207.7235 [hep-ex] .
  * Lee (1973) T. D. Lee, Phys. Rev. D 8, 1226 (1973).
  * Cárcamo Hernández  _et al._ (2021) A. E. Cárcamo Hernández, C. Espinoza, J. Carlos Gómez-Izquierdo, and M. Mondragón, (2021), arXiv:2104.02730 [hep-ph] .
  * Arcadi  _et al._ (2022) G. Arcadi, A. S. de Jesus, T. B. de Melo, F. S. Queiroz, and Y. S. Villamizar, Nucl. Phys. B 982, 115882 (2022), arXiv:2104.04456 [hep-ph] .
  * Bell  _et al._ (2017) N. F. Bell, G. Busoni, and I. W. Sanderson, JCAP 03, 015 (2017), arXiv:1612.03475 [hep-ph] .
  * Dorsch  _et al._ (2017) G. C. Dorsch, S. J. Huber, K. Mimasu, and J. M. No, JHEP 12, 086 (2017), arXiv:1705.09186 [hep-ph] .
  * Camargo  _et al._ (2019) D. A. Camargo, M. D. Campos, T. B. de Melo, and F. S. Queiroz, Phys. Lett. B795, 319 (2019), arXiv:1901.05476 [hep-ph] .
  * Ivanov (2017) I. P. Ivanov, Prog. Part. Nucl. Phys. 95, 160 (2017), arXiv:1702.03776 [hep-ph] .
  * Branco  _et al._ (2012) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, Phys. Rept. 516, 1 (2012), arXiv:1106.0034 [hep-ph] .
  * Grzadkowski  _et al._ (2014) B. Grzadkowski, O. M. Ogreid, and P. Osland, JHEP 01, 105 (2014), arXiv:1309.6229 [hep-ph] .
  * Boto  _et al._ (2020) R. Boto, T. V. Fernandes, H. E. Haber, J. C. Romão, and J. P. Silva, Phys. Rev. D 101, 055023 (2020), arXiv:2001.01430 [hep-ph] .
  * Grzadkowski  _et al._ (2016) B. Grzadkowski, O. M. Ogreid, and P. Osland, Phys. Rev. D 94, 115002 (2016), arXiv:1609.04764 [hep-ph] .
  * Sokolowska  _et al._ (2008) D. Sokolowska, K. A. Kanishev, and M. Krawczyk, PoS CHARGED2008, 016 (2008), arXiv:0812.0296 [hep-ph] .
  * Crivellin  _et al._ (2013) A. Crivellin, A. Kokulu, and C. Greub, Phys. Rev. D87, 094031 (2013), arXiv:1303.5877 [hep-ph] .
  * Kim  _et al._ (2015) C. S. Kim, Y. W. Yoon, and X.-B. Yuan, JHEP 12, 038 (2015), arXiv:1509.00491 [hep-ph] .
  * Branco  _et al._ (1996) G. C. Branco, W. Grimus, and L. Lavoura, Phys. Lett. B380, 119 (1996), arXiv:hep-ph/9601383 [hep-ph] .
  * Ordell  _et al._ (2020) A. Ordell, R. Pasechnik, and H. Serôdio, Phys. Rev. D 102, 035016 (2020), arXiv:2006.08676 [hep-ph] .
  * Nottensteiner  _et al._ (2019) F. Nottensteiner, A. Ordell, R. Pasechnik, and H. Serôdio, Phys. Rev. D100, 115038 (2019), arXiv:1909.05548 [hep-ph] .
  * Campos  _et al._ (2017) M. D. Campos, D. Cogollo, M. Lindner, T. Melo, F. S. Queiroz, and W. Rodejohann, JHEP 08, 092 (2017), arXiv:1705.05388 [hep-ph] .
  * Aad  _et al._ (2020a) G. Aad _et al._ (ATLAS), JHEP 10, 061 (2020a), arXiv:2005.05138 [hep-ex] .
  * Aaboud  _et al._ (2018a) M. Aaboud _et al._ (ATLAS), JHEP 10, 180 (2018a), arXiv:1807.11471 [hep-ex] .
  * Aad  _et al._ (2021a) G. Aad _et al._ (ATLAS), JHEP 11, 209 (2021a), arXiv:2108.13391 [hep-ex] .
  * Sirunyan  _et al._ (2020a) A. M. Sirunyan _et al._ (CMS), JHEP 05, 033 (2020a), arXiv:1911.03947 [hep-ex] .
  * Sirunyan  _et al._ (2021) A. M. Sirunyan _et al._ (CMS), JHEP 07, 208 (2021), arXiv:2103.02708 [hep-ex] .
  * Sirunyan  _et al._ (2018a) A. M. Sirunyan _et al._ (CMS), JHEP 04, 073 (2018a), arXiv:1802.01122 [hep-ex] .
  * Staub (2014) F. Staub, Comput. Phys. Commun. 185, 1773 (2014), arXiv:1309.7223 [hep-ph] .
  * Porod and Staub (2012) W. Porod and F. Staub, Comput. Phys. Commun. 183, 2458 (2012), arXiv:1104.1573 [hep-ph] .
  * Bechtle  _et al._ (2015) P. Bechtle, S. Heinemeyer, O. Stal, T. Stefaniak, and G. Weiglein, Eur. Phys. J. C 75, 421 (2015), arXiv:1507.06706 [hep-ph] .
  * Bechtle  _et al._ (2014) P. Bechtle, S. Heinemeyer, O. Stål, T. Stefaniak, and G. Weiglein, Eur. Phys. J. C 74, 2711 (2014), arXiv:1305.1933 [hep-ph] .
  * Straub (2018) D. M. Straub, (2018), arXiv:1810.08132 [hep-ph] .
  * Ordell (2020) A. Ordell,  _Phenomenology in multi-scalar extensions of the Standard Model_ , Ph.D. thesis (2020), Lund University.
  * Botella  _et al._ (2011) F. J. Botella, G. C. Branco, M. Nebot, and M. N. Rebelo, JHEP 10, 037 (2011), arXiv:1102.0520 [hep-ph] .
  * Das  _et al._ (2021) D. Das, P. M. Ferreira, A. P. Morais, I. Padilla-Gay, R. Pasechnik, and J. P. Rodrigues, JHEP 11, 079 (2021), arXiv:2106.06425 [hep-ph] .
  * Zyla  _et al._ (2020a) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020a).
  * Aaboud  _et al._ (2018b) M. Aaboud _et al._ (ATLAS), Phys. Rev. D 98, 052008 (2018b), arXiv:1808.02380 [hep-ex] .
  * Khachatryan  _et al._ (2016a) V. Khachatryan _et al._ (CMS), Phys. Lett. B 755, 217 (2016a), arXiv:1510.01181 [hep-ex] .
  * Khachatryan  _et al._ (2016b) V. Khachatryan _et al._ (CMS), Phys. Lett. B 759, 369 (2016b), arXiv:1603.02991 [hep-ex] .
  * Aad  _et al._ (2020b) G. Aad _et al._ (ATLAS), Phys. Rev. D 102, 032004 (2020b), arXiv:1907.02749 [hep-ex] .
  * Aad  _et al._ (2021b) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 81, 332 (2021b), arXiv:2009.14791 [hep-ex] .
  * Sirunyan  _et al._ (2020b) A. M. Sirunyan _et al._ (CMS), JHEP 03, 055 (2020b), arXiv:1911.03781 [hep-ex] .
  * Zyla  _et al._ (2020b) P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020b).
  * Aad  _et al._ (2020c) G. Aad _et al._ (ATLAS), Phys. Rev. Lett. 125, 051801 (2020c), arXiv:2002.12223 [hep-ex] .
  * Aaboud  _et al._ (2018c) M. Aaboud _et al._ (ATLAS), Eur. Phys. J. C 78, 24 (2018c), arXiv:1710.01123 [hep-ex] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Aaboud  _et al._ (2019) M. Aaboud _et al._ (ATLAS), JHEP 01, 030 (2019), arXiv:1804.06174 [hep-ex] .
  * Sirunyan  _et al._ (2019) A. M. Sirunyan _et al._ (CMS), JHEP 01, 040 (2019), arXiv:1808.01473 [hep-ex] .
  * Sirunyan  _et al._ (2018b) A. M. Sirunyan _et al._ (CMS), Phys. Lett. B 781, 244 (2018b), arXiv:1710.04960 [hep-ex] .