# Global analysis and LHC study of a vector-like extension of the Standard
Model with extra scalars

A. E. Cárcamo Hernández<EMAIL_ADDRESS>Universidad Técnica Federico
Santa María,
Casilla 110-V, Valparaíso, Chile  Centro Científico-Tecnológico de Valparaíso,
Casilla 110-V, Valparaíso, Chile  Millennium Institute for Subatomic Physics
at High-Energy Frontier (SAPHIR), Fernández Concha 700, Santiago, Chile
Kamila Kowalska<EMAIL_ADDRESS>National Centre for Nuclear
Research,
Pasteura 7, 02-093 Warsaw, Poland  Huchan Lee<EMAIL_ADDRESS>National
Centre for Nuclear Research,
Pasteura 7, 02-093 Warsaw, Poland  Daniele Rizzo<EMAIL_ADDRESS>National Centre for Nuclear Research,
Pasteura 7, 02-093 Warsaw, Poland

###### Abstract

We perform a global analysis of a vector-like extension of the Standard Model,
which also features additional doublet and singlet scalars. The usual Yukawa
interactions are forbidden in this setup by an extra U(1) global symmetry and
the masses of the second and third family quarks and leptons are generated via
the mixing with the vector-like sector. We identify three best-fit benchmark
scenarios which satisfy the constraints imposed by the stability of the scalar
potential, the perturbativity of the coupling constants, the measurement of
the muon anomalous magnetic moment and the non-observation of the flavor
violating tau decays. We show that dominant contributions to the muon $(g-2)$
originate in this model from the charged Higgs/neutral lepton one-loop
diagrams, thus correcting an inaccurate statement than can be found in the
literature. We also perform a detailed LHC analysis of the benchmark
scenarios. We investigate the experimental constraints stemming from direct
searches for vector-like quarks, vector-like leptons and exotic scalars. While
we show that the model is not currently tested by any collider experiment, we
point out that decays of a heavy Higgs boson into two tau leptons may offer a
smoking gun signature for the model verification in upcoming runs at the LHC.

###### Contents

  1. I Introduction
  2. II Generation of fermion masses and mixing
    1. II.1 Hierarchy of masses
    2. II.2 CKM mixing matrix
  3. III Scalar potential constraints
    1. III.1 Scalar masses in the alignment limit
    2. III.2 Bounded-from-below limits
    3. III.3 Vacuum stability
  4. IV Flavor physics constraints
    1. IV.1 Muon anomalous magnetic moment
    2. IV.2 Lepton flavor violating decays
    3. IV.3 CKM anomaly
  5. V Perturbativity constraints
  6. VI Numerical analysis and benchmark scenarios
    1. VI.1 Scanning methodology
    2. VI.2 Benchmark scenarios
  7. VII LHC study of the benchmark scenarios
    1. VII.1 Vector-like quarks
    2. VII.2 Vector-like leptons
    3. VII.3 Exotic scalars
  8. VIII Conclusions
  9. A Fermion mass matrices
    1. A.1 Charged leptons
    2. A.2 Up-type quarks
    3. A.3 Down-type quarks
    4. A.4 Neutrino sector
  10. B Scalar mass matrices
  11. C Derivation of the bounded-from-below conditions
  12. D Renormalization group equations

##  I Introduction

The origin of the flavor structure of the Standard Model (SM), i.e. the
observed hierarchy between fermion masses and mixing angles of the Cabibbo-
Kobayashi-Maskawa (CKM) matrix, is one of the greatest mysteries of particle
physics that still lacks a convincing and commonly accepted explanation. A
number of New Physics (NP) ideas have been put forward in recent decades to
address the flavor puzzle, among which the Froggatt-Nielsen (FN) mechanism
Froggatt and Nielsen (1979) and extra dimensions Arkani-Hamed and Schmaltz
(2000); Gherghetta and Pomarol (2000); Huber and Shafi (2001) are those that
admittedly received the most attention and applications. The underlying
concept is to introduce a new quantity that in some sense would be “larger”
than the electroweak symmetry breaking (EWSB) scale. This could be the vacuum
expectation value (vev) of a flavon field, or a distance of a fermion field
from the infra-red brane. Such a hierarchy of scales can be then translated
into a hierarchy of masses and mixing angles of the SM quarks and leptons. A
similar idea gave rise to the famous seesaw mechanism of neutrino mass
generation Minkowski (1977); Gell-Mann et al. (1979); Yanagida (1979); Glashow
(1980); Mohapatra and Senjanovic (1981); Schechter and Valle (1982, 1980),
where tiny values of the SM neutrino masses arise as a result of suppression
of the EWSB scale by a very large Majorana mass.

In Ref. King (2018) a FN-inspired model was proposed to explain the observed
masses and mixing patterns of the SM fermions. The SM Yukawa interactions are
forbidden in this setup by an extra abelian symmetry U$(1)_{X}$, which could
be either global or local. The particle content of the model corresponds to
the two-Higgs-doublet model (2HDM) extended by a full family of vector-like
(VL) fermions charged under U$(1)_{X}$, and one U$(1)_{X}$-breaking singlet
scalar which plays the role of a FN flavon. The large third-family Yukawa
couplings are then effectively generated via mixing of the SM quarks and
leptons with the SU$(2)_{L}$ doublet VL fermions, while the Yukawa couplings
of the second family emerge from a seesaw-like construction, mediated by the
heavy VL SU$(2)_{L}$ singlets.

The rich structure of the model introduced in Ref. King (2018) makes it a
perfect framework for providing a combined explanation both for the flavor
pattern of the SM and for the miscellaneous anomalies which emerged in recent
years in collider experiments. In this context, lepton-flavor violating
anomalies in the rare semi-leptonic decays of the $B$ mesons were analyzed in
Refs. King (2018); Lee and Cárcamo Hernández (2022), $Z$-mediated Flavour
Changing Neutral Currents in Ref. Cárcamo Hernández et al. (2022), and a
deviation from the SM prediction in the measured value of the anomalous
magnetic moment of muon in Refs. Cárcamo Hernández et al. (2020); Hernández et
al. (2021); Lee and Cárcamo Hernández (2022). In the latter study, in which
the extra U$(1)_{X}$ symmetry was assumed to be global, five benchmark points
were identified that could account for the muon $(g-2)$ anomaly and, at the
same time, give rise to the mass and mixing patterns of the SM fermions. The
scenarios pinpointed in Ref. Hernández et al. (2021) were characterized by
relatively low ($\sim 200$ GeV) masses of the VL lepton doublets and large
($\sim 10$) quartic couplings of the scalar potential, which may indicate a
loss of perturbativity at scales very close to the typical scale set by the
masses of the NP particles in the analyzed model.

In this study, we reassess the findings of Ref. Hernández et al. (2021)
improving and extending its analysis in several different directions. Firstly,
we thoroughly discuss the impact of the most recent bounds from direct NP
searches at the Large Hadron Collider (LHC) on the allowed parameter space of
the model, a topic which was not addressed in detail in Refs. Cárcamo
Hernández et al. (2020, 2022); Hernández et al. (2021); Lee and Cárcamo
Hernández (2022). While we show that the model is not currently tested by any
collider experiment, we point out that decays of a heavy Higgs boson into two
tau leptons may offer a smoking gun signature for the detection of the model
in the upcoming runs of the LHC.

Secondly, we demonstrate that the quartic and Yukawa couplings of the model
are subject to strong constraints from their renormalization group (RG)
running. In Ref. Hernández et al. (2021) it was required that all the
dimensionless parameters of the lagrangian remain perturbative (in a loose
sense of being smaller than $\sqrt{4\pi}$ for the gauge/Yukawa and smaller
than $4\pi$ for the scalar potential couplings) at the characteristic energy
scale of the model. We argue that such a simplistic implementation of the
perturbativity bounds should be taken with a grain of salt. The breakdown of
perturbativity usually calls for an extension of the theoretical setup by
extra degrees of freedom in order to cure pathological behavior of the running
couplings, or/and for an inclusion of non-perturbative effects (like bound-
state formation). If any of those arose at the scale specific to the original
NP model, they would most likely affect its phenomenological predictions.
Therefore, it is more correct to apply the perturbativity bounds to the
running couplings evaluated at an energy scale which is high enough that the
phenomenology of the specific NP model can be trusted. Once this improvement
had been implemented in our study, we discovered that all the benchmark points
found previously in Ref. Hernández et al. (2021) were disfavored.

Last but not least, we refine the derivation of the stability conditions for
the scalar potential which in Refs. Hernández et al. (2021); Lee and Cárcamo
Hernández (2022) was simplified to the 2HDM case by integrating out the
singlet flavon field. In the current work we derive all the relevant stability
conditions in the full three-scalar setup, obtaining additional constraints on
the quartic couplings.

With all the improvements in place, we identify three benchmark scenarios that
satisfy our theoretical and experimental requirements. While these best-fit
points emerge from a random numerical scan, they present features that are
generic for the model in study. Most importantly, we point out that a charged
Higgs/heavy neutrino loop is a dominant contribution to the muon $(g-2)$
anomaly. This results from the fact that the competing neutral scalar/heavy
charged lepton contributions are governed by the same Yukawa coupling that
determines the tree-level muon mass and is thus required to be small. Once
more, this finding is qualitatively different from the conclusions obtained in
Refs. Hernández et al. (2021); Lee and Cárcamo Hernández (2022), where only
the charged lepton loops were considered.

The structure of the paper is the following. In Sec. II we briefly review the
field content of the model. We also show how the SM fermion masses and the CKM
matrix are generated in this framework. Sec. III is dedicated to the scalar
sector of the theory. Tree-level scalar masses in the alignment limit are
presented, as well as three-field potential stability conditions. Experimental
constraints from the flavor physics observables (muon $(g-2)$, rare $\tau$
decays, CKM anomaly) are examined in Sec. IV. In Sec. V we discuss the RG flow
of the model couplings and we derive the corresponding perturbativity bounds.
Sec. VI comprises the numerical analysis of the model. We discuss the setup of
our numerical scan and we identify three benchmark scenarios that satisfy all
the theoretical and phenomenological constraints. In Sec. VII we present a
detailed analysis of the LHC searches that may test the parameter space of the
model. We summarize our findings in Sec. VIII. Appendices feature,
respectively, explicit forms of the fermion (Appendix A) and scalar (Appendix
B) mass matrices, derivation of the bounded-from-below constraints (Appendix
C), and the RG equations (Appendix D).

##  II Generation of fermion masses and mixing

We begin our study by reviewing the structure and the main properties of the
model introduced in Ref. King (2018). In the following, we focus mostly on
these features of the model which play a pivotal role in the subsequent
phenomenological analysis. Technical details of the model, including the
analytical diagonalization of the fermion mass matrices and the derivation of
the interaction vertices in the mass basis, can be found in Refs. King (2018);
Cárcamo Hernández et al. (2022); Hernández et al. (2021); Lee and Cárcamo
Hernández (2022).

The particle content of the model is summarized in Table 1. The SM fermion
sector, collectively denoted as $\psi_{i}$
($\psi_{i}=Q_{iL},\,u_{iR},\,d_{iR},\,L_{iL},\,e_{iR}$ and $i=1,2,3$ stands
for a generation index) is extended by one full family of VL fermions,
indicated collectively as $(\psi_{4},\,\widetilde{\psi}_{4})$. We adopt the
convention of using the left-chiral two-component Weyl spinors, therefore the
subscripts $L,R$ indicate the names of the fermions, not the chiralities. The
scalar sector contains, besides the usual SU(2)L Higgs doublet dubbed as
$H_{u}$, an extra scalar doublet $H_{d}$ and a scalar singlet $\phi$. Note
that all the NP particles and the Higgs doublet $H_{u}$ are charged under an
extra global gauge symmetry U$(1)_{X}$, while the SM fermions are U$(1)_{X}$
singlets. As a result, the ordinary SM Yukawa interactions are forbidden.

Field | $Q_{iL}$ | $u_{iR}$ | $d_{iR}$ | $L_{iL}$ | $e_{iR}$ | $Q_{4L}$ | $u_{4R}$ | $d_{4R}$ | $L_{4L}$ | $e_{4R}$ | $\nu_{4R}$ | $\widetilde{Q}_{4R}$ | $\widetilde{u}_{4L}$ | $\widetilde{d}_{4L}$ | $\widetilde{L}_{4R}$ | $\widetilde{e}_{4L}$ | $\widetilde{\nu}_{4L}$ | $\phi$ | $H_{u}$ | $H_{d}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
SU$(3)_{C}$ | $\mathbf{3}$ | $\mathbf{\bar{3}}$ | $\mathbf{\bar{3}}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{3}$ | $\mathbf{\bar{3}}$ | $\mathbf{\bar{3}}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{\bar{3}}$ | $\mathbf{3}$ | $\mathbf{3}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$
SU$(2)_{L}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{1}$ | $\mathbf{2}$ | $\mathbf{2}$
U$(1)_{Y}$ | $\frac{1}{6}$ | $-\frac{2}{3}$ | $\frac{1}{3}$ | $-\frac{1}{2}$ | $1$ | $\frac{1}{6}$ | $-\frac{2}{3}$ | $\frac{1}{3}$ | $-\frac{1}{2}$ | $1$ | $0$ | $-\frac{1}{6}$ | $\frac{2}{3}$ | $-\frac{1}{3}$ | $\frac{1}{2}$ | $-1$ | $0$ | $0$ | $\frac{1}{2}$ | $-\frac{1}{2}$
U$(1)_{X}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$ | $-1$ | $-1$ | $-1$ | $-1$ | $-1$ | $-1$ | $1$ | $-1$ | $-1$
Table 1: The particle content of the NP model considered in this study.

All the renormalizable Yukawa interactions between the SM and NP fermions
which are allowed by the extended gauge symmetry can be schematically written
as:

$\mathcal{L}_{\mathop{\mathrm{r}en}}^{\mathop{\mathrm{Y}ukawa}}=y_{i4}^{\psi}\psi_{iL}H\psi_{4R}+y_{4j}^{\psi}\psi_{4L}H\psi_{jR}+x_{i4}^{\psi}\psi_{iL}\phi\,\widetilde{\psi}_{4R}+x_{4j}^{\psi}\widetilde{\psi}_{4L}\phi\,\psi_{jR}\\\
+M_{4}^{\psi}\,\psi_{4L}\widetilde{\psi}_{4R}+M_{4}^{\widetilde{\psi}}\,\widetilde{\psi}_{4L}\psi_{4R}+\mathop{\mathrm{h}.c.},$
(1)

where $H$ is either $H_{u}$ or $H_{d}$ and $M_{4}^{\psi}$
($M_{4}^{\widetilde{\psi}}$) denotes the VL doublet (singlet) mass parameter.
Note that with the U$(1)_{X}$ charges given in Table 1 the scalar $H_{u}$ only
couples to the up-type quarks, while $H_{d}$ to the down-type quarks and
charged leptons, reminiscent of the 2HDM Type-II model.

###  II.1 Hierarchy of masses

Once the neutral components of the scalar fields develop their vevs, the
$5\times 5$ fermions mass matrices are generated. Since their upper $3\times
3$ blocks contain only zeros (we recall that the SM Yukawa couplings are
forbidden by the U$(1)_{X}$ symmetry), one has the freedom to rotate the first
three families. It can easily be shown King (2018) that this allows one to
choose a flavor basis in which the fermion mass matrices acquire the following
form:

$\mathcal{M}_{\psi}=\left(\begin{array}[]{c|ccccc}&\psi_{1R}&\psi_{2R}&\psi_{3R}&\psi_{4R}&\widetilde{\psi}_{4R}\\\
\hline\cr\psi_{1L}&0&0&0&(y^{\psi}_{14}\langle H^{0}\rangle)&0\\\
\psi_{2L}&0&0&0&y^{\psi}_{24}\langle H^{0}\rangle&0\\\
\psi_{3L}&0&0&0&y^{\psi}_{34}\langle
H^{0}\rangle&x^{\psi}_{34}\langle\phi\rangle\\\
\psi_{4L}&0&0&y^{\psi}_{43}\langle H^{0}\rangle&0&M_{4}^{\psi}\\\
\widetilde{\psi}_{4L}&0&x^{\psi}_{42}\langle\phi\rangle&x^{\psi}_{43}\langle\phi\rangle&M_{4}^{\widetilde{\psi}}&0\\\
\end{array}\right)\,.$ (2)

In the above, the term in parentheses assumes a non-zero value in the mass
matrix of the down-type quarks, while it is zero for the up-type quarks and
charged leptons. The exact forms of the matrices $\mathcal{M}_{u}$,
$\mathcal{M}_{d}$ and $\mathcal{M}_{e}$ are presented in Appendix A.

In order to calculate the masses of the physical quarks and leptons, the
$5\times 5$ matrices of Eq. (2) need to be diagonalized. Due to a large number
of free parameters in the Yukawa sector one may expect that the resulting
functional dependence of the eigenvalues of $\mathcal{M}_{\psi}$ on the
couplings $y^{\psi}_{i4},\,y^{\psi}_{43},\,x^{\psi}_{4i},\,x^{\psi}_{34}$ and
the masses $M_{4}^{\psi(\widetilde{\psi})}$ is highly nontrivial. It turns
out, however, that it is not necessarily the case and that simplified
expressions for the fermion masses can be derived. Denoting the scalar vevs as

$\langle H_{u}^{0}\rangle=v_{u}/\sqrt{2},\qquad\langle
H_{d}^{0}\rangle=v_{d}/\sqrt{2},\qquad\langle\phi\rangle=v_{\phi}/\sqrt{2}$
(3)

and defining $\tan\beta=v_{u}/v_{d}$, the masses of the third and second
family quarks and leptons are approximately given by (see also Ref. King
(2018) for a related derivation)

$\displaystyle m_{t}$ $\displaystyle\approx$
$\displaystyle\frac{1}{\sqrt{2}}\frac{y^{u}_{43}x^{Q}_{34}v_{\phi}v_{u}}{\sqrt{(x^{Q}_{34}v_{\phi})^{2}+2(M_{4}^{Q})^{2}}},\qquad
m_{c}\approx\frac{y^{u}_{24}x^{u}_{42}v_{\phi}v_{u}}{2\,M_{4}^{u}}$ (4)
$\displaystyle m_{b}$ $\displaystyle\approx$
$\displaystyle\frac{1}{\sqrt{2}}\frac{y^{d}_{43}x^{Q}_{34}v_{\phi}v_{d}}{\sqrt{(x^{Q}_{34}v_{\phi})^{2}+2(M_{4}^{Q})^{2}}},\qquad
m_{s}\approx\frac{y^{d}_{24}x^{d}_{42}v_{\phi}v_{d}}{2\,M_{4}^{d}}$ (5)
$\displaystyle m_{\tau}$ $\displaystyle\approx$
$\displaystyle\frac{1}{\sqrt{2}}\frac{y^{e}_{43}x^{L}_{34}v_{\phi}v_{d}}{\sqrt{(x^{L}_{34}v_{\phi})^{2}+2(M_{4}^{L})^{2}}},\qquad
m_{\mu}\approx\frac{y^{e}_{24}x^{e}_{42}v_{\phi}v_{d}}{2\,M_{4}^{e}}\,.$ (6)

While Eqs. (4)-(6) allow determination of the SM fermion masses with an
accuracy within a factor of $2-3$ only, they can be used to gain intuition of
which NP Yukawa couplings play a dominant role in establishing the correct
masses of particular fermions. For example, large $x^{Q}_{34}$ and
$y^{u}_{43}$ are expected to fit $m_{t}$, while $y^{u}_{24}\ll 1$ or
$x^{u}_{42}\ll 1$ would be required to suppress the charm mass. Similarly,
large $y^{d}_{43}$ is needed to generate $m_{b}=4.18$ GeV. Additionally, in
order to obtain the correct value of the top quark mass, the singlet scalar
vev $v_{\phi}$ should be of the same order as the VL mass parameter
$M_{4}^{Q}$. Note, however, that in our phenomenological analysis we always
perform the numerical diagonalization of the mass matrices (67), (70) and
(73).

One important observation which can be deduced from Eqs. (4) and (5) is that
the ratio of the top and bottom masses, $m_{t}/m_{b}\approx 34$, puts relevant
constraints on the allowed parameter space of the model. In fact, we have

$\frac{m_{t}}{m_{b}}\approx\frac{y_{43}^{u}}{y_{43}^{d}}\tan\beta\,.$ (7)

The relation (7) leads to two distinct classes of solutions. In the first one,
with both the Yukawa couplings of order one, $\tan\beta\sim\mathcal{O}(10)$ is
required. In the other one, with $\tan\beta\sim\mathcal{O}(1)$, a large
hierarchy between the up and down sector couplings, $y^{d}_{43}\ll
y^{u}_{43}$, must be imposed.

The masses of VL fermions are given, to a very good approximation, by the
corresponding VL mass parameters with small contributions stemming from their
mixing with the second and the third family,

$\displaystyle M_{U_{1}}$ $\displaystyle\approx$
$\displaystyle\sqrt{(M^{Q}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{Q}_{34})^{2}-\frac{(M_{4}^{Q}y^{u}_{43}v_{u})^{2}}{(x^{Q}_{34}v_{\phi})^{2}+2(M_{4}^{Q})^{2}}}$
(8) $\displaystyle M_{U_{2}}$ $\displaystyle\approx$
$\displaystyle\sqrt{(M^{u}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{u}_{43})^{2}+\frac{1}{2}(v_{\phi}x^{u}_{42})^{2}+\frac{2(M^{u}_{4}y^{u}_{43}v_{u})^{2}}{2(M^{u}_{4})^{2}+(v_{\phi}x^{u}_{43})^{2}+(v_{\phi}x^{u}_{42})^{2}}}$
(9) $\displaystyle M_{D_{1}}$ $\displaystyle\approx$
$\displaystyle\sqrt{(M^{Q}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{Q}_{34})^{2}},\qquad
M_{D_{2}}=\sqrt{(M^{d}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{d}_{43})^{2}+\frac{1}{2}(v_{\phi}x^{d}_{42})^{2}}$
(10) $\displaystyle M_{E_{1}}$ $\displaystyle\approx$
$\displaystyle\sqrt{(M^{L}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{L}_{34})^{2}},\qquad
M_{E_{2}}=\sqrt{(M^{e}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{e}_{43})^{2}+\frac{1}{2}(v_{\phi}x^{e}_{42})^{2}}.$
(11)

In the neutrino sector, the corresponding mass matrix is $7\times 7$ and its
explicit form can be found in Eq. (77). The resulting masses of the heavy
neutrinos read

$M_{N_{1}}=M_{N_{2}}\approx M^{\nu}_{4},\qquad
M_{N_{3}}=M_{N_{4}}\approx\sqrt{(M^{L}_{4})^{2}+\frac{1}{2}(v_{\phi}x^{L}_{34})^{2}}.$
(12)

By comparing Eq. (12) with Eq. (11), we can pinpoint two generic features of
the model considered in this study: heavy neutrinos $N_{1,2}$ are the lightest
VL leptons in the spectrum, while the pair $N_{3,4}$ is mass-degenerate (at
the tree level) with the charged VL lepton $E_{1}$. We will later see that
this mass pattern has important consequences for the resulting phenomenology.

As a final remark, let us notice that one complete VL family allows us to give
masses to the second and third family of the SM fermions only. To generate the
masses for the first family as well, one extra VL family is required (for an
example of such a construction, see Ref. Hernández et al. (2021)). Since such
an extension would only increase the number of free parameters in the model
without affecting any phenomenological findings, in this study we limit
ourselves to its most economical version.

###  II.2 CKM mixing matrix

The full $5\times 5$ mixing matrix takes the following form Cárcamo Hernández
et al. (2022):

$V_{\mathop{\mathrm{m}ixing}}=V_{L}^{u}.\mathop{\mathrm{d}iag}\left(1,1,1,1,0\right).V_{L}^{d\dagger}$
(13)

where $V_{L}^{u}$ and $V_{L}^{d}$ are the left-handed mixing matrices of Eqs.
(71) and (74) which diagonalize the up- and down-type quark mass matrices
$\mathcal{M}_{u}$ and $\mathcal{M}_{d}$. The zero element of the matrix (13)
indicates the fact that the singlet VL quarks do not interact with the SM
gauge bosons $W^{\pm}$. Following the strategy of Ref. Cárcamo Hernández et
al. (2022) and working under the assumption that $v_{u,d}/M_{4}^{Q,u,d}\ll 1$,
we can approximate the $3\times 3$ CKM matrix as

$V_{\mathop{\mathrm{C}KM}}^{3\times
3}\approx\left(\begin{array}[]{ccc}1-x_{ud}^{2}/2&x_{ud}&x_{ud}x_{d}\\\
-x_{ud}&1-x_{ud}^{2}/2&x_{d}-x_{u}\\\ -x_{u}x_{ud}&x_{u}-x_{d}&1\\\
\end{array}\right)\,,$ (14)

where

$x_{d}=\frac{y^{d}_{24}x^{d}_{43}M_{4}^{Q}}{y^{d}_{43}x^{Q}_{34}M_{4}^{d}}\,,\qquad\qquad\qquad
x_{u}=\frac{y^{u}_{24}x^{u}_{43}M_{4}^{Q}}{y^{u}_{43}x^{Q}_{34}M_{4}^{u}}\,,\qquad\qquad\qquad
x_{ud}=\frac{y^{d}_{14}}{y^{d}_{24}}\,.$ (15)

Based on the conclusions from Sec. II.1 one expects $x_{u},x_{d}\ll 1$. Note
also that:

  * •

The element $V_{us}$ of the CKM matrix is given by $y_{14}^{d}/y^{d}_{24}$ in
our model. The presence of a non-zero coupling $y_{14}^{d}$ is thus crucial to
generate the Cabibbo angle of the right size. We also expect
$y^{d}_{14}\approx 0.22\,y^{d}_{24}.$

  * •

The correct value of the element $V_{ud}$ is generated automatically once the
Cabibbo angle is set.

  * •

To reproduce the correct value of the element $V_{ub}$, one needs
$x_{d}\approx 0.017$. It then follows that $x_{u}\approx-0.023$ is required in
order to fit the element $V_{cb}$ (it also implies $x^{u}_{43}$ of order one).

  * •

The only element of the CKM matrix that can not be accurately reproduced is
$V_{td}$.

To analzye this issue more quantitatively, it is convenient to rewrite the CKM
matrix (14) in terms of the Wolfenstein parameters Wolfenstein (1983).
Defining, for example,

$x_{d}=A\,\lambda^{3}\,\sqrt{\eta^{2}+\rho^{2}}\,,\qquad\qquad
x_{u}=A\,\lambda^{2}(\sqrt{\eta^{2}+\rho^{2}}\,-1)\,,\qquad\qquad
x_{ud}=\frac{x_{u}\,x_{d}}{\lambda}\,,$ (16)

one obtains

$|V_{\mathop{\mathrm{C}KM}}^{3\times
3}|=\left(\begin{array}[]{ccc}1-\lambda^{2}/2&\lambda&A\,\lambda^{3}\,\sqrt{\eta^{2}+\rho^{2}}\\\
\lambda&1-\lambda^{2}/2&A\,\lambda^{2}\\\
A\,\lambda^{3}\left(1-\sqrt{\eta^{2}+\rho^{2}}\right)&A\,\lambda^{2}&1\\\
\end{array}\right)+\mathcal{O}(\lambda^{4})\,.$ (17)

Plugging the Wolfenstein parameters extracted from the global fit Workman et
al. (2022) into Eq. (17) and comparing it with the experimental determination
of the CKM matrix elements reported in Ref. Workman et al. (2022), one can
estimate to what extent the measured structure of the CKM matrix can be
reproduced in our model. One obtains

$\frac{|V_{\mathop{\mathrm{C}KM}}^{\mathop{\mathrm{e}xp}}|-|V_{\mathop{\mathrm{C}KM}}^{3\times
3}|}{\delta|V_{\mathop{\mathrm{C}KM}}^{\mathop{\mathrm{e}xp}}|}=\left(\begin{array}[]{ccc}0&0&0\\\
0&0.04&0\\\ 8.88&0.23&0.01\\\ \end{array}\right)\,.$ (18)

It results from Eq. (18) that in the framework of our model we may not be able
to correctly reproduce all the elements of the CKM matrix (this observation
will be later confirmed by our numerical scan). 111Note that modifying the
definitions of the parameters $x_{u}$, $x_{d}$ and $x_{ud}$ in Eq. (16), one
could be able to fit better the element $V_{td}$, but at the price of losing
the accuracy in reproducing $V_{cb}$. Once more, this issue could be solved by
introducing an extra VL family.

To conclude this section, we would like to stress again that the approximation
adopted in the foregoing discussion hinges on the assumption of the specific
mass hierarchy in the NP sector, which may not be entirely fulfilled.
Therefore, in the phenomenological analysis we will be always calculating all
the elements of the CKM matrix numerically.

##  III Scalar potential constraints

In this section, we discuss the constraints stemming from the scalar potential
of the model. In particular, we define the alignment limit of the SM-like
Higgs boson, we derive the conditions for the scalar potential to be bounded
from below in the presence of three independent scalar fields, and we verify
whether the electroweak (EW) vacuum is stable.

###  III.1 Scalar masses in the alignment limit

In the interaction basis, the most generic renormalizable scalar potential of
the model defined in Table 1 takes the form Hernández et al. (2021):

$\begin{split}V&=\mu_{u}^{2}(H_{u}^{\dagger}H_{u})+\mu_{d}^{2}(H_{d}^{\dagger}H_{d})+\mu_{\phi}^{2}(\phi^{*}\phi)-\frac{1}{2}\mu_{\mathop{\mathrm{s}b}}^{2}\left(\phi^{2}+\phi^{*2}\right)\\\
&+\frac{1}{2}\lambda_{1}(H_{u}^{\dagger}H_{u})^{2}+\frac{1}{2}\lambda_{2}(H_{d}^{\dagger}H_{d})^{2}+\lambda_{3}(H_{u}^{\dagger}H_{u})(H_{d}^{\dagger}H_{d})+\lambda_{4}(H_{u}^{\dagger}H_{d})(H_{d}^{\dagger}H_{u})\\\
&-\frac{1}{2}\lambda_{5}(\epsilon_{ij}H_{u}^{i}H_{d}^{j}\phi^{2}+\mathop{\mathrm{H}.c.})+\frac{1}{2}\lambda_{6}(\phi^{*}\phi)^{2}+\lambda_{7}(\phi^{*}\phi)(H_{u}^{\dagger}H_{u})+\lambda_{8}(\phi^{*}\phi)(H_{d}^{\dagger}H_{d}),\end{split}$
(19)

where $\mu_{u,d,\phi}^{2}$ are dimensionful mass parameters,
$\lambda_{1,2,\cdots,8}$ denote dimensionless quartic coupling constants, and
$\mu^{2}_{\mathop{\mathrm{s}b}}$ is an extra mass term which softly violates
the global U(1)X symmetry. The main reason to introduce the latter is to
prevent a massless Goldstone boson of the spontaneously broken U(1)X to appear
in the spectrum. As we will see below, the soft-breaking term does not affect
the CP-even and the charged scalar masses since it only enters the mass matrix
of the pseudoscalars.

Expanding the fields $H_{u}$, $H_{d}$ and $\phi$ around their vacuum states,

$\displaystyle H_{u}=\begin{pmatrix}H_{u}^{+}\\\
\frac{1}{\sqrt{2}}\left(v_{u}+\mathop{\mathrm{R}e}H_{u}^{0}+i\mathop{\mathrm{I}m}H_{u}^{0}\right)\end{pmatrix},\qquad
H_{d}=\begin{pmatrix}\frac{1}{\sqrt{2}}\left(v_{d}+\mathop{\mathrm{R}e}H_{d}^{0}+i\mathop{\mathrm{I}m}H_{d}^{0}\right)\\\
H_{d}^{-}\end{pmatrix},$
$\displaystyle\phi=\frac{1}{\sqrt{2}}\left(v_{\phi}+\mathop{\mathrm{R}e}\phi+i\mathop{\mathrm{I}m}\phi\right),\qquad\qquad\qquad\qquad\qquad\qquad$
(20)

where the vevs are defined in Eq. (3), one can use the minimization conditions
for the scalar potential (19) to express the dimensionful mass parameters in
terms of the quartic couplings and the vevs,

$\begin{split}\mu_{u}^{2}&=-\frac{1}{2}\left(\lambda_{1}v_{u}^{2}+\lambda_{3}v_{d}^{2}+\lambda_{7}v_{\phi}^{2}\right)-\frac{1}{4}\lambda_{5}\left(\frac{v_{d}}{v_{u}}\right)v_{\phi}^{2},\\\
\mu_{d}^{2}&=-\frac{1}{2}\left(\lambda_{2}v_{d}^{2}+\lambda_{3}v_{u}^{2}+\lambda_{8}v_{\phi}^{2}\right)-\frac{1}{4}\lambda_{5}\left(\frac{v_{u}}{v_{d}}\right)v_{\phi}^{2},\\\
\mu_{\phi}^{2}&=-\frac{1}{2}\left(\lambda_{6}v_{\phi}^{2}+\lambda_{5}v_{d}v_{u}+\lambda_{7}v^{2}_{u}+\lambda_{8}v_{d}^{2}\right)+\mu_{\mathop{\mathrm{s}b}}^{2}\,.\end{split}$
(21)

One must have

$\mu^{2}_{u}<0\,,\qquad\mu_{d}^{2}<0\,,\qquad\mu_{\phi}^{2}<0\,.$ (22)

in order to generate the non-zero vevs for all the scalar fields.

The explicit forms of the scalar mass matrices derived from the potential (19)
are collected in Appendix B. The real parts of the scalar fields,
$\mathop{\mathrm{R}e}H_{u}^{0}$, $\mathop{\mathrm{R}e}H_{d}^{0}$ and
$\mathop{\mathrm{R}e}\phi$, account for three CP-even Higgs bosons. The
corresponding mass matrix $\mathbf{M}_{\mathop{\mathrm{C}P-even}}^{2}$ (see
Eq. (79)) can be diagonalized with a mixing matrix $R_{h}$ defined in Eq.
(80). The masses of three physical neutral scalars, $h_{1}$, $h_{2}$ and
$h_{3}$, correspond to the eigenvalues of
$\mathbf{M}_{\mathop{\mathrm{C}P-even}}^{2}$,

${\mathrm{d}iag}\\{M_{h_{1}}^{2},M_{h_{2}}^{2},M_{h_{3}}^{2}\\}=R_{h}(\mathbf{M}_{\mathop{\mathrm{C}P-even}}^{2})R_{h}^{T}\,.$
(23)

In the following, we will want to identify the SM Higgs boson with the
lightest neutral scalar $h_{1}$. To this end, we choose to work in the so-
called alignment limit, defined as a set of constraints on the quartic
couplings $\lambda_{i}$ under which $h_{1}$ features the same tree-level
couplings with the SM particles as the SM Higgs. We show in Appendix B that
this assumption requires

$\displaystyle\lambda_{8}\,\cos^{2}\beta+\lambda_{7}\,\sin^{2}\beta+\lambda_{5}\,\sin\beta\cos\beta$
$\displaystyle=$ $\displaystyle 0$ (24)
$\displaystyle\lambda_{2}\,\cos^{2}\beta-\lambda_{1}\,\sin^{2}\beta-\lambda_{3}(\cos^{2}\beta-\sin^{2}\beta)$
$\displaystyle=$ $\displaystyle 0\,,$ (25)

where the equality imposes a perfect alignment condition. The masses of the
CP-even scalars in the alignment limit read

$\displaystyle M^{2}_{h_{1}}$ $\displaystyle=$ $\displaystyle
v^{2}\left(\lambda_{1}\,\sin^{2}\beta+\lambda_{3}\,\cos^{2}\beta\right)$ (26)
$\displaystyle M^{2}_{h_{2}}$ $\displaystyle=$
$\displaystyle\lambda_{6}\,v^{2}_{\phi}-\frac{1}{8\sin\beta\cos\beta}\left(B_{23}+\sqrt{4A_{23}^{2}+B_{23}^{2}}\right)$
(27) $\displaystyle M^{2}_{h_{3}}$ $\displaystyle=$
$\displaystyle\lambda_{6}\,v^{2}_{\phi}-\frac{1}{8\sin\beta\cos\beta}\left(B_{23}-\sqrt{4A_{23}^{2}+B_{23}^{2}}\right)\,,$
(28)

with $A_{23}$ and $B_{23}$ defined in Eq. (90) and Eq. (91), respectively, and
$v=\sqrt{v_{u}^{2}+v_{d}^{2}}=246$ GeV.

The CP-odd scalar mass matrix in the basis
$\left(\mathop{\mathrm{I}m}H_{u}^{0},\,\mathop{\mathrm{I}m}H_{d}^{0},\,\mathop{\mathrm{I}m}\phi\right)$,
$\mathbf{M}_{\mathop{\mathrm{C}P-odd}}^{2}$, is defined in Eq. (92). After the
diagonalization, the physical CP-odd spectrum consists of one massless
Goldstone boson and two massive pseudoscalars, $a_{1}$ and $a_{2}$,

${\mathrm{d}iag}\\{0,M_{a_{1}}^{2},M_{a_{2}}^{2}\\}=R_{a}(\mathbf{M}_{\mathop{\mathrm{C}P-odd}}^{2})R_{a}^{T}\,,$
(29)

with the masses given by

$\displaystyle M^{2}_{a_{1}}$ $\displaystyle=$
$\displaystyle-\frac{\lambda_{5}}{2\sin
2\beta}\left(v^{2}\sin^{2}2\beta+v^{2}_{\phi}\right)$ (30) $\displaystyle
M^{2}_{a_{2}}$ $\displaystyle=$ $\displaystyle
2\,\mu_{\mathop{\mathrm{s}b}}^{2}\,.$ (31)

Note that $\lambda_{5}<0$ and $\mu_{\mathop{\mathrm{s}b}}^{2}>0$ are required
to guarantee the positivity of $M^{2}_{a_{1}}$ and $M^{2}_{a_{2}}$.

Finally, the charged scalar mass matrix in the basis
$\left(H_{u}^{\pm},H_{d}^{\pm}\right)$,
$\mathbf{M}_{\mathop{\mathrm{C}harged}}^{2}$, is defined in Eq. (93). After
the diagonalization with a mixing matrix $R_{\beta}$, one is left with a
massless charged Goldstone boson and a charged Higgs boson,

${\mathrm{d}iag}\\{0,M_{h^{\pm}}^{2}\\}=R_{\beta}(\mathbf{M}_{\mathop{\mathrm{C}harged}}^{2})R_{\beta}^{T}\,.$
(32)

The corresponding mass reads in this case

$M^{2}_{h^{\pm}}=\frac{\lambda_{4}v^{2}}{2}-\frac{\lambda_{5}v_{\phi}^{2}}{2\sin
2\beta}\,.$ (33)

As a closing remark, let us notice that the alignment condition (25) indicates

$\lambda_{2}=\lambda_{3}+\tan^{2}\beta(\lambda_{1}-\lambda_{3})\,.$ (34)

In order to preserve the perturbativity of $\lambda_{2}$ (more on this in Sec.
V), the term in parentheses needs to be fine-tuned with a precision
$\mathcal{O}(1/\tan^{2}\beta)$ or better, effectively fixing
$\lambda_{3}\approx\lambda_{1}$ with the same accuracy. On the other hand, Eq.
(26) implies that we can identify $\lambda_{1}$ with the quartic coupling of
the SM, $\lambda_{1}=0.258$, as long as $\tan\beta\gtrsim 3$. Similarly, the
alignment condition (24) gives

$\lambda_{8}=-\tan\beta\left(\lambda_{7}\tan\beta+\lambda_{5}\right)\,.$ (35)

Perturbativity of $\lambda_{8}$ then requires
$\lambda_{7}\sim\mathcal{O}(1/\tan^{2}\beta)$ and
$\lambda_{5}\sim\mathcal{O}(1/\tan\beta)$.

###  III.2 Bounded-from-below limits

To guarantee that the minimum around which we expand the scalar potential (19)
is physically meaningful, we must ensure that the potential is bounded from
below, which means that it cannot tend to negative infinity along any
direction in the field space. This requirement puts additional restrictions on
the allowed values of the couplings $\lambda_{i}$. To derive the ‘bounded-
from-below’ constraints, one should analyze all possible directions along
which the scalar fields $H_{u}$, $H_{d}$ and $\phi$ can flow towards
arbitrarily large values. The details of our derivation are presented in
Appendix C. Here we summarize our findings in the form of inequality
conditions which need to be satisfied by the quartic couplings of the
potential (19):

$\begin{split}\lambda_{8}+\sqrt{\lambda_{2}\lambda_{6}}>0\\\
\lambda_{7}+\sqrt{\lambda_{1}\lambda_{6}}>0\\\
\lambda_{3}+\sqrt{\lambda_{2}\lambda_{1}}>0\\\
\lambda_{3}+\lambda_{4}+\sqrt{\lambda_{2}\lambda_{1}}>0\\\
-\frac{1}{4}\frac{(\mathop{\mathrm{R}e}\lambda_{5})^{2}+(\mathop{\mathrm{I}m}\lambda_{5})^{2}}{\lambda_{a}}+\lambda_{4}>0\\\
4\lambda_{b}^{2}-(\mathop{\mathrm{R}e}\lambda_{5})^{2}+\mathop{\mathrm{R}e}\lambda_{5}\mathop{\mathrm{I}m}\lambda_{5}>0\\\
4\lambda_{b}^{2}-(\mathop{\mathrm{I}m}\lambda_{5})^{2}+\mathop{\mathrm{R}e}\lambda_{5}\mathop{\mathrm{I}m}\lambda_{5}>0\end{split}$
(36)

where
$\lambda_{a}=\frac{3}{2}\lambda_{6}+\lambda_{3}\frac{\lambda_{6}}{\sqrt{\lambda_{1}\lambda_{2}}}+\lambda_{7}\sqrt{\frac{\lambda_{6}}{\lambda_{1}}}+\lambda_{8}\sqrt{\frac{\lambda_{6}}{\lambda_{2}}}$
and $\lambda_{b}=\sqrt{\lambda_{a}\lambda_{4}}$. Since in this study we do not
investigate the CP violation, we assume that all the parameters of the
lagrangian are real, indicating $\mathop{\mathrm{I}m}\lambda_{5}=0$. Note also
that several novel conditions w.r.t. the findings of Refs. Hernández et al.
(2021); Lee and Cárcamo Hernández (2022) are identified in Eq. (36).

###  III.3 Vacuum stability

In theories which feature an extended scalar sector, the scalar potential can
easily develop more than one local minimum. As a result, the theory may tunnel
from one minimum to another. In principle, color and charge breaking minima
deeper than the EWSB minimum of Eq. (3) can arise in our model (see, e.g.
Muhlleitner et al. (2017)). Moreover, several charge and color conserving
minima can coexist, in which case we do not know a priori which of them
corresponds to the desired EWSB minimum.

The strong vacuum stability condition for the scalar potential requires that
the EWSB vacuum corresponds to a global minimum. In such a case the potential
is said to be stable. If, on the other hand, the EWSB minimum is a local
minimum but the tunneling time to a true global minimum exceeds the age of the
Universe, the potential is said to be metastable. In this study we employ the
publicly available numerical package Vevacious++ Camargo-Molina and O’Leary
(2014) (the C++ version of Camargo-Molina et al. (2013)) to find all tree- and
one-loop level minima of the scalar potential defined in Eq. (19) and to
calculate the tunneling time from the EWSB minimum to the deepest minimum
found.

##  IV Flavor physics constraints

In this section, we review additional constraints which may affect the allowed
parameter space of the analyzed model. These extra restrictions come from the
experimental measurements of several flavor observables, including the
anomalous magnetic moment of the muon, the lepton flavor violating decays of
the tau lepton, and the elements of the CKM matrix. We discuss them in the
following one by one.

###  IV.1 Muon anomalous magnetic moment

The discrepancy between the SM prediction Czarnecki et al. (2003); Melnikov
and Vainshtein (2004); Aoyama et al. (2012); Kurz et al. (2014); Davier et al.
(2011); Gnendiger et al. (2013); Colangelo et al. (2014); Davier et al.
(2017); Masjuan and Sanchez-Puertas (2017); Colangelo et al. (2017);
Hoferichter et al. (2018); Keshavarzi et al. (2018); Colangelo et al. (2019);
Hoferichter et al. (2019); Davier et al. (2020); Keshavarzi et al. (2020);
Gérardin et al. (2019); Bijnens et al. (2019); Colangelo et al. (2020); Blum
et al. (2020); Aoyama et al. (2019, 2020) and the experimental measurement of
the anomalous magnetic moment of the muon has been confirmed separately by the
Brookhaven National Laboratory Bennett et al. (2006) and the Fermilab
experimental groups Abi et al. (2021); Aguillard et al. (2023), giving rise to
the combined $5.1\,\sigma$ anomaly:

$\Delta
a_{\mu}=a_{\mu}^{\mathop{\mathrm{e}xp}}-a_{\mu}^{\mathop{\mathrm{S}M}}=\left(2.49\pm
0.48\right)\times 10^{-9}.$ (37)

In a generic NP model which features heavy scalars $\phi_{i}$ and fermions
$\psi_{j}$ coupled to the SM muons via the Yukawa-type interactions
$y_{L}^{ij}\phi_{i}\,\bar{\psi}_{j}P_{L}\,\mu$ and
$y_{R}^{ij}\phi_{i}\,\bar{\psi}_{j}P_{R}\,\mu$ (where
$P_{L,R}=(1\mp\gamma^{5})/2$ are the usual projection operators), a well-known
one-loop contribution to the muon anomalous magnetic moment reads

$\Delta
a_{\mu}=\sum_{i,j}\left\\{-\frac{m_{\mu}^{2}}{16\pi^{2}M_{\phi_{i}}^{2}}\left(|y^{ij}_{L}|^{2}+|y^{ij}_{R}|^{2}\right)\left[Q_{j}\mathcal{F}_{1}\left(x_{ij}\right)-Q_{i}\mathcal{G}_{1}\left(x_{ij}\right)\right]\right.\\\
\left.-\,\frac{m_{\mu}\,M_{\psi_{j}}}{16\pi^{2}M_{\phi_{i}}^{2}}\textrm{Re}\left(y^{ij}_{L}y_{R}^{ij\ast}\right)\left[Q_{j}\mathcal{F}_{2}\left(x_{ij}\right)-Q_{i}\mathcal{G}_{2}\left(x_{ij}\right)\right]\right\\},$
(38)

where $M_{\phi_{i}}$ is the physical mass of a heavy scalar, $M_{\psi_{j}}$ is
the physical mass of a heavy fermion,
$x_{ij}=M_{\psi_{j}}^{2}/M_{\phi_{i}}^{2}$, and the electric charges of
$\phi_{i}$ and $\psi_{j}$ are related as $Q_{i}+Q_{j}=-1$. The loop functions
are defined in the following way:

$\begin{split}\mathcal{F}_{1}(x)&=\frac{1}{6\left(1-x\right)^{4}}\left(2+3x-6x^{2}+x^{3}+6x\ln
x\right)\\\
\mathcal{F}_{2}(x)&=\frac{1}{\left(1-x\right)^{3}}\left(-3+4x-x^{2}-2\ln
x\right)\\\
\mathcal{G}_{1}(x)&=\frac{1}{6\left(1-x\right)^{4}}\left(1-6x+3x^{2}+2x^{3}-6x^{2}\ln
x\right)\\\
\mathcal{G}_{2}(x)&=\frac{1}{\left(1-x\right)^{3}}\left(1-x^{2}+2x\ln
x\right)\,.\end{split}$ (39)

The first addend in Eq. (38) captures the loop chirality-conserving
contributions to $\Delta a_{\mu}$. These are known to be generically too small
to account for the anomaly (37) when the most recent LHC bounds on the NP
masses are taken into account Kowalska and Sessolo (2017); Athron et al.
(2021). We will thus focus on the second addend in Eq. (38), which corresponds
to the loop chirality-flipping contributions to $\Delta a_{\mu}$.

In the framework of the model defined in Table 1, two classes of contributions
to the anomalous magnetic moment of the muon can arise, induced by one-loop
diagrams with an exchange of neutral (pseudo)scalars and charged VL leptons,
as shown in Fig. 1(a), or charged scalars and neutral VL leptons, as shown in
Fig. 1(b). In the first case, the chirality-flipping contributions to $\Delta
a_{\mu}$ read

$\Delta
a_{\mu}^{Eh^{0}}\approx\frac{1}{16\pi^{2}}\sum_{j=1}^{2}\sum_{i=1}^{3}\left[\frac{m_{\mu}\,M_{E_{j}}}{M_{h^{0}_{i}}^{2}}\textrm{Re}\left(c_{L}c_{R}^{\ast}\right)^{E_{j},h_{i}^{0}}\mathcal{F}_{2}\left(M_{E_{j}}^{2}/M^{2}_{h^{0}_{i}}\right)\right]$
(40)

for the CP-even scalars and

$\Delta
a_{\mu}^{Ea}\approx\frac{1}{16\pi^{2}}\sum_{j=1}^{2}\sum_{i=1}^{2}\left[\frac{m_{\mu}\,M_{E_{j}}}{M_{a_{i}}^{2}}\textrm{Re}\left(c_{L}c_{R}^{\ast}\right)^{E_{j},a_{i}}\mathcal{F}_{2}\left(M_{E_{j}}^{2}/M^{2}_{a_{i}}\right)\right]$
(41)

for the CP-odd scalars. The one-loop contributions to $\Delta a_{\mu}$ from
the neutral leptons and charged scalars are given by

$\Delta
a_{\mu}^{Nh^{\pm}}\approx-\frac{1}{16\pi^{2}}\sum_{j=1}^{4}\left[\frac{m_{\mu}\,M_{N_{j}}}{M_{h^{\pm}}^{2}}\textrm{Re}\left(c_{L}c_{R}^{\ast}\right)^{N_{j},h^{\pm}}\mathcal{G}_{2}\left(M_{N_{j}}^{2}/M^{2}_{h^{\pm}}\right)\right].$
(42)

The parameters $c_{L/R}$ denote the effective couplings arising from the
muon-(pseudo)scalar-VL fermion vertices in the mass basis. They depend on the
lepton Yukawa couplings of Eqs. (1) and (76), as well as on the elements of
the mixing matrices $R_{h}$ (Eq. (23)), $R_{a}$ (Eq. (29)), and $V^{e}_{L/R}$
(Eq. (68)). The explicit forms of $c_{L/R}$ are rather complex and we refrain
from showing them here. Note, however, that in our numerical analysis we are
going to compute all the contributions to $\Delta a_{\mu}$ with the numerical
package SPheno Porod (2003); Porod and Staub (2012).

(a)

(b)

Figure 1: The one-loop chirality-flipping contributions to $\Delta a_{\mu}$
mediated by (a) a neutral (pseudo)scalar/charged lepton exchange, and (b) a
charged scalar/neutral lepton exchange.

###  IV.2 Lepton flavor violating decays

Due to the non-zero mixing between the second and the third generation of
fermions, charged lepton flavour violating processes may occur. The
$\tau\rightarrow\mu\gamma$ decay receives contributions from the one-loop
diagrams analogous to those of $\Delta a_{\mu}$. The corresponding branching
ratio (BR) is given by Mandal and Pich (2019)

${\mathrm{B}R}(\tau\to\mu\gamma)=\frac{\alpha_{\mathrm{em}}m_{\tau}^{3}}{4\,\Gamma_{\tau}}\sum_{i,j}\left(|A_{L}^{ij}|^{2}+|A_{R}^{ij}|^{2}\right)\,,$
(43)

where $\Gamma_{\tau}=2.3\times 10^{-12}$ Workman et al. (2022) indicates the
total decay width of the tau, $\alpha_{\mathrm{em}}$ is the fine structure
constant, and the decay amplitude $A_{L}^{ij}$ reads

$A_{L}^{ij}=\frac{1}{32\pi^{2}M_{\phi_{i}}^{2}}\left\\{m_{\tau}\left(y^{ij}_{\tau,L}y^{ij\ast}_{\mu,L}\right)\left[Q_{j}\mathcal{F}_{1}\left(x_{ij}\right)-Q_{i}\mathcal{G}_{1}\left(x_{ij}\right)\right]\right.\\\
\left.+M_{\psi_{j}}\left(y^{ij}_{\tau,L}y^{ij\ast}_{\mu,R}\right)\left[Q_{j}\mathcal{F}_{2}\left(x_{ij}\right)-Q_{i}\mathcal{G}_{2}\left(x_{ij}\right)\right]\right\\}\,.$
(44)

The corresponding amplitude $A_{R}^{ij}$ is obtained from Eq. (44) by
replacing $L\leftrightarrow R$. Just like it was in the $\Delta a_{\mu}$ case,
the main contribution to ${\mathrm{B}R}(\tau\to\mu\gamma)$ originates from the
second addend in Eq. (44). The current experimental 90% confidence level
(C.L.) upper bound on
$\mathop{\mathrm{B}R}\left(\tau\rightarrow\mu\gamma\right)$ from the Belle
collaboration reads Abdesselam et al. (2021):

$\mathop{\mathrm{B}R}\left(\tau\rightarrow\mu\gamma\right)_{\mathrm{exp}}<4.2\times
10^{-8}\,.$ (45)

The $\tau\to 3\,\mu$ decay can proceed through the one-loop penguin and box
diagrams. The latter are subdominant in our model as they do not receive the
chiral enhancement. The corresponding formulae for the penguin-diagram BRs are
lengthy and not particularly enlightening. They can be found, for example, in
Eq. (37) of Ref. Mandal and Pich (2019). The 90% C.L. upper bound on
$\mathop{\mathrm{B}R}\left(\tau\rightarrow 3\,\mu\right)$ by the Belle
collaboration reads Hayasaka et al. (2010):

$\mathop{\mathrm{B}R}\left(\tau\rightarrow
3\,\mu\right)_{\mathrm{exp}}<2.1\times 10^{-8}\,.$ (46)

###  IV.3 CKM anomaly

Among the experimental puzzles which are not explained by the SM we should
also mention various tensions between three different determinations of the
Cabibbo angle. This observable can be extracted from the short distance
radiative corrections to the $\beta$ decay, from the experimental data on kaon
decays, and from the lattice calculations Grossman et al. (2020); Belfatto and
Trifinopoulos (2023); Crivellin and Hoferichter (2020); Czarnecki et al.
(2020). All these measurements are in tension with each other, giving rise to
two interesting anomalies.

The first anomaly is related to the violation of the CKM matrix unitarity when
one compares the values of $\bigl{|}V_{ud}\bigr{|}$ and
$\bigl{|}V_{us}\bigr{|}$ resulting from the $\beta$ decay and from the kaon
decays. The second anomaly originates from two different measurements of
$\bigl{|}V_{us}\bigr{|}$: from the semileptonic $K\rightarrow\pi l\nu$ and the
leptonic $K\rightarrow\mu\nu$ decay, respectively.

The experimental upper bound on the CKM deviation from the unitarity reads
Workman et al. (2022)

$\Delta_{\mathop{\mathrm{C}KM}}=\sqrt{1-V_{ud}^{2}-V_{us}^{2}-V_{ub}^{2}}<0.05.$
(47)

To explain the anomaly of Eq. (47), one can consider extensions of the SM in
which the fermion sector is enlarged by VL quarks mixing at the tree level
with the SM quarks Belfatto and Trifinopoulos (2023); Cheung et al. (2020);
Branco et al. (2021); Albergaria and Branco (2023). In such a setting
deviations from the unitarity of the three-dimensional CKM matrix can arise
quite naturally. Since the model defined in Table 1 contains all the necessary
ingredients to account for the CKM anomaly, we include it in our list of
constraints.

##  V Perturbativity constraints

The model defined in Table 1 is intended as a phenomenological scenario which
correctly describes the physics around the energy scale determined by the
typical masses in the NP sector. Nevertheless, it is important to understand
what is the range of validity of such a model or, in other words, what is the
energy scale at which the model can not be trusted anymore and should be
embedded in some more fundamental UV completion. While such a “cut-off” scale
lacks a truly rigorous definition, one can try to estimate it by simply
requiring that whatever extra degrees of freedom emerge in the theory above
this scale to make the model UV complete, they do not affect its
phenomenological predictions.

As an example, let us consider the muon anomalous magnetic moment operator,
which in the low-energy effective filed theory (EFT) reads

$\frac{e}{2\,m_{\mu}}\Delta
a_{\mu}\left(\bar{\mu}\,\sigma_{\mu\nu}\,F^{\mu\nu}\mu\right)\equiv\frac{C}{\Lambda}\left(\bar{\mu}\,\sigma_{\mu\nu}\,F^{\mu\nu}\mu\right)\,.$
(48)

Here $\Lambda$ is a cut-off scale of the examined EFT while $C$ denotes a
generic Wilson coefficient. Note that since the operator in Eq. (48) is
chirality flipping, it is more convenient to define
$C=\tilde{C}\,m_{\mu}/\Lambda$. One can now derive from Eq. (48) rough
estimates of the energy scale associated with a hypothetical NP contributing
to $\Delta a_{\mu}$ at different loop orders,

$\displaystyle\textrm{tree level}:$ $\displaystyle\tilde{C}\approx
1,\qquad\quad\,\,\,\,$ $\displaystyle\Lambda\approx 3000\,\,\textrm{GeV}$ (49)
$\displaystyle\textrm{1 loop}:$ $\displaystyle\tilde{C}\approx
1/16\pi^{2},\quad\,$ $\displaystyle\Lambda\approx 230\,\,\textrm{GeV}$ (50)
$\displaystyle\textrm{2 loop}:$
$\displaystyle\tilde{C}\approx(1/16\pi^{2})^{2},$ $\displaystyle\Lambda\approx
20\,\,\textrm{GeV}$ (51)

and so on.

Going beyond the EFT approximation, let us investigate a one-loop chirality
flipping contribution to $\Delta a_{\mu}$ like the one in Eq. (38). Assuming
that it arises from an unspecified UV completion of our model above the scale
$\Lambda$, it can be estimated by the corresponding UV mass $m_{\Lambda}$ and
the UV Yukawa couplings $y_{L/R}(\Lambda)$ as

$\Delta
a_{\mu}^{\Lambda}\sim\frac{1}{16\pi^{2}}\frac{m_{\mu}\,v}{m_{\Lambda}^{2}}y_{L}(\Lambda)\,y_{R}(\Lambda).$
(52)

By demanding that the new contribution (52) does not shift our
phenomenological predictions for $\Delta a_{\mu}$ by more than $3\,\sigma$, we
can derive a lower bound on the UV mass,

$m_{\Lambda}\gtrsim\sqrt{y_{L}(\Lambda)y_{R}(\Lambda)}\,15\,\textrm{TeV}\,.$
(53)

For the Yukawa couplings at the upper edge of perturbativity,
$y_{L}(\Lambda)=y_{R}(\Lambda)=\sqrt{4\pi}$, Eq. (53) translates into a
conservative estimation of the scale of validity of our phenomenological
model,

$m_{\Lambda}\gtrsim 50\,\,\textrm{TeV}.$ (54)

In other words, the model can not be UV completed below $m_{\Lambda}$.

An immediate consequence of Eq. (54) is that the commonly employed
perturbativity bounds, which read $\lesssim\sqrt{4\pi}$ for the gauge/Yukawa
and $\lesssim 4\pi$ for the quartic couplings, need to be imposed on the
running parameters of the model evaluated at the scale $\Lambda$ rather than
on the bare couplings of the lagrangian (as it was done, for example, in Ref.
Hernández et al. (2021)).

To implement the RG-based perturbativity constraints, we follow the RG flow of
all the coupling constants from the scale $\mu_{0}=1.5$ TeV, which is a proxy
for the NP scale in our model, to $\Lambda=50$ TeV. The one-loop RG equations
(RGEs) were computed using the publicly available numerical code SARAH Staub
(2014, 2015) and are summarized in Appendix D. Due to a large number of Yukawa
and quartic interactions in our model, it is not possible to perform the
perturbativity analysis in a generic way as the RGEs are non-linear
differential equations that can not be solved analytically. On the other hand,
the perturbativity bounds are expected to be relevant only for those couplings
whose values must be of order $1$ (or larger) for phenomenological reasons.
This observation allows us to reduce the RGE system and to simplify the
analysis.

In the Yukawa sector, the couplings of interest are $x^{Q}_{34}$,
$y^{u}_{43}$, $y^{u}_{34}$ and $x^{u}_{43}$ (see Sec. II for the discussion).
We find that the modulus of their value cannot exceed $1.4$ at $\mu_{0}$ if
they are to remain perturbative up to 50 TeV. This conclusion is derived under
the assumption that all the other couplings (but two) are set to $1$ at the
initial scale $\mu_{0}$. The two exceptions are $y^{d}_{14}$ and $y^{e}_{24}$
(expected to be much smaller than 1 as the Yukawas of the second generation),
whose values at $\mu_{0}$ are set to $0.7$.

In the scalar sector, the perturbativity bounds are presumably most relevant
for the couplings $\lambda_{1}$, $\lambda_{6}$ and $\lambda_{7}$, whose RGEs
feature a power-four dependence on the large Yukawa couplings $y^{u}_{43}$,
$x^{u}_{43}$ and $x^{Q}_{34}$ (cf. Eq. (120), Eq. (125) and Eq. (126),
respectively). In Fig. 2 we illustrate the RG running of $\lambda_{1}$,
$\lambda_{6}$ and $\lambda_{7}$ for a randomly chosen benchmark point which
satisfies all the constraints discussed in Secs. II and III. The running of
all the remaining quartic couplings is very slow in the considered energy
range and does not pose any danger from the point of view of their
perturbativity. Once the whole system is analyzed with the alignment
conditions (24) and (25) in place, it turns out that the perturbativity
requires the modulus of the quartic couplings to be smaller than 2. A
straightforward consequence of this result is that all the benchmark points
found previously in Ref. Hernández et al. (2021) are disfavored.

Figure 2: The RG running of the quartic couplings $\lambda_{1}$, $\lambda_{6}$
and $\lambda_{7}$ for a randomly chosen benchmark point which satisfies all
the constraints discussed in Secs. II and III. The renormalization scale $\mu$
ranges from $1.5$ TeV to $1000\mathop{\mathrm{T}eV}$. $\mu_{0}=1.5$ TeV is a
reference scale. We do not show the RG evolution of other quartic and Yukawa
couplings as it is very slow in the considered energy range.

Finally, let us comment on another constraint which may arise in our model,
the so called perturbative unitarity. Although the $S$-matrix for a scattering
process must be unitary in the full theory, it may happen that at some order
in the perturbative expansion the unitarity is violated, signaling the
breakdown of the expansion. This is usually related to some of the couplings
becoming too large. The perturbative unitarity translates into conditions for
the partial wave amplitudes, which have to be smaller than 1/2. To examine
such constraints in our model we use SPheno, which computes the maximal
eigenvalue of a $2\rightarrow 2$ scattering matrix at the tree-level. On the
other hand, since we already require all the quartic couplings to remain
perturbative up to the energy scale of 50 TeV, we may suspect that the
perturbative unitarity bounds are automatically satisfied. As we will see in
the next section, this is indeed the case.

##  VI Numerical analysis and benchmark scenarios

In this section we perform a global numerical analysis of the model. We begin
by discussing the employed scanning methodology, the definition of the chi-
square ($\chi^{2}$) statistics and the initial ranges for all the model’s
parameters. Next, we present three best-fit benchmark scenarios which arise
from the minimization of the $\chi^{2}$ function. Finally, we provide a
discussion of some experimental signatures that these benchmark scenarios
could produce.

###  VI.1 Scanning methodology

Scalar sector $\tan\beta$ $\left[2,50\right]$ $v_{\phi}$
$\left[1000,1500\right]$ $\mu_{sb}^{2}$ $\left[4,64\right]\times 10^{4}$
$\lambda_{2}$ $\left[-2.0,+2.0\right]$ $\lambda_{3}$ $\left[0.24,0.28\right]$
$\lambda_{4}$ $\left[-2.0,+2.0\right]$ $\lambda_{5}$ $\left[-0.2,0.0\right]$
$\lambda_{6}$ $\left[-2.0,+2.0\right]$ $\lambda_{7}$
$\left[-0.01,+0.01\right]$ $\lambda_{8}$ $\left[-1.0,+1.0\right]$ Lepton
sector $y_{24}^{e}$ $\left[-0.7,+0.7\right]$ $y_{43}^{e}$
$\left[-1.0,+1.0\right]$ $y_{14}^{\nu}$ $\left[-1.0,+1.0\right]\times
10^{-10}$ $y_{14}^{\prime\nu}$ $\left[-1.0,+1.0\right]$ $M_{4}^{e}$
$\pm\left[200,1000\right]$ $y_{34}^{e}$ $\left[-1.0,+1.0\right]$ $x_{42}^{e}$
$\left[-1.0,+1.0\right]$ $y_{24}^{\nu}$ $\left[-1.0,+1.0\right]\times
10^{-10}$ $y_{24}^{\prime\nu}$ $\left[-1.0,+1.0\right]$ $M_{4}^{\nu}$
$\pm\left[200,1000\right]$ $x_{34}^{L}$ $\left[-1.0,+1.0\right]$ $x_{43}^{e}$
$\left[-1.0,+1.0\right]$ $y_{34}^{\nu}$ $\left[-1.0,+1.0\right]\times
10^{-10}$ $y_{34}^{\prime\nu}$ $\left[-1.0,+1.0\right]$ $M_{4}^{L}$
$\pm\left[200,1000\right]$ Quark sector $y_{24}^{u}$ $\left[-1.0,+1.0\right]$
$y_{43}^{u}$ $\left[-1.4,+1.4\right]$ $y_{14}^{d}$ $\left[-0.7,+0.7\right]$
$y_{43}^{d}$ $\left[-1.0,+1.0\right]$ $M_{4}^{d}$ $\pm\left[1200,4000\right]$
$y_{34}^{u}$ $\left[-1.4,+1.4\right]$ $x_{42}^{u}$ $\left[-1.0,+1.0\right]$
$y_{24}^{d}$ $\left[-1.0,+1.0\right]$ $x_{42}^{d}$ $\left[-1.0,+1.0\right]$
$M_{4}^{u}$ $\pm\left[1200,4000\right]$ $x_{34}^{Q}$ $\left[-1.0,+1.0\right]$
$x_{43}^{u}$ $\left[-1.4,+1.4\right]$ $y_{34}^{d}$ $\left[-1.0,+1.0\right]$
$x_{43}^{d}$ $\left[-1.0,+1.0\right]$ $M_{4}^{Q}$ $\pm\left[1200,4000\right]$

Table 2: Scanning ranges for the input parameters of the model defined in
Table 1. The alignment limit (cf. Sec. III), the RGE perturbativity
constraints (cf. Sec. V) and a tentative lower bound on the VL mass parameters
(see the text) are imposed. In the Yukawa sector only the non-zero couplings
are shown. Dimensionful quantities are given in GeV and GeV2.

In Table 2 we summarize the scanning ranges for all the parameters of the
model. These include the quartic couplings and the soft-breaking term of the
scalar potential (19), the non-zero Yukawa couplings and the mass parameters
of the lagrangian (1), the vev of the singlet scalar, and $\tan\beta$.

In the scalar sector the alignment conditions (24) and (25) are imposed,
leading to the limited scanning ranges for $\lambda_{3}$, $\lambda_{5}$ and
$\lambda_{7}$ (cf. Eqs. (34) and (35)). For all the other quartic couplings
the perturbativity bounds discussed in Sec. V are enforced. Similarly, the
Yukawa couplings are scanned in the ranges consistent with their RGE
perturbativity constraints. Finally, small values of some of the neutrino
coupling constants are necessary to generate tiny masses for the SM neutrinos.

A tentative lower bound of $1200\mathop{\mathrm{G}eV}$ is imposed on the VL
quark mass parameters. This is a rough (and conservative) approximation of the
constraints from the direct NP searches at the LHC, which will be discussed in
more details in Sec. VII.1. The scanning range for $v_{\phi}$ then follows
from the requirement of reproducing the correct mass of the top quark, as
discussed in Sec. II.1. Similarly, we adopt $200$ GeV lower bounds on the VL
lepton mass parameters in order to be roughly consistent with the
corresponding LHC constraints, which we examine in Sec. VII.2. Finally, the
range for $\mu^{2}_{sb}$ was chosen to make sure that the mass of the
associated CP-odd state (cf. Eq. (31)) is not excluded by the current
experimental searches Workman et al. (2022).

The experimental constraints employed in our numerical scan are listed in
Table 3. The central values and the experimental errors for the quark and
lepton masses and for the CKM matrix elements are quoted after the PDG report
Workman et al. (2022). Since the uncertainties for $m_{\mu}$ and $m_{\tau}$
are very small, rendering the fitting procedure numerically challenging, we
adopt an error of 10% for these two observables. The experimental constraints
from the flavor physics were discussed in Sec. IV.

We construct the $\chi^{2}$-statistic function as

$\chi^{2}=\sum_{i}\frac{\left(\mathcal{O}_{i}^{\mathop{\mathrm{m}odel}}-\mathcal{O}_{i}^{\mathop{\mathrm{c}en}}\right)^{2}}{(\mathcal{O}_{i}^{\mathop{\mathrm{e}rr}})^{2}}\,,$
(55)

where $\mathcal{O}_{i}^{\mathop{\mathrm{m}odel}}$ indicates the value of an
observable calculated in our model, $\mathcal{O}_{i}^{\mathop{\mathrm{c}en}}$
is the central value of its experimental measurement,
$\mathcal{O}_{i}^{\mathop{\mathrm{e}rr}}$ is the corresponding experimental
error, and the sum runs over all the measured observables listed in Table 3.
The upper bounds, corresponding to the last three rows of Table 3, are not
included in the $\chi^{2}$ function, but applied as hard-cuts instead (a point
in the parameter space is rejected if such a condition is not satisfied).

Measurement |  Central Value |  Exp. Error |  Measurement |  Central Value |  Exp. Error
---|---|---|---|---|---
$m_{\mu}$ |  $0.10566$ |  $10\%$ |  $|V_{ud}|$ |  $0.97370$ |  $0.00014$
$m_{\tau}$ |  $1.77686$ |  $10\%$ |  $|V_{us}|$ |  $0.22450$ |  $0.00080$
$m_{c}$ |  $1.270$ |  $0.020$ |  $|V_{ub}|$ |  $0.00382$ |  $0.00024$
$m_{s}$ |  $0.0934$ |  $0.0034$ |  $|V_{cd}|$ |  $0.22100$ |  $0.00400$
$m_{b}$ |  $4.18$ |  $0.02$ |  $|V_{cs}|$ |  $0.98700$ |  $0.01100$
$m_{t}$ |  $172.76$ |  $0.30$ |  $|V_{cb}|$ |  $0.04100$ |  $0.00140$
$\Delta a_{\mu}$ |  $2.49\times 10^{-9}$ |  $0.48\times 10^{-9}$ |  $|V_{td}|$ |  $0.00800$ |  $0.00030$
|  |  |  $|V_{ts}|$ |  $0.03880$ |  $0.00110$
|  |  |  $|V_{tb}|$ |  $1.01300$ |  $0.03000$
Measurement |  Upper bound |  |  |  |
$\mathop{\mathrm{B}R}\left(\tau\rightarrow\mu\gamma\right)$ |  $<4.2\times 10^{-8}$ |  |  |  |
$\mathop{\mathrm{B}R}\left(\tau\rightarrow 3\mu\right)$ |  $<2.1\times 10^{-8}$ |  |  |  |
$\Delta_{\mathop{\mathrm{C}KM}}$ |  $<$0.05 |  |  |  |
Table 3: The experimental measurements which we employ in our numerical scan.
Masses are in GeV.

To minimize the $\chi^{2}$ function, we adopt the following strategy. First,
we perform an initial scan of the parameter space consistent with Table 2. As
a result, we obtain a seed which is then used to minimize the $\chi^{2}$
function by iterating a random walk algorithm with an adaptive step function.
The step function is chosen such that at each iteration all input parameters
are updated by less than $\kappa\%$, and $\kappa$ reduces with an exponential
decay law throughout the minimization procedure. During each iteration, we
discard all the points that do not satisfy the upper bounds on
$\Delta_{\text{CKM}}$, $\textrm{BR}(\tau\rightarrow\mu\gamma)$ and
$\textrm{BR}(\tau\rightarrow 3\mu)$, as well as the boundedness constraints on
the scalar potential given in Eq. (36) and the perturbative unitarity.
Moreover, we investigate the vacuum stability with Vevacious++ Camargo-Molina
and O’Leary (2014) and we keep only those points whose vacuum is identified as
“stable”.

###  VI.2 Benchmark scenarios

Scalar sector BP1 BP2 BP3 BP1 BP2 BP3 $\tan\beta$ $13$ $8$ $12$ $\lambda_{1}$
$0.258$ $0.258$ $0.258$ $v_{u}$ $245.3$ $244.3$ $245.2$ $\lambda_{2}$ $0.514$
$0.153$ $0.623$ $v_{d}$ $18.9$ $30.5$ $20.4$ $\lambda_{3}$ $0.257$ $0.260$
$0.256$ $v_{\phi}$ $1015$ $1077$ $1012$ $\lambda_{4}$ $0.552$ $0.304$ $0.167$
$\mu_{u}^{2}$ $-7.8\times 10^{3}$ $-6.6\times 10^{3}$ $-7.6\times 10^{3}$
$\lambda_{5}$ $-0.039$ $-0.072$ $-0.061$ $\mu_{d}^{2}$ $-8.2\times 10^{3}$
$-8.6\times 10^{4}$ $-3.4\times 10^{4}$ $\lambda_{6}$ $0.370$ $0.487$ $0.663$
$\mu_{\phi}^{2}$ $-4.9\times 10^{4}$ $-9.4\times 10^{4}$ $-2.3\times 10^{5}$
$\lambda_{7}$ $0.001$ $0.002$ $0.002$ $\mu_{\mathop{\mathrm{s}b}}^{2}$
$1.4\times 10^{5}$ $1.9\times 10^{5}$ $1.1\times 10^{5}$ $\lambda_{8}$ $0.254$
$0.423$ $0.417$ Quark sector Lepton sector BP1 BP2 BP3 BP1 BP2 BP3
$y_{24}^{u}$ $-0.051$ $-0.049$ $0.050$ $y_{24}^{e}$ $0.028$ $-0.015$ $0.022$
$y_{34}^{u}$ $-0.980$ $1.185$ $-1.024$ $y_{34}^{e}$ $-0.895$ $0.612$ $0.790$
$x_{34}^{Q}$ $0.924$ $-0.842$ $-0.877$ $x_{34}^{L}$ $0.616$ $-0.729$ $0.724$
$y_{43}^{u}$ $1.382$ $1.093$ $-1.337$ $y_{43}^{e}$ $-0.223$ $0.144$ $-0.191$
$x_{42}^{u}$ $0.550$ $0.821$ $-0.595$ $x_{42}^{e}$ $0.156$ $0.165$ $0.188$
$x_{43}^{u}$ $1.286$ $1.261$ $1.263$ $x_{43}^{e}$ $-0.168$ $0.228$ $-0.205$
$y_{14}^{d}$ $-0.022$ $0.035$ $0.026$ $y_{14}^{\nu}$ $-2\times 10^{-11}$
$5\times 10^{-11}$ $3\times 10^{-11}$ $y_{24}^{d}$ $0.096$ $0.151$ $-0.113$
$y_{24}^{\nu}$ $3\times 10^{-11}$ $8\times 10^{-12}$ $6\times 10^{-11}$
$y_{34}^{d}$ $-0.684$ $0.274$ $0.267$ $y_{34}^{\nu}$ $-5\times 10^{-11}$
$9\times 10^{-11}$ $9\times 10^{-11}$ $y_{43}^{d}$ $-0.672$ $-0.489$ $0.656$
$y_{14}^{\prime\nu}$ $-0.824$ $-0.674$ $-0.674$ $x_{42}^{d}$ $-0.371$ $-0.110$
$0.225$ $y_{24}^{\prime\nu}$ $-0.895$ $-0.874$ $-0.896$ $x_{43}^{d}$ $-0.160$
$0.072$ $-0.127$ $y_{34}^{\prime\nu}$ $0.701$ $0.744$ $-0.812$ Mass parameters
BP1 BP2 BP3 BP1 BP2 BP3 $M_{4}^{u}$ $-1317$ $1405$ $1334$ $M_{4}^{e}$ $-517$
$-575$ $533$ $M_{4}^{d}$ $-3644$ $3068$ $-2882$ $M_{4}^{\nu}$ $204$ $-212$
$217$ $M_{4}^{Q}$ $-1384$ $1443$ $1322$ $M_{4}^{L}$ $-206$ $-222$ $-202$

Table 4: Input parameters for three best-fit benchmark scenarios. Dimensionful
quantities are given in GeV and GeV2.

In Table 4 we present input parameters for three best-fit benchmark scenarios
identified by performing the numerical scan discussed in Sec. VI.1. The
corresponding mass spectra are summarized in Table 5 while the breakdown of
individual contributions to the $\chi^{2}$ function is shown in Table 6.

SM fermions BP1 BP2 BP3 BP1 BP2 BP3 $m_{c}$ $1.262$ $1.282$ $1.259$ $m_{\mu}$
$0.110$ $0.110$ $0.110$ $m_{t}$ $172.7$ $172.8$ $172.6$ $m_{\tau}$ $1.864$
$1.756$ $1.765$ $m_{s}$ $0.089$ $0.093$ $0.091$
$m_{\nu_{2}}\left[10^{-10}\right]$ $4.659$ $6.587$ $0.252$ $m_{b}$ $4.169$
$4.196$ $4.175$ $m_{\nu_{3}}\left[10^{-10}\right]$ $8.253$ $18.38$ $20.95$ NP
fermions Quark sector Lepton sector BP1 BP2 BP3 BP1 BP2 BP3 $M_{U_{1}}$ $1495$
$1561$ $1440$ $M_{E_{1}}$ $487$ $596$ $554$ $M_{U_{2}}$ $1708$ $1842$ $1704$
$M_{E_{2}}$ $543$ $615$ $570$ $M_{D_{1}}$ $1534$ $1579$ $1464$ $M_{N_{1,2}}$
$205$ $214$ $218$ $M_{D_{2}}$ $3655$ $3070$ $2888$ $M_{N_{3,4}}$ $488$ $598$
$556$ Scalars BP1 BP2 BP3 BP1 BP2 BP3 $M_{h_{1}}$ $125$ $125$ $125$
$M_{a_{1}}$ $362$ $411$ $433$ $M_{h_{2}}$ $362$ $412$ $435$ $M_{a_{2}}$ $532$
$614$ $469$ $M_{h_{3}}$ $617$ $752$ $824$ $M_{h^{\pm}}$ $384$ $423$ $440$

Table 5: Mass spectra for three best-fit benchmark scenarios. All masses are
in GeV.

In general, the three benchmark scenarios demonstrate quite similar features,
both in terms of the input parameters and of the resulting NP spectra. This is
largely due to the fact that we aim at reproducing masses and mixings of the
SM fermions and this, as we discussed in Sec. II, puts strong constraints on
(some of) the model’s parameters.

Let us first notice that the masses of all the SM fermions of the third and
second generation can be fitted very precisely. Each individual contribution
to the $\chi^{2}$ function is smaller than $0.7$, with an exception of $m_{s}$
in BP1, in which case we have $\chi^{2}_{b}=1.6$. We can also observe that, as
we anticipated in Sec. II, the Yukawa couplings which link the VL sector with
the SM fermions of the third generation are in general larger than those
associated with the second generation.

On the other hand, fitting the CKM matrix is a little bit more tricky and the
bulk of the total $\chi^{2}$ stems from this very sector. As anticipated in
Sec. II.2, the main contribution to the $\chi^{2}$ function is given by the
element $|V_{td}|$, with the corresponding $\chi_{V_{td}}^{2}$ ranging from
$16$ for BP1 to $25$ for BP3. Smaller yet still relevant contributions come
from the entries $|V_{ub}|$ and $|V_{ts}|$. Finally, an order 10 contribution
to the $\chi^{2}$ function from the element $|V_{us}|$ is mainly due to a very
small experimental error associated with this particular observable. All other
elements of the CKM matrix are fitted within their $1\sigma$ experimental
ranges. As an illustration, we present below the full $5\times 5$ CKM matrix
for the benchmark scenario BP1,

$|V_{\mathop{\mathrm{C}KM}}|^{(\textrm{BP1})}=\left(\begin{array}[]{ccc|cc}0.97394&0.22679&0.00298&1.4\times
10^{-7}&0.00008\\\ 0.22671&0.97301&0.04296&0.00003&0.00042\\\
0.00681&0.04236&0.99821&0.00968&0.00221\\\ \hline\cr
0.00054&0.00270&0.02853&0.86532&0.00114\\\
0.00082&0.00390&0.02996&0.50113&0.00076\\\ \end{array}\right)\,.$ (56)

The two remaining best-fit points follow the same pattern. Incidentally, note
that the CKM anomaly is $\mathcal{O}(10^{-4})$ in our setup, well below the
experimental upper bound of Eq. (47).

Interestingly, in all three cases each quartic (Yukawa) coupling remains
smaller than $4\pi$ ($\sqrt{4\pi}$) up to $1000\mathop{\mathrm{T}eV}$. We can
therefore conclude that the validity range of our model extends well beyond
the putative scale of $50\mathop{\mathrm{T}eV}$. We also checked that the
maximal eigenvalue of the scattering matrix computed by SPheno is
$\mathcal{O}(10^{-2})$ for all the benchmark scenarios, indicating that the
perturbative unitarity bound is satisfied as well.

Quarks masses CKM elements BP1 BP2 BP3 BP1 BP2 BP3 $\chi_{c}^{2}$ $0.154$
$0.360$ $0.280$ $\chi_{V_{us}}^{2}$ $8.225$ $8.290$ $5.986$ $\chi_{t}^{2}$
$0.022$ $0.018$ $0.119$ $\chi_{V_{ub}}^{2}$ $12.33$ $10.36$ $9.327$
$\chi_{s}^{2}$ $1.569$ $0.014$ $0.450$ $\chi_{V_{td}}^{2}$ $15.69$ $18.30$
$24.58$ $\chi_{b}^{2}$ $0.330$ $0.640$ $0.052$ $\chi_{V_{ts}}^{2}$ $10.45$
$9.796$ $6.559$ $\chi_{Q}^{2}$ $2.075$ $1.031$ $0.901$ $\chi_{V}^{2}$ $55.45$
$55.26$ $55.92$ Charged leptons masses $\Delta a_{\mu}$ $\chi_{\mu}^{2}$
$0.210$ $0.170$ $0.207$ $\chi_{\Delta a_{\mu}}^{2}$ $0.328$ $0.657$ $0.375$
$\chi_{\tau}^{2}$ $0.241$ $0.014$ $0.004$ Total $\chi_{L}^{2}$ $0.451$ $0.183$
$0.211$ $\chi^{2}_{\textrm{TOT}}$ $58.30$ $57.31$ $57.41$

Table 6: Breakdown of the $\chi^{2}$ contributions from various observables
implemented in the $\chi^{2}$ function of Eq. (55). The CKM contributions
which are smaller than 3 are not shown. $\chi^{2}_{Q}$, $\chi^{2}_{L}$ and
$\chi^{2}_{V}$ indicate total $\chi^{2}$ contributions from the quark masses,
lepton masses, and the CKM matrix elements, respectively.
$\chi^{2}_{\textrm{TOT}}$ stands for the total $\chi^{2}$ function of each
best-fit scenario.

Contributions to $\Delta a_{\mu}\times 10^{9}$ Charged scalars CP-even scalars
Loop BP1 BP2 BP3 Loop BP1 BP2 BP3 $h^{\pm},N_{1,2}$ $-1.076\,\,\,\,$
$-0.792\,\,\,\,$ $-0.942\,\,\,\,$ $h_{1},E_{1}$ $-0.003\,\,\,\,$
$-0.001\,\,\,\,$ $-0.009\,\,\,\,$ $h^{\pm},N_{3,4}$ $3.300$ $2.898$ $3.153$
$h_{1},E_{2}$ $0.003$ $0.001$ $0.009$ $h^{\pm},N_{\text{tot}}$ $2.225$ $2.106$
$2.211$ $h_{2},E_{1}$ $-0.409\,\,\,\,$ $-0.520\,\,\,\,$ $-0.969\,\,\,\,$ CP-
odd scalars $h_{2},E_{2}$ $0.437$ $0.548$ $0.994$ $a_{1},E_{1}$ $0.425$
$0.528$ $0.938$ $h_{3},E_{1}$ $0.018$ $0.115$ $0.076$ $a_{1},E_{2}$
$-0.544\,\,\,\,$ $-0.611\,\,\,\,$ $-1.529\,\,\,\,$ $h_{3},E_{2}$
$-0.017\,\,\,\,$ $-0.127\,\,\,\,$ $-0.076\,\,\,\,$ $a_{2},E_{1}$
$-0.033\,\,\,\,$ $-0.135\,\,\,\,$ $-0.071\,\,\,\,$ $h,E_{\mathrm{tot}}$
$0.032$ $0.027$ $0.025$ $a_{2},E_{2}$ $0.110$ $0.196$ $0.621$ Total
$a,E_{\mathrm{tot}}$ $-0.015\,\,\,\,$ $-0.023\,\,\,\,$ $-0.041\,\,\,\,$
$\Delta a_{\mu}$ $2.215$ $2.101$ $2.196$

Table 7: Contributions to $\Delta a_{\mu}$ from the individual one-loop
diagrams shown in Fig. (1). The subscript ”tot” indicates the sum of all the
contributions of a given type.

Masses of the NP leptons are determined, to a large extent, by correctly
fitting the experimental value of $\Delta a_{\mu}$ (an overall contribution
from this observable to the total $\chi^{2}$ function does not exceed 0.7 in
all the benchmark scenarios). Contributions to $\Delta a_{\mu}$ from the
individual one-loop diagrams of Fig. (1) are summarized in Table 7. We present
separately fractions of $\Delta a_{\mu}$ generated by the charged scalars
$h^{\pm}$ and the neutral leptons $N_{1,2,3,4}$, by the CP-odd scalars
$a_{1,2}$ and the charged leptons $E_{1,2}$, and by the CP-even scalars
$h_{1,2,3}$ and the charged leptons $E_{1,2}$. We also show the sum of all the
contributions of a given type, indicated by a subscript ”tot”.

We observe that the largest contributions to $\Delta a_{\mu}$ arise from the
charged scalar/heavy neutrino loops. We thus disprove the conclusions of Refs.
Hernández et al. (2021); Lee and Cárcamo Hernández (2022) where it was assumed
that the charged lepton loops were the only NP contributions to muon $(g-2)$
present in the model. As we show in our analysis, all possible one-loop
diagrams contributing to $\Delta a_{\mu}$ should be treated at equal footing
and none of them should be discarded a priori.

Even more interestingly, the observed dominance of the heavy neutrino
contributions to the anomalous magnetic moment of the muon seems to be a
generic feature of the model which does not pertain exclusively to the
identified benchmark scenarios. The charged scalar/heavy neutrino loops are
determined, among other parameters, by combinations of the $y^{\prime\nu}$
couplings which are not constrained by any SM masses and mixing and thus can
become relatively large. Contrarily, the same Yukawa coupling is responsible
for the generation of the neutral (pseudo)scalar/charged lepton loops and for
the correct tree-level mass of the muon. It is thus required to be small and
the corresponding contributions to $\Delta a_{\mu}$ are suppressed.
Additionally, one also observes cancellations between the individual
contributions to $\Delta a_{\mu}$ stemming from the (pseudo)scalar diagrams
with different VL leptons, which is a known and common feature of many NP
models with the VL fermions (see, e.g., Darmé et al. (2018); Kowalska and
Sessolo (2021) for a discussion).

Finally, we should mention the size of the BRs for the lepton flavor violating
decays $\tau\rightarrow\mu\gamma$ and $\tau\to 3\mu$ in our model, which are
of the order $(3-4)\times 10^{-8}$ for the former and $(6-9)\times 10^{-10}$
for the latter. Given that in the future the Belle-II collaboration is
expected to improve their current exclusion bounds by an order of magnitude or
more Kou et al. (2019); Aggarwal et al. (2022), it may turn out that the tau
leptonic decays offer the best experimental way of verifying the predictions
of the NP model analyzed in this study.

##  VII LHC study of the benchmark scenarios

In this section we confront the benchmark scenarios identified in Sec. VI with
the null results of the direct NP searches at the LHC. We analyze, one by one,
the constraints originating from considering the production of VL quarks, VL
leptons, and exotic scalars.

###  VII.1 Vector-like quarks

The VL quarks (VLQs) can be copiously produced at the LHC, either in pairs
through the strong interactions or singly through an exchange of the EW gauge
bosons. In the former case, the dominant production channels at the leading
order are gluon fusion and quark-antiquark annihilation, whose production
cross sections depend on the VLQ mass and its SU(3)C quantum numbers only (see
Refs. Kowalska and Kumar (2019); Olivas et al. (2022) for analytical
formulae). Therefore, the experimental lower bounds on VLQ masses are expected
to be, to a large extent, model independent, baring only a slight dependence
on the relative strength of the individual VLQ decay channels.

Figure 3: Pair production of the VLQs $T$ and $B$ via the gluon fusion at the
LHC considered by the ATLAS collaboration in Ref. Aad et al. (2023a).

The most recent analysis from ATLAS, based on the data from proton–proton
collisions at a centre-of-mass energy of $\sqrt{s}=13$ TeV, corresponding to
an integrated luminosity of $139$ fb-1 Aad et al. (2023a), considered the pair
production of VL top partners $T$ and VL bottom partners $B$ with the decay
channels $T\to Zt,\,ht,\,Wb$ and $B\to Zb,\,hb,\,Wt$ and with large missing
transverse momentum. The corresponding Feyman diagrams are depicted in Fig. 3.
The strongest 95% C.L. lower bounds on the VLQ mass derived in Ref. Aad et al.
(2023a) read

$M_{T/B}>1.41\mathop{\mathrm{T}eV}$ (57)

for the EW doublets222For the VLQ mass larger than 800 GeV this indicates
BR$(T\to Zt)$=BR$(T\to ht)=50\%$ and BR$(B\to Wt)=100\%$ Aad et al. (2023a).
and

$M_{T}>1.26\mathop{\mathrm{T}eV},\qquad M_{B}>1.33\mathop{\mathrm{T}eV}$ (58)

for the EW singlets.333For the VLQ mass larger than 800 GeV this indicates
BR$(T\to Zt)=25\%$, BR$(T\to ht)=25\%$, BR$(T\to Wb)=50\%$, BR$(B\to
Zb)=25\%$, BR$(B\to hb)=25\%$ and BR$(B\to Wt)=50\%$ Aad et al. (2023a). The
analogous results from CMS can be found in Ref. Tumasyan et al. (2023).
Similarly, by assuming that at least one of the VLQs decays into a $Z$ boson
with the BR=100%, the $13$ TeV ATLAS search Aad et al. (2023b) obtained even
stronger bounds,

$M_{T}>1.60\mathop{\mathrm{T}eV},\qquad M_{B}>1.42\mathop{\mathrm{T}eV}\,.$
(59)

At first glance it may seem that all our benchmark scenarios are consistent
with the VLQ exclusion bounds. On the other hand, in the model considered in
this study the couplings of the physical heavy quarks $U_{1,2}$ and $D_{1,2}$
with the third-generation SM quarks and the EW gauge (Higgs) bosons are
generated via tree-level mixing after the EW and U(1)X symmetries are
spontaneously broken (cf. Sec. II and Appendix A). Since there is a priori no
reason for the resulting BRs to correspond to any of the benchmark cases
considered by ATLAS and CMS in their analyses (exotic decays to the charged
Higgs are possible, for example), we need to reexamine the experimental
results in the framework of our model.

To this end, we calculate with MadGraph5 MC@NLO Alwall et al. (2014) the cross
sections for the pair production of $U_{1}$, $U_{2}$, $D_{1}$ and $D_{2}$. The
results are presented in Table 8. By comparing these numbers with the observed
experimental 95% C.L. upper bounds on the signal cross section from Ref. Aad
et al. (2023a) (to give an example,
$\sigma^{\textrm{exp}}_{95\%}(p\,p\to\bar{T}T)=4\times 10^{-3}$ pb for
$M_{\textrm{VLQ}}=1.5$ TeV) we conclude that our benchmark scenarios are
indeed not excluded by the current LHC searches for the VLQs, irrespectively
of the actual sizes of their BRs.

$M_{U_{1}}$ $\sigma(pp\rightarrow U_{1}U_{1})$ $M_{U_{2}}$
$\sigma(pp\rightarrow U_{2}U_{2})$ $M_{D_{1}}$ $\sigma(pp\rightarrow
D_{1}D_{1})$ $M_{D_{2}}$ $\sigma(pp\rightarrow D_{2}D_{2})$ BP1 $1495$
$1.3\times 10^{-3}$ $1708$ $3.9\times 10^{-4}$ $1534$ $1.0\times 10^{-3}$
$3655$ $4.5\times 10^{-9}$ BP2 $1561$ $8.9\times 10^{-4}$ $1842$ $1.9\times
10^{-4}$ $1579$ $8.1\times 10^{-4}$ $3070$ $1.6\times 10^{-7}$ BP3 $1440$
$1.8\times 10^{-3}$ $1704$ $4.0\times 10^{-4}$ $1464$ $1.6\times 10^{-3}$
$2888$ $5.0\times 10^{-7}$

Table 8: Cross sections (in pb) for the pair production of the VLQs for our
three benchmark scenarios. Masses are in GeV.

It is instructive to investigate the prospects of testing our model in future
runs at the LHC. The total cross section for the VLQ pair production, followed
by a decay into the third generation quarks and the EW gauge/Higgs bosons, can
be expressed using the narrow width approximation (NWA) as

$\widetilde{\sigma}\left(pp\rightarrow\overline{Q}Q\rightarrow
f\overline{f}\,VV\right)\approx\sigma\left(pp\rightarrow
Q\overline{Q}\right)\mathop{\mathrm{B}R}\left(Q\rightarrow
fV\right)\mathop{\mathrm{B}R}\left(\overline{Q}\rightarrow\overline{f}V\right),$
(60)

where $Q=U_{1,2},D_{1,2}$, $f=t,b$ and $V=W,Z,h_{1}$. Under the assumption
that $\mathop{\mathrm{B}R}\left(Q\rightarrow
fV\right)=\mathop{\mathrm{B}R}\left(Q\rightarrow
h_{1}t/b\right)+\mathop{\mathrm{B}R}\left(Q\rightarrow
Zt/b\right)+\mathop{\mathrm{B}R}\left(Q\rightarrow Wb/t\right)=1$, the cross
section (60) reduces to the signal cross section constrained by the
experimental collaborations,
$\sigma^{\textrm{exp}}_{95\%}\approx\sigma\left(pp\rightarrow
Q\overline{Q}\right)$. If, on the other hand, the three BRs do not sum to one,
we expect the resulting exclusion bounds to be weaker than the bounds reported
by ATLAS and CMS.

The lightest VLQ in our model, $U_{1}$, is characterized by the following BRs:

$\displaystyle\mathop{\mathrm{B}P1}:\quad$
$\displaystyle\textrm{BR}\left(U_{1}\rightarrow
h_{1}t\right)=0.188,\quad\textrm{BR}\left(U_{1}\rightarrow
Zt\right)=0.146,\quad\textrm{BR}\left(U_{1}\rightarrow Wb\right)=0.040$
$\displaystyle\mathop{\mathrm{B}P2}:\quad$
$\displaystyle\textrm{BR}\left(U_{1}\rightarrow
h_{1}t\right)=0.135,\quad\textrm{BR}\left(U_{1}\rightarrow
Zt\right)=0.101,\quad\textrm{BR}\left(U_{1}\rightarrow Wb\right)=0.044$ (61)
$\displaystyle\mathop{\mathrm{B}P3}:\quad$
$\displaystyle\textrm{BR}\left(U_{1}\rightarrow
h_{1}t\right)=0.209,\quad\textrm{BR}\left(U_{1}\rightarrow
Zt\right)=0.165,\quad\textrm{BR}\left(U_{1}\rightarrow Wb\right)=0.037$

with the resulting cross sections $\widetilde{\sigma}^{\textrm{BP1}}=1.8\times
10^{-4}$ pb, $\widetilde{\sigma}^{\textrm{BP2}}=7.0\times 10^{-5}$ pb and
$\widetilde{\sigma}^{\textrm{BP3}}=3.0\times 10^{-4}$ pb. The second-to-the-
lightest VLQ, $D_{1}$, has

$\displaystyle\mathop{\mathrm{B}P1}:\quad$
$\displaystyle\textrm{BR}\left(D_{1}\rightarrow
h_{1}b\right)=0.001,\quad\textrm{BR}\left(D_{1}\rightarrow
Zb\right)=0.001,\quad\textrm{BR}\left(D_{1}\rightarrow Wt\right)=0.375$
$\displaystyle\mathop{\mathrm{B}P2}:\quad$
$\displaystyle\textrm{BR}\left(D_{1}\rightarrow
h_{1}b\right)=0.001,\quad\textrm{BR}\left(D_{1}\rightarrow
Zb\right)=0.001,\quad\textrm{BR}\left(D_{1}\rightarrow Wt\right)=0.221$ (62)
$\displaystyle\mathop{\mathrm{B}P3}:\quad$
$\displaystyle\textrm{BR}\left(D_{1}\rightarrow
h_{1}b\right)=0.001,\quad\textrm{BR}\left(D_{1}\rightarrow
Zb\right)=0.001,\quad\textrm{BR}\left(D_{1}\rightarrow Wt\right)=0.387$

with $\widetilde{\sigma}^{\textrm{BP1}}=1.4\times 10^{-4}$ pb,
$\widetilde{\sigma}^{\textrm{BP2}}=4.1\times 10^{-5}$ pb and
$\widetilde{\sigma}^{\textrm{BP3}}=2.4\times 10^{-4}$ pb. We can thus conclude
that in order to probe the VL masses featured by our benchmark scenarios, at
least one order of magnitude enhancement of the experimental sensitivity in
the VLQ searches is required.

Finally, we analyze the possibility of testing the model via processes in
which the VLQs are produced one at the time. The single VL $T$ quark
production was analyzed by ATLAS in Refs. Aad et al. (2023c, d, 2022), while
the single VL $B$ quark production in Ref. Aad et al. (2023e). The
corresponding Feynman diagrams are shown in Fig. 4. The 95% C.L. experimental
upper bounds on the relevant signal cross sections are of the order
$(10^{-2}-10^{-1})$ pb.

We calculated the cross sections for the VLQ single productions of our three
benchmark points using the NWA. The hadronic cross sections were obtained with
MadGraph5 MC@NLO and the BRs with SPheno. We found that the cross section for
a single production of the VLQs $U_{1}$ is $\mathcal{O}(10^{-5})$ pb, while
for $D_{1}$ it amounts to $\mathcal{O}(10^{-7})$ pb. We can thus conclude that
in our model the single production is a less promising search strategy than
the pair production. This was to be expected as the single production is in
general less competitive than the pair production for the Yukawa couplings
smaller than 1, see e.g. Faroughy et al. (2017); Kowalska et al. (2019).

Figure 4: Single production of the VLQs $T$ and $B$ via the EW gauge boson
exchange at the LHC, as considered by the ATLAS collaboration in Refs. Aad et
al. (2023c, d, 2022) for $T$ and in Ref. Aad et al. (2023e) for $B$.

###  VII.2 Vector-like leptons

At the tree level, the VL leptons (VLLs) are pair produced at the LHC via the
Drell-Yan processes. The corresponding cross sections for our three benchmark
scenarios are collected in Table 9. The analysis of all the possible
experimental signatures is in this case much more involved than for the VLQs
as the lepton decay BRs strongly depend on the presence in the spectrum of the
exotic scalars lighter than the VLLs. The following mass hierarchies are
observed:

$\displaystyle\mathop{\mathrm{B}P1}:$ $\displaystyle
M_{N_{1,2}}<M_{h_{2}},M_{a_{1}},M_{h^{\pm}}<M_{E_{1}},M_{N_{3,4}}<M_{a_{2}}<M_{E_{2}}<M_{h_{3}}$
$\displaystyle\mathop{\mathrm{B}P2}:$ $\displaystyle
M_{N_{1,2}}<M_{h_{2}},M_{a_{1}},M_{h^{\pm}}<M_{E_{1}},M_{N_{3,4}}<M_{a_{2}}<M_{E_{2}}<M_{h_{3}}$
(63) $\displaystyle\mathop{\mathrm{B}P3}:$
$\displaystyle\,\,M_{N_{1,2}}<M_{h_{2}},M_{a_{1}},M_{h^{\pm}}<M_{a_{2}}<M_{E_{1}},M_{N_{3,4}}<M_{E_{2}}<M_{h_{3}}\,.$

$M_{E_{2}}$ $\sigma(pp\rightarrow\bar{E}_{2}E_{2})$ $M_{E_{1}}$
$\sigma(pp\rightarrow\bar{E}_{1}E_{1})$ $M_{N_{3,4}}$ $\sigma(pp\rightarrow
N_{3,4}N_{3,4})$ $\sigma(pp\rightarrow E_{1}N_{3,4})$ BP1 $543$ $1.5\times
10^{-3}$ $487$ $5.7\times 10^{-3}$ $488$ $2.7\times 10^{-5}$ $3.6\times
10^{-3}$ BP2 $615$ $8.0\times 10^{-4}$ $596$ $1.8\times 10^{-3}$ $598$
$8.5\times 10^{-6}$ $1.2\times 10^{-3}$ BP3 $570$ $1.2\times 10^{-3}$ $554$
$2.8\times 10^{-3}$ $556$ $2.8\times 10^{-5}$ $1.8\times 10^{-3}$

Table 9: Cross sections (in pb) for the pair production of the VLLs for our
three benchmark scenarios. Masses are in GeV.

In all the cases the lightest VLLs, neutrinos $N_{1,2}$, originate
predominantly from the SU(2)L singlets and their production cross section at
the LHC is suppressed, $\mathcal{O}(10^{-5})$ pb.

The second-to-the-lightest VLLs, $E_{1}$ and heavy neutrinos $N_{3,4}$, come
from the same SU(2)L doublets and are almost degenerate in mass. Therefore,
three production channels should be considered simultaneously:
$p\,p\rightarrow Z/\gamma\rightarrow E_{1}\bar{E}_{1}$, $p\,p\rightarrow
Z/\gamma\rightarrow N_{3,4}N_{3,4}$ and $p\,p\rightarrow W^{\pm}\rightarrow
E_{1}N_{3,4}$. The dominant branching ratios for the subsequent decays of
$E_{1}$ and $N_{3,4}$, evaluated with SPheno, are collected in Table 10. In
all three cases the VLLs decay predominantly to the SM muons, which is a
direct consequence of the fact that we impose $\Delta a_{\mu}$ as a constraint
in our likelihood function and the largish muon-lepton-scalar Yukawa couplings
are preferred.

BR BP1 BR BP2 BR BP3 $E_{1}\to\mu\,a_{1}$ 37% $E_{1}\to\mu\,a_{1}$ 23%
$E_{1}\to\mu\,a_{1}$ 21% $E_{1}\to\mu\,h_{2}$ 37% $E_{1}\to\mu\,h_{2}$ 24%
$E_{1}\to\mu\,h_{2}$ 25% $E_{1}\to N_{1,2}\,W^{\pm}$ 26% $E_{1}\to\tau\,a_{1}$
11% BR BP1 BR BP2 BR BP3 $N_{3,4}\to\mu\,h^{\pm}$ 70% $N_{3,4}\to\mu\,h^{\pm}$
51% $N_{3,4}\to\mu\,h^{\pm}$ 56% $N_{3,4}\to N_{1,2}\,Z$ 12% $N_{3,4}\to
N_{1,2}\,h_{1}$ 11%

Table 10: Dominant BRs for the decays of $E_{1}$ and $N_{3,4}$.

A closer look at Table 10 reveals that the relative strengths of various VLL
decay channels are, to some extent, scenario dependent. Moreover, the final
experimental signatures hinge on the subsequent decay channels of the scalar
particles, which are also pretty complex (we discuss it in more details in
Sec. VII.3). As an example, let us consider a process $p\,p\rightarrow
E_{1}\,\bar{E}_{1}\to\mu\,\bar{\mu}\,a_{1}\,a_{1}$ for the benchmark scenario
BP1. The lightest pseudoscalar can decay in this case either to a $b\,\bar{b}$
pair (with the BR of 28%) or to $\nu\,N_{1,2}$ (with the BR of 69%). The decay
of the heavy neutrinos then proceeds as $N_{1,2}\to
e^{\pm}/\mu^{\pm}\,W^{\pm}$ with the BR of 56%, or $N_{1,2}\to\nu\,Z$ with the
BR of 28%. We can thus expect the following distinctive experimental
signatures emerging from the $p\,p\rightarrow
E_{1}\,\bar{E}_{1}\to\mu\,\bar{\mu}\,a_{1}\,a_{1}$ process: a) 2 muons +
(b)-jets, b) multileptons + missing energy, c) multileptons + jets + missing
energy, with the total signal cross section reduced w.r.t. the production
cross section reported in Table 9 by the product of the subsequent BRs.

To make the things even worse, there are not many LHC analysis that would
explicitly look for the VLLs. The only dedicated ATLAS search, based on the
139 $\mathop{\mathrm{f}b}^{-1}$ of data from the 13 TeV run Aad et al.
(2023f), looks for the VLLs coupled predominantly to taus. The analogous CMS
analysis based on the 77.4 $\mathop{\mathrm{f}b}^{-1}$ of data can be found in
Ref. Sirunyan et al. (2019a). The decay chains considered by the two
collaborations are shown in Fig. 5. In both cases, the total cross sections
for the VLL production can be probed down to $10^{-3}$ pb (of course, the
actual value is mass dependent).

Figure 5: Pair production of the VLLs $E_{1}$ and $N_{3,4}$ at the LHC, as
considered by the ATLAS Aad et al. (2023f) and CMS Sirunyan et al. (2019a)
collaborations.

In our model all three benchmark scenarios feature very low BRs for $E_{1}$
and $N_{3,4}$ decaying to taus, which do not exceed 10%. We can thus expect a
strong suppression of the resulting signal w.r.t. the experimental analysis.
Indeed, using the NWA the total cross section for the process considered in
Refs. Aad et al. (2023f) and Sirunyan et al. (2019a) can be written as
follows:

$\displaystyle\sigma\left(p\,p\rightarrow\tau^{-}\tau^{+}l^{-}l^{+}q\,\overline{q}\right)$
$\displaystyle\approx$ $\displaystyle\sigma\left(p\,p\rightarrow
E_{1}\overline{E}_{1}\right)\mathop{\mathrm{B}R}\left(E_{1}\rightarrow\tau^{-}l^{-}l^{+}\right)\mathop{\mathrm{B}R}\left(\overline{E}_{1}\rightarrow\tau^{+}q\,\overline{q}\right)$
(64) $\displaystyle+$ $\displaystyle\sigma\left(p\,p\rightarrow
E_{1}N_{3,4}\right)\mathop{\mathrm{B}R}\left(E_{1}\rightarrow\tau^{-}l^{-}l^{+}\right)\mathop{\mathrm{B}R}\left(N_{3,4}\rightarrow\tau^{+}q\,\overline{q}\right)\,.$

Combining the cross sections from Table 9 with the relevant BRs calculated
with SPheno, we obtain

$\displaystyle\mathop{\mathrm{B}P1}:$
$\displaystyle\sigma\left(p\,p\rightarrow\tau^{-}\tau^{+}l^{-}l^{+}q\,\overline{q}\right)$
$\displaystyle=8.3\times 10^{-7}\,\mathop{\mathrm{p}b}$
$\displaystyle\mathop{\mathrm{B}P2}:$
$\displaystyle\sigma\left(p\,p\rightarrow\tau^{-}\tau^{+}l^{-}l^{+}q\,\overline{q}\right)$
$\displaystyle=8.1\times 10^{-6}\,\mathop{\mathrm{p}b}$ (65)
$\displaystyle\mathop{\mathrm{B}P3}:$
$\displaystyle\sigma\left(p\,p\rightarrow\tau^{-}\tau^{+}l^{-}l^{+}q\,\overline{q}\right)$
$\displaystyle=5.3\times 10^{-7}\,\mathop{\mathrm{p}b}.$

If we now compare the predictions of Eq. (VII.2) with the corresponding
experimental 95% C.L. exclusion cross sections from Ref. Aad et al. (2023f),
we can conclude that the benchmark scenarios identified in Sec. VI are not
excluded by the dedicated LHC searches for the VLLs.444In principle, some of
the SUSY searches looking for the chargino/neutralino production, e.g. Ref.
Aad et al. (2023g), analyze signatures that could be generated in our model.
However, the resulting cross section for such a process is way too low to
allow the derivation of any constraints. Moreover, it may also be challenging
to test the VLL sector of our model in future runs at the LHC, if no dedicated
experimental strategies for the muon final state signatures are proposed.

###  VII.3 Exotic scalars

Process |  BP1 |  BP2 |  BP3
---|---|---|---
$a_{1}\to\nu\,N_{1,2}$ |  69% |  82% |  72%
$a_{1}\to\bar{b}\,b$ |  28% |  14% |  25%
$a_{2}\to\bar{t}\,t$ |  73% |  73% |  36%
$a_{2}\to\bar{b}\,b$ |  |  |  15%
$a_{2}\to\nu\,N_{1,2}$ |  |  |  46%
$h_{2}\to\nu\,N_{1,2}$ |  69% |  84% |  72%
$h_{2}\to\bar{b}\,b$ |  28% |  14% |  25%
$h_{3}\to\bar{t}\,t$ |  50% |  43% |  34%
$h_{3}\to\bar{\tau}\,E_{1}$ |  10% |  9% |  20%
$h^{\pm}\to\bar{\mu}\,N_{1,2}$ |  41% |  51% |  48%
$h^{\pm}\to\bar{e}\,N_{1,2}$ |  34% |  30% |  26%
$h^{\pm}\to\bar{b}\,t$ |  19% |  13% |  19%
Table 11: Dominant BRs ($>5\%$) for the decays of the exotic scalars.

Finally, we investigate the possibility of testing the predictions of our
model via the LHC searches in the scalar sector. There is a plethora of
experimental analyses, both by ATLAS and CMS, that look for the non-SM Higgs
bosons (see Ref. Kling et al. (2020) for a recent review in the framework of
the 2HDM). At the same time, in our benchmark scenarios the exotic scalars can
decay through a variety of channels (the dominant BRs, obtained with SPheno,
are collected in Table 11).

To facilitate the analysis, we use the publicly available code HiggsTools Bahl
et al. (2023), a toolbox for evaluating bounds from the direct searches for
the exotic scalar particles at LEP and the LHC, whose database contains 258
different limits. We find that all our benchmark scenarios are tagged as
“allowed” by HiggsTools.

It is instructive to take a closer look at the output of HiggsTools, as it
indicates which searches are most sensitive to the spectra featured by our
best-fit scenarios. This is quantified by a parameter called “observed ratio”,
$R_{\textrm{obs}}$, which is the ratio of the predicted cross section and the
experimental limit at the 95% C.L. The point in the parameter space is
excluded if $R_{\textrm{obs}}>1$. We observe that the highest values of
$R_{\textrm{obs}}$ (0.6 for BP1 and BP3, 0.14 for BP2) are reached for the
$h_{2}\to\tau^{+}\tau^{-}$ and $a_{1}\to\tau^{+}\tau^{-}$ decays constrained
by the ATLAS $139\mathop{\mathrm{f}b}^{-1}$ analysis Aad et al. (2020a),
despite very low decay BRs in this channel.

To investigate it in more details, we calculated the
$a_{1}/h_{2}\to\tau^{+}\tau^{-}$ cross sections with MadGraph5 MC@NLO. The
results are reported in the last three columns of Table 12. These are to be
compared with the 95% C.L. experimental lower bounds on the cross section
reported in the third column of Table 12. We find a very good agreement with
the output of HiggsTools in terms of the parameter $R_{\textrm{obs}}$, thus
confirming that the decays of the exotic scalars into taus are going to be the
most promising way of testing the predictions of the model at the LHC.

In Table 12 we also present other decay channels of $a_{1}$ and $h_{2}$ that
feature the high sensitivity. While the current experimental bounds on those
searches are weaker that those relative to the $\tau^{-}\tau^{+}$ final state,
they may offer complementary signatures of the model in future LHC runs.

Incidentally, note that the BRs for the decays of $h_{2}$ and $a_{1}$ to the
EW gauge bosons, $\gamma\gamma$, $ZZ$ and $WW$, are $\mathcal{O}(10^{-8})$ and
$\mathcal{O}(10^{-9})$, respectively, which is orders of magnitude below the
current experimental bounds.

Finally, let us comment on the possibility of testing the
$a_{1}/h_{2}\rightarrow t^{+}t^{-}$ decay through the measurement of an
effective coupling $g_{a_{1}/h_{2}tt}$. This scenario was investigated by CMS
in Ref. Sirunyan et al. (2020a). Comparing the values of the
$g_{a_{1}/h_{2}tt}$ coupling evaluated with SPheno (0.084 for BP2, 0.096 for
BP3) with the experimental 95% C.L. upper bounds (0.80 for BP2 and 0.70 for
BP3)555The CMS analysis Sirunyan et al. (2020a) does not cover the scalar
masses below 400 GeV. we conclude that no additional constraint on our model
arises from this particular search.

Channel |  Experiment |  $\sigma^{\mathop{\mathrm{e}xp}}_{95\%}$ (BP1, BP2, BP3) |  $\sigma^{\textrm{BP1}}$ |  $\sigma^{\textrm{BP2}}$ |  $\sigma^{\textrm{BP3}}$
---|---|---|---|---|---
$a_{1}/h_{2}\rightarrow\tau^{+}\tau^{-}$ |  CMS Sirunyan et al. (2018a) |  $0.060\quad 0.030\quad 0.020$ |  $0.037$ |  $0.004$ |  $0.011$
|  ATLAS Aad et al. (2020a) |  $0.050\quad 0.020\quad 0.016$ |  |  |
$a_{1}/h_{2}\rightarrow\mu^{+}\mu^{-}$ |  CMS Sirunyan et al. (2019b) |  $0.007\quad 0.006\quad 0.005$ |  $1.3\times 10^{-4}$ |  $1.4\times 10^{-5}$ |  $4.4\times 10^{-5}$
|  ATLAS Aaboud et al. (2019) |  $0.009\quad 0.004\quad 0.003$ |  |
$a_{1}/h_{2}\rightarrow b^{+}b^{-}$ |  CMS Sirunyan et al. (2018b) |  $6.0\qquad 3.5\qquad 3.0$ |  $0.554$ |  $0.061$ |  $0.061$
|  ATLAS Aad et al. (2020b) |  $-\qquad-\qquad-$ |  |  |
$h^{\pm}\rightarrow\bar{t}b$ |  CMS Sirunyan et al. (2020b) |  $0.40\qquad 0.30\qquad 0.25$ |  $7.6\times 10^{-3}$ |  $2.1\times 10^{-3}$ |  $4.5\times 10^{-3}$
|  ATLAS Aad et al. (2021) |  $0.45\qquad 0.30\qquad 0.25$ |  |  |
Table 12: An overview of the LHC scalar searches which present the highest
sensitivity to the benchmark scenarios identified in Sec. VI. The columns
show, respectively, the decay channel, the experimental analyses investigating
this channel, the experimental 95% C.L. upper bound on the cross section for
the mass corresponding to the mass of the scalar in a benchmark scenario, the
actual cross section calculated in the benchmark scenario.

##  VIII Conclusions

In this study, we performed a global analysis of an extension of the SM which
contains one full family of VL fermions, an extra SU(2)L scalar doublet and an
SU(2)L scalar singlet. It also features a U(1)X global symmetry spontaneously
broken by the singlet scalar vev. This scenario was originally proposed in
Ref. King (2018) to generate the masses of the third and the second family of
the SM fermions, as well as to account correctly for their mixing patterns.

In our analysis we confronted the model with the experimental bounds from the
flavor physics observables, which include the anomalous magnetic moment of the
muon and the rare decays of the tau lepton. Additionally, the model was
subjected to the theoretical constraints stemming from the stability of the
scalar potential and from the perturbativity of the renormalized couplings.
Importantly, we revisited and corrected the bounded-from-below and the
alignment limits, which in the context of the same model were previously
discussed in Refs. Hernández et al. (2021); Lee and Cárcamo Hernández (2022).
In particular, we showed that additional constraints on the quartic couplings
arise if the full three-scalar potential is considered. We also argued that
the perturbativity bounds should not be imposed on the low-scale parameters of
the lagrangian but on the running couplings evaluated at the renormalization
scale which sets an upper limit of the model’s validity. These RG-based
perturbativity conditions require the low-scale scalar couplings to be smaller
than 2 and the Yukawa couplings smaller than 1.5.

With all the constraints in place, we performed a numerical scan of the
model’s parameter space and we identified three benchmark scenarios that
satisfied all the theoretical and experimental requirements. One distinctive
feature of these solutions is that the charged scalar/heavy neutrino loops
provide dominant contributions to the observable $\Delta a_{\mu}$. This
finding is qualitatively different from the conclusions obtained in Refs.
Hernández et al. (2021); Lee and Cárcamo Hernández (2022) where only the
charged lepton loops were considered. We would like to emphasize that the
dominance of the heavy neutrino contribution to $\Delta a_{\mu}$ is a generic
characteristic of the model and not a mere artifact of the specific benchmark
scenarios. The main reason behind this feature is that the same coupling which
generates the neutral scalar/charged lepton loops is also responsible for the
correct tree-level mass of the muon and thus it is required to be small.

We also performed a detailed LHC analysis of our three best-fit scenarios. We
investigated the experimental constraints stemming from the direct searches
for VLQs, VLLs and exotic scalars. We found that none of the currently
available exclusion bounds can test the spectra featured by the benchmark
scenarios. This provides a proof of concept that the model in study is
feasible as an explanation of both the SM masses and mixings and of the
relevant experimental phenomena.

Regarding future prospects for experimental verification of the model, several
observations can be made. Firstly, both charged and neutral VLLs decay
predominantly to muons in our framework. On the other hand, all currently
available LHC analyses focus on taus in the final state, for which the cross
sections obtained in our model are several orders of magnitude below the
experimental sensitivity. Therefore, we would like to encourage the
experimental collaborations to provide dedicated analyses of the VLLs coupled
to the second family of the SM fermions. Such a study would not only allow to
test the predictions of our model, but it would prove very useful in any
phenomenological research that aims at explaining the muon $(g-2)$ anomaly in
a NP framework with the VLLs.

Secondly, the experimental searches for the VLQs can become a fruitful testing
ground for our model already in the current run of the LHC. The cross sections
for the pair production of VLQs featured by the benchmark scenarios are one
order of magnitude smaller than the current experimental upper bounds and
should be in reach of the dedicated VLQs searches based on the larger data
samples.

Finally, we observed that the most constraining decay channel for the exotic
scalars is $a_{1}/h_{2}\rightarrow\tau^{+}\tau^{-}$, for which the ratio of
the predicted to the experimental cross sections is close to 1. It may thus
provide complementary signatures of our model in future runs at the LHC.

However, the ultimate verification of the NP scenario considered in this study
may come from the flavor physics. The Belle-II collaboration plans on
improving, by at least one order of magnitude, their experimental bounds on
the rare leptonic decays of the tau lepton. As the corresponding branching
ratios featured by our three benchmark scenarios are very close to the current
90% C.L. exclusion limits, the rare decays could be the first experimental
signatures to be tested.

## Acknowledgements

K.K. would like to thank Enrico Sessolo for discussions and comments on the
manuscript. A.E.C.H is supported by ANID-Chile FONDECYT 1210378, ANID
PIA/APOYO AFB180002 and Milenio-ANID-ICN2019_044. K.K., H.L. and D.R. are
supported by the National Science Centre (Poland) under the research Grant No.
2017/26/E/ST2/00470. The use of the CIS computer cluster at the National
Centre for Nuclear Research in Warsaw is gratefully acknowledged.

##  Appendix A Fermion mass matrices

###  A.1 Charged leptons

The mass matrix for the charged leptons, $\mathcal{M}_{e}$, can be derived
from Eq. (1) after identifying the generic fermions $\psi$ with the
corresponding lepton fields from Table 1 and the generic scalar $H$ with
$H_{d}$. One thus has

$\psi_{iR}=e_{iR},\qquad\psi_{iL}=L_{iL},\qquad\psi_{4R}=e_{4R},\qquad\psi_{4L}=L_{4L},\qquad\widetilde{\psi}_{4R}=\widetilde{L}_{4R},\qquad\widetilde{\psi}_{4L}=\widetilde{e}_{4L}\,$
(66)

with the following components of the SU(2)L doublets:
$L_{iL}=(\nu_{iL},e_{iL})^{T}$, $L_{4L}=(\nu_{4L},e_{4L})^{T}$ and
$\widetilde{L}_{4R}=(\widetilde{\nu}_{4R},\widetilde{e}_{4R})^{T}$. As a
result, the mass matrix reads

$\mathcal{M}_{e}=\left(\begin{array}[]{c|ccccc}&e_{1R}&e_{2R}&e_{3R}&e_{4R}&\widetilde{e}_{4R}\\\\[2.15277pt]
\hline\cr e_{1L}&0&0&0&0&0\\\\[4.30554pt]
e_{2L}&0&0&0&y_{24}^{e}\frac{v_{d}}{\sqrt{2}}&0\\\\[4.30554pt]
e_{3L}&0&0&0&y_{34}^{e}\frac{v_{d}}{\sqrt{2}}&-x_{34}^{L}\frac{v_{\phi}}{\sqrt{2}}\\\\[4.30554pt]
e_{4L}&0&0&y_{43}^{e}\frac{v_{d}}{\sqrt{2}}&0&-M_{4}^{L}\\\\[4.30554pt]
\widetilde{e}_{4L}&0&x_{42}^{e}\frac{v_{\phi}}{\sqrt{2}}&x_{43}^{e}\frac{v_{\phi}}{\sqrt{2}}&M_{4}^{e}&0\\\
\end{array}\right)\,,$ (67)

where $M_{4}^{L}$ ($M_{4}^{e}$) denotes the mass of the VL lepton doublet
(singlet) and $x^{L}_{34}\equiv x^{e}_{34}$. To facilitate the comparison with
the corresponding mass matrix defined in SARAH, we adopt the sign convention
used in the code. Note, however, that such a choice does not affect the
conclusions drawn in our study as we allow all the Yukawa couplings and all
the VL mass parameters to assume both positive and negative values in our
numerical scan.

The $5\times 5$ charged lepton mass matrix $\mathcal{M}_{e}$ can be
diagonalized by means of two unitary transformations $V_{L}^{e}$ and
$V_{R}^{e}$,

$V_{L}^{e}\mathcal{M}_{e}V_{R}^{e\dagger}=\mathop{\mathrm{d}iag}\left(0,m_{\mu},m_{\tau},M_{E_{1}},M_{E_{2}}\right).$
(68)

In Sec. II the approximate expressions for the eigenvalues $m_{\mu}$ and
$m_{\tau}$ were provided in Eq. (6), whereas the analogous formulae for the
eigenvalues $M_{E_{1}}$ and $M_{E_{2}}$ were given in Eq. (11). While those
equations are very useful to get a general idea on which lagrangian parameters
are relevant for generating the physical charged lepton masses, in our
numerical analysis we diagonalize all the fermion mass matrices numerically,
employing the SPheno code generated by SARAH.

###  A.2 Up-type quarks

In analogy to the charged lepton sector, the mass matrix for the up-type
quarks, $\mathcal{M}_{u}$, can be derived from Eq. (1) after taking $H=H_{u}$
and making the following identification:

$\psi_{iR}=u_{iR},\qquad\psi_{iL}=Q_{iL},\qquad\psi_{4R}=u_{4R},\qquad\psi_{4L}=Q_{4L},\qquad\widetilde{\psi}_{4R}=\widetilde{Q}_{4R},\quad\,\,\widetilde{\psi}_{4L}=\widetilde{u}_{4L}\,.$
(69)

In Eq. (69) the SU(2)L doublets have the following components:
$Q_{iL}=(u_{iL},d_{iL})^{T}$, $Q_{4L}=(u_{4L},d_{4L})^{T}$ and
$\widetilde{Q}_{4R}=(\widetilde{u}_{4R},\widetilde{d}_{4R})^{T}$. The
corresponding mass matrix with the SARAH sign convention reads

$\mathcal{M}_{u}=\left(\begin{array}[]{c|ccccc}&u_{1R}&u_{2R}&u_{3R}&u_{4R}&\widetilde{u}_{4R}\\\\[2.15277pt]
\hline\cr u_{1L}&0&0&0&0&0\\\\[4.30554pt]
u_{2L}&0&0&0&-y_{24}^{u}\frac{v_{u}}{\sqrt{2}}&0\\\\[4.30554pt]
u_{3L}&0&0&0&-y_{34}^{u}\frac{v_{u}}{\sqrt{2}}&x_{34}^{Q}\frac{v_{\phi}}{\sqrt{2}}\\\\[4.30554pt]
u_{4L}&0&0&-y_{43}^{u}\frac{v_{u}}{\sqrt{2}}&0&-M_{4}^{Q}\\\\[4.30554pt]
\widetilde{u}_{4L}&0&x_{42}^{u}\frac{v_{\phi}}{\sqrt{2}}&x_{43}^{u}\frac{v_{\phi}}{\sqrt{2}}&M_{4}^{u}&0\\\
\end{array}\right)\,,$ (70)

with $x^{Q}_{34}\equiv x^{u}_{34}$. The up-type quark mass matrix
$\mathcal{M}_{u}$ can be diagonalized via the mixing matrices $V_{L}^{u}$ and
$V_{R}^{u}$ as

$V_{L}^{u}\mathcal{M}_{u}V_{R}^{u\dagger}=\mathop{\mathrm{d}iag}\left(0,m_{c},m_{t},M_{U_{1}},M_{U_{2}}\right).$
(71)

The approximate expressions for the eigenvalues $m_{c}$ and $m_{t}$ can be
found in Eq. (4), and for the eigenvalues $M_{U_{1}}$ and $M_{U_{2}}$ in Eqs.
(8) and (9), respectively.

###  A.3 Down-type quarks

The down-type quark lagrangian can be obtained from the generic lagrangian (1)
by taking $H=H_{d}$ and making the following replacements of the fermion
fields,

$\psi_{iR}=d_{iR},\qquad\psi_{iL}=Q_{iL},\qquad\psi_{4R}=d_{4R},\qquad\psi_{4L}=Q_{4L},\qquad\widetilde{\psi}_{4R}=\widetilde{Q}_{4R},\quad\,\,\widetilde{\psi}_{4L}=\widetilde{d}_{4L}\,.$
(72)

The corresponding down-type quark mass matrix with the SARAH sign convention
reads

$\mathcal{M}_{d}=\left(\begin{array}[]{c|ccccc}&d_{1R}&d_{2R}&d_{3R}&d_{4R}&\widetilde{d}_{4R}\\\\[2.15277pt]
\hline\cr d_{1L}&0&0&0&y_{14}^{d}\frac{v_{d}}{\sqrt{2}}&0\\\\[4.30554pt]
d_{2L}&0&0&0&y_{24}^{d}\frac{v_{d}}{\sqrt{2}}&0\\\\[4.30554pt]
d_{3L}&0&0&0&y_{34}^{d}\frac{v_{d}}{\sqrt{2}}&-x_{34}^{Q}\frac{v_{\phi}}{\sqrt{2}}\\\\[4.30554pt]
d_{4L}&0&0&y_{43}^{d}\frac{v_{d}}{\sqrt{2}}&0&M_{4}^{Q}\\\\[4.30554pt]
\widetilde{d}_{4L}&0&x_{42}^{d}\frac{v_{\phi}}{\sqrt{2}}&x_{43}^{d}\frac{v_{\phi}}{\sqrt{2}}&M_{4}^{d}&0\\\
\end{array}\right)\,.$ (73)

Note that, unlike in the case of the up-type quarks and charged leptons, it is
impossible to rotate away the $(1,4)$ element of the matrix $\mathcal{M}_{d}$.
The reason is that the mixing between the SM doublets $Q_{1L}$ and $Q_{2L}$
has already been used in the up-quark sector to rotate away the corresponding
entry of $\mathcal{M}_{u}$ King (2018). As a result, the Yukawa coupling
$y_{14}^{d}$ is present in $\mathcal{M}_{d}$. The down-type quark mass matrix
can be diagonalized by the unitary matrices $V_{L}^{d}$ and $V_{R}^{d}$,

$V_{L}^{d}\mathcal{M}_{d}V_{R}^{d\dagger}=\mathop{\mathrm{d}iag}\left(0,m_{s},m_{b},M_{D_{1}},M_{D_{2}}\right).$
(74)

The approximate formulae for the eigenvalues $m_{s}$ and $m_{b}$ can be found
in Eq. (5), and for the eigenvalues $M_{D_{1}}$ and $M_{D_{2}}$ in Eq. (10).

Incidentally, the presence of the matrix element $y_{14}^{d}v_{d}$ has
important consequences for the phenomenology of the model defined in Table 1.
As it was discussed in Sec. II, the first generation of the SM fermions
remains massless if only one complete VL family is added to the spectrum. On
the other hand, the mixing of the $d$ quark with the strange and bottom quarks
is mediated by $y_{14}^{d}v_{d}$. As a result, the full CKM matrix can be
generated in this setup and one needs to include its elements in the global
fit.

###  A.4 Neutrino sector

Finally, we discuss the neutrino mass matrix which emerges from the particle
content given in Table 1. The corresponding lagrangian can be deduced from Eq.
(1) after the following identification:

$\psi_{iL}=L_{iL},\qquad\psi_{4R}=\nu_{4R},\qquad\psi_{4L}=L_{4L},\qquad\widetilde{\psi}_{4R}=\widetilde{L}_{4R},\qquad\widetilde{\psi}_{4L}=\widetilde{\nu}_{4L},\qquad
H=H_{u}\,,$ (75)

where the SU(2)L doublets $L_{iL}$, $L_{4L}$ and $\widetilde{L}_{4R}$ are
defined in Sec. A.1. Note that since there is no $\nu_{iR}$ field in our
model, the couplings $y^{\nu}_{4j}$ and $x^{\nu}_{4j}$ vanish. On the other
hand, the VL neutrino $\nu_{4R}$ is a singlet under the SM gauge symmetry, so
an extra term with $H_{d}^{\ast}$ replacing $H_{u}$ arises. In the end, the
neutrino lagrangian reads:

$\mathcal{L}_{\mathop{\mathrm{r}en},\nu}^{\mathop{\mathrm{Y}ukawa}}=y_{i4}^{\nu}L_{iL}H_{u}\nu_{4R}+x_{i4}^{L}L_{iL}\phi\widetilde{L}_{4R}+y_{i4}^{\prime\nu}L_{iL}H_{d}^{\ast}\widetilde{\nu}_{4L}+M_{4}^{L}L_{4L}\widetilde{L}_{4R}+M_{4}^{\nu}\widetilde{\nu}_{4L}\nu_{4R}+\mathop{\mathrm{h}.c.}\,.$
(76)

Eq. (76) defines a mixed Majorana-Dirac neutrino sector, which after the EWSB
gives rise to a $7\times 7$ Majorana neutrino mass matrix

$\mathcal{M}_{\nu}=\left(\begin{array}[]{c|ccccccc}&\nu_{1L}&\nu_{2L}&\nu_{3L}&\nu_{4L}&\nu_{4R}&\widetilde{\nu}_{4L}&\widetilde{\nu}_{4R}\\\\[2.15277pt]
\hline\cr\nu_{1L}&0&0&0&0&-y_{14}^{\nu}\frac{v_{u}}{\sqrt{2}}&y_{14}^{\prime\nu}\frac{v_{d}}{\sqrt{2}}&0\\\\[4.30554pt]
\nu_{2L}&0&0&0&0&-y_{24}^{\nu}\frac{v_{u}}{\sqrt{2}}&y_{24}^{\prime\nu}\frac{v_{d}}{\sqrt{2}}&0\\\\[4.30554pt]
\nu_{3L}&0&0&0&0&-y_{34}^{\nu}\frac{v_{u}}{\sqrt{2}}&y_{34}^{\prime\nu}\frac{v_{d}}{\sqrt{2}}&x_{34}^{L}\frac{v_{\phi}}{\sqrt{2}}\\\\[4.30554pt]
\nu_{4L}&0&0&0&0&0&0&M_{4}^{L}\\\\[4.30554pt]
\nu_{4R}&-y_{14}^{\nu}\frac{v_{u}}{\sqrt{2}}&-y_{24}^{\nu}\frac{v_{u}}{\sqrt{2}}&-y_{34}^{\nu}\frac{v_{u}}{\sqrt{2}}&0&0&M_{4}^{\nu}&0\\\\[4.30554pt]
\widetilde{\nu}_{4L}&y_{14}^{\prime\nu}\frac{v_{d}}{\sqrt{2}}&y_{24}^{\prime\nu}\frac{v_{d}}{\sqrt{2}}&y_{34}^{\prime\nu}\frac{v_{d}}{\sqrt{2}}&0&M_{4}^{\nu}&0&0\\\\[4.30554pt]
\widetilde{\nu}_{4R}&0&0&x_{34}^{L}\frac{v_{\phi}}{\sqrt{2}}&M_{4}^{L}&0&0&0\\\
\end{array}\right),$ (77)

where once again we chose to work with the SARAH sign convention.

The neutrino mass matrix is symmetric, it can thus be diagonalized via an
orthogonal mixing matrix $V^{\nu}$,

$V^{\nu}\mathcal{M}_{\nu}V^{\nu\dagger}=\mathop{\mathrm{d}iag}\left(0,m_{\nu_{2}},m_{\nu_{3}},M_{N_{1}},M_{N_{2}},M_{N_{3}},M_{N_{4}}\right).$
(78)

##  Appendix B Scalar mass matrices

In this Appendix we collect the explicit formulae for the scalar mass matrices
derived from the scalar potential (19) under the spontaneous symmetry breaking
conditions (3).

The CP-even scalar mass matrix in the basis
$(\mathop{\mathrm{R}e}H_{u}^{0},\,\mathop{\mathrm{R}e}H_{d}^{0},\,\mathop{\mathrm{R}e}\phi)$
evaluated at the vacuum reads

$\mathbf{M}_{\mathop{\mathrm{C}P-even}}^{2}=\left(\begin{array}[]{ccc}\lambda_{1}v_{u}^{2}-\lambda_{5}\frac{v_{d}v_{\phi}^{2}}{4v_{u}}&\lambda_{3}v_{u}v_{d}+\lambda_{5}\frac{v^{2}_{\phi}}{4}&\lambda_{7}v_{u}v_{\phi}+\lambda_{5}\frac{v_{d}v_{\phi}}{2}\\\
\lambda_{3}v_{u}v_{d}+\lambda_{5}\frac{v^{2}_{\phi}}{4}&\lambda_{2}v_{d}^{2}-\lambda_{5}\frac{v_{u}v_{\phi}^{2}}{4v_{d}}&\lambda_{8}v_{d}v_{\phi}+\lambda_{5}\frac{v_{u}v_{\phi}}{2}\\\
\lambda_{7}v_{u}v_{\phi}+\lambda_{5}\frac{v_{d}v_{\phi}}{2}&\lambda_{8}v_{d}v_{\phi}+\lambda_{5}\frac{v_{u}v_{\phi}}{2}&\lambda_{6}v_{\phi}^{2}\\\
\end{array}\right).$ (79)

The matrix (79) can be diagonalized by an orthogonal matrix $R_{h}$
parameterized with three mixing angles. We will denote them as $\alpha_{12}$
for the $(H_{u},H_{d})$ mixing, $\alpha_{13}$ for the $(H_{u},\phi)$ mixing,
and $\alpha_{23}$ for the $(H_{d},\phi)$ mixing. In this parametrization, the
mixing matrix $R_{h}$ is given by

$R_{h}=\left(\begin{array}[]{ccc}c_{12}c_{13}&s_{12}c_{13}&s_{13}\\\
-s_{12}c_{23}-c_{12}s_{13}s_{23}&c_{12}c_{23}-s_{12}s_{13}s_{23}&c_{13}s_{23}\\\
s_{12}s_{23}-c_{12}s_{13}c_{23}&-c_{12}s_{23}-s_{12}s_{13}c_{23}&c_{13}c_{23}\end{array}\right),$
(80)

with the standard notation $s_{ij}=\sin\alpha_{ij}$ and
$c_{ij}=\cos\alpha_{ij}$.

The elements of the matrix $R_{h}$ determine the couplings of the physical
Higgs bosons with the SM particles. It is convenient to define a reduced
coupling as the ratio between the coupling of the physical Higgs scalar
$h_{i}$ and the corresponding coupling of the SM Higgs,

$c_{h_{i}XX}=\frac{g_{h_{i}XX}}{g_{h_{\mathrm{SM}}XX}}\,,$ (81)

where $X$ stands for the SM fermions and gauge bosons. For the model defined
in Table 1, the reduced couplings to quarks and charged leptons are given by

$c_{h_{i}tt}=\frac{(R_{h})_{i1}}{\sin\beta},\qquad
c_{h_{i}bb}=\frac{(R_{h})_{i2}}{\cos\beta}\qquad
c_{h_{i}\tau\tau}=\frac{(R_{h})_{i2}}{\cos\beta}\,,$ (82)

while the reduced couplings to the EW gauge bosons read

$c_{h_{i}ZZ}=c_{h_{i}WW}=(R_{h})_{i1}\sin\beta+(R_{h})_{i2}\cos\beta\,.$ (83)

In this study we choose to work in the alignment limit, which is defined as a
set of constraints on the quartic couplings $\lambda_{i}$ under which the
lightest CP-even scalar $h_{1}$ has the same tree-level couplings with the SM
particles as the SM Higgs. This means that the reduced couplings to fermions
should be very close to 1,

$\frac{\cos\alpha_{12}\cos\alpha_{13}}{\sin\beta}\approx
1,\qquad\quad\frac{\sin\alpha_{12}\cos\alpha_{13}}{\cos\beta}\approx 1\,.$
(84)

It can be easily verify that Eq. (84) leads to the following conditions on the
CP-even scalars mixing angles,

$\alpha_{12}+\beta=\frac{\pi}{2}+n\,\pi,\qquad\alpha_{13}=2\,n\pi\,,\quad\textrm{with}\,\,n=0,1,2\dots$
(85)

indicating no mixing between the doublet $H_{u}$ and the singlet $\phi$. In
this setting, the two SU(2)L scalar doublets mix with the mixing angle
$\frac{\pi}{2}-\beta$, while the doublet $H_{d}$ mixes with the singlet $\phi$
with the mixing angle $\alpha_{23}$. The CP-even scalars mixing matrix thus
reduces to

$R_{h}^{\text{alignment}}=\left(\begin{array}[]{ccc}s_{\beta}&c_{\beta}&0\\\
-c_{\beta}c_{23}&s_{\beta}c_{23}&s_{23}\\\
c_{\beta}s_{23}&-s_{\beta}s_{23}&c_{23}\end{array}\right)=\left(\begin{array}[]{ccc}1&0&0\\\
0&c_{23}&s_{23}\\\
0&-s_{23}&c_{23}\end{array}\right)\times\left(\begin{array}[]{ccc}s_{\beta}&c_{\beta}&0\\\
-c_{\beta}&s_{\beta}&0\\\ 0&0&1\end{array}\right)\,.$ (86)

The alignment conditions (84) translates into the non-trivial relations
between the scalar potential couplings,

$\displaystyle\lambda_{8}\,c^{2}_{\beta}+\lambda_{7}\,s^{2}_{\beta}+\lambda_{5}\,s_{\beta}c_{\beta}$
$\displaystyle=$ $\displaystyle 0$ (87)
$\displaystyle\lambda_{2}\,c^{2}_{\beta}-\lambda_{1}\,s^{2}_{\beta}-\lambda_{3}(c^{2}_{\beta}-s^{2}_{\beta})$
$\displaystyle=$ $\displaystyle 0\,.$ (88)

One can also express the mixing angle $\alpha_{23}$ in terms of the parameters
of the scalar potential,

$\cos(2\alpha_{23})=-\frac{B_{23}}{\sqrt{4A_{23}^{2}+B_{23}^{2}}}\,,\qquad\sin(2\alpha_{23})=-\frac{A_{23}}{\sqrt{4A_{23}^{2}+B_{23}^{2}}}\,,$
(89)

where we define

$\displaystyle A_{23}$ $\displaystyle=$ $\displaystyle
2\,v\,v_{\phi}s_{\beta}(c_{\beta}\lambda_{5}+2\lambda_{7}s_{\beta})$ (90)
$\displaystyle B_{23}$ $\displaystyle=$
$\displaystyle\lambda_{5}v_{\phi}^{2}+4s_{\beta}c_{\beta}\left(\lambda_{6}v_{\phi}^{2}-(\lambda_{1}-\lambda_{3})v^{2}s^{2}_{\beta}\right)\,.$
(91)

The CP-odd mass matrix in the basis
$\left(\mathop{\mathrm{I}m}H_{u}^{0},\mathop{\mathrm{I}m}H_{d}^{0},\mathop{\mathrm{I}m}\phi\right)$
reads

$\mathbf{M}_{\mathop{\mathrm{C}P-odd}}^{2}=-\lambda_{5}\left(\begin{array}[]{ccc}\frac{v_{d}v_{\phi}^{2}}{4v_{u}}&\frac{v^{2}_{\phi}}{4}&\frac{v_{d}v_{\phi}}{2}\\\
\frac{v^{2}_{\phi}}{4}&\frac{v_{u}v_{\phi}^{2}}{4v_{d}}&\frac{v_{u}v_{\phi}}{2}\\\
\frac{v_{d}v_{\phi}}{2}&\frac{v_{u}v_{\phi}}{2}&v_{u}v_{d}-2\frac{\mu_{\mathop{\mathrm{s}b}}^{2}}{\lambda_{5}}\\\
\end{array}\right)\,.$ (92)

Finally, the charged scalar mass matrix is given by

$\mathbf{M}_{\mathop{\mathrm{C}harged}}^{2}=\left(\begin{array}[]{cc}\lambda_{4}\frac{v_{d}^{2}}{2}-\lambda_{5}\frac{v_{d}v^{2}_{\phi}}{4v_{u}}&\lambda_{4}\frac{v_{u}v_{d}}{2}-\lambda_{5}\frac{v_{\phi}^{2}}{4}\\\
\lambda_{4}\frac{v_{u}v_{d}}{2}-\lambda_{5}\frac{v_{\phi}^{2}}{4}&\lambda_{4}\frac{v_{u}^{2}}{2}-\lambda_{5}\frac{v_{u}v^{2}_{\phi}}{4v_{d}}\end{array}\right).$
(93)

##  Appendix C Derivation of the bounded-from-below conditions

In this Appendix, we derive the scalar potential bounded-from-below conditions
shown in Eq. (36). In doing so, we follow the approach of Ref. Bhattacharyya
and Das (2016) and we extend it to the three-field case.

In order to determine the shape of the scalar potential (19) in the limit of
the large fields, it is enough to investigate the behavior of the quartic
terms,

$\begin{split}V_{4}&=\frac{1}{2}\lambda_{1}(H_{u}^{\dagger}H_{u})^{2}+\frac{1}{2}\lambda_{2}(H_{d}^{\dagger}H_{d})^{2}+\lambda_{3}(H_{u}^{\dagger}H_{u})(H_{d}^{\dagger}H_{d})+\lambda_{4}(H_{u}^{\dagger}H_{d})(H_{d}^{\dagger}H_{u})\\\
&-\frac{1}{2}\lambda_{5}(\epsilon_{ij}H_{u}^{i}H_{d}^{j}\phi^{2}+\mathop{\mathrm{H}.c.})+\frac{1}{2}\lambda_{6}(\phi^{*}\phi)^{2}+\lambda_{7}(\phi^{*}\phi)(H_{u}^{\dagger}H_{u})+\lambda_{8}(\phi^{*}\phi)(H_{d}^{\dagger}H_{d}).\end{split}$
(94)

It is convenient to parameterize each quartic term in the following way,

$\begin{split}a&=H_{u}^{\dagger}H_{u}\\\ b&=H_{d}^{\dagger}H_{d}\\\
c&=\phi^{*}\phi\\\ d&=\mathop{\mathrm{R}e}H_{u}^{\dagger}H_{d}\\\
e&=\mathop{\mathrm{I}m}H_{u}^{\dagger}H_{d}\\\
f&=\mathop{\mathrm{R}e}\epsilon_{ij}H_{u}^{i}H_{d}^{j}\phi^{2}\\\
g&=\mathop{\mathrm{I}m}\epsilon_{ij}H_{u}^{i}H_{d}^{j}\phi^{2}\,.\end{split}$
(95)

To make our results more general, we allow $\lambda_{5}$ to be complex. Note
that $a,b,c\geq 0$ by definition, and

$\begin{split}a\,b&\geq d^{2}+e^{2}\\\ a\,b\,c^{2}&\geq f^{2}+g^{2}\geq
2fg.\end{split}$ (96)

In terms of the new parameters, the scalar potential (94) can be rewritten as

$\begin{split}V_{4}&=\frac{1}{4}\left(\sqrt{\lambda_{1}}a-\sqrt{\lambda_{2}}b\right)^{2}+\left(\frac{1}{2}\sqrt{\lambda_{1}\lambda_{2}}+\lambda_{3}\right)a\,b\\\
&+\frac{1}{4}\left(\sqrt{\lambda_{1}}a-\sqrt{\lambda_{6}}c\right)^{2}+\left(\frac{1}{2}\sqrt{\lambda_{1}\lambda_{6}}+\lambda_{7}\right)a\,c\\\
&+\frac{1}{4}\left(\sqrt{\lambda_{2}}b-\sqrt{\lambda_{6}}c\right)^{2}+\left(\frac{1}{2}\sqrt{\lambda_{2}\lambda_{6}}+\lambda_{8}\right)b\,c\\\
&+\lambda_{4}\left(d^{2}+e^{2}\right)-\left(\mathop{\mathrm{R}e}\lambda_{5}f-\mathop{\mathrm{I}m}\lambda_{5}g\right)\,.\end{split}$
(97)

We are now ready to analyze the asymptotic behaviour of the potential (97) in
different field directions.

$\mathbf{a=0}$.

The parameters $d,e,f,g$ automatically vanish, see Eq. (96), and the global
potential reduces to

$V_{4}\left(a=d=e=f=g=0\right)=\frac{1}{2}\left(\sqrt{\lambda_{2}}b-\sqrt{\lambda_{6}}c\right)^{2}+\left(\lambda_{8}+\sqrt{\lambda_{2}\lambda_{6}}\right)bc\,,$
(98)

giving rise to the condition

$\lambda_{8}+\sqrt{\lambda_{2}\lambda_{6}}>0.$ (99)

$\mathbf{b=0}$.

In analogy to the previous case, one obtains

$V_{4}\left(b=d=e=f=g=0\right)=\frac{1}{2}\left(\sqrt{\lambda_{1}}a-\sqrt{\lambda_{6}}c\right)^{2}+\left(\lambda_{7}+\sqrt{\lambda_{1}\lambda_{6}}\right)ac,$
(100)

which gives

$\lambda_{7}+\sqrt{\lambda_{1}\lambda_{6}}>0.$ (101)

$\mathbf{c=0}$.

This time, only the parameters $f$ and $g$ vanish and the reduced scalar
potential reads

$V_{4}\left(c=f=g=0\right)=\frac{1}{2}\left(\sqrt{\lambda_{1}}a-\sqrt{\lambda_{2}}b\right)^{2}+\left(\lambda_{3}+\sqrt{\lambda_{1}\lambda_{2}}\right)ab+\lambda_{4}\left(d^{2}+e^{2}\right)\,.$
(102)

In order to determine the fate of the scalar potential $V_{4}$ at the large
field values, we need to analyze additional directions in the field space. We
first choose a direction along which
$a=\sqrt{\frac{\lambda_{2}}{\lambda_{1}}}b$ and $d=e=0$. Inserting these
expressions into Eq. (102), we arrive to the following condition

$\lambda_{3}+\sqrt{\lambda_{1}\lambda_{2}}>0.$ (103)

Choosing another direction, $a=\sqrt{\frac{\lambda_{2}}{\lambda_{1}}}b$ and
$ab=d^{2}+e^{2}$, we obtain

$\lambda_{3}+\lambda_{4}+\sqrt{\lambda_{1}\lambda_{2}}>0.$ (104)

$\mathbf{a=\sqrt{\frac{\lambda_{6}}{\lambda_{1}}}c}$,
$\mathbf{b=\sqrt{\frac{\lambda_{6}}{\lambda_{2}}}c}$.

Under this assumption the scalar potential (97) reduces to

$\begin{split}V_{4}&=\lambda_{a}c^{2}+\lambda_{4}\left(d^{2}+e^{2}\right)-\left(\mathop{\mathrm{R}e}\lambda_{5}f-\mathop{\mathrm{I}m}\lambda_{5}g\right),\end{split}$
(105)

where one defines

$\lambda_{a}=\frac{3}{2}\lambda_{6}+\lambda_{3}\frac{\lambda_{6}}{\sqrt{\lambda_{1}\lambda_{2}}}+\lambda_{7}\frac{\lambda_{6}}{\lambda_{1}}+\lambda_{8}\frac{\lambda_{6}}{\lambda_{2}}.$
(106)

From Eq. (96) one has

$c^{2}\geq\frac{f^{2}+g^{2}}{d^{2}+e^{2}}\,,$ (107)

leading to

$\begin{split}V_{4}\geq\lambda_{a}\frac{f^{2}+g^{2}}{d^{2}+e^{2}}+\lambda_{4}\left(d^{2}+e^{2}\right)-\left(\mathop{\mathrm{R}e}\lambda_{5}\,f-\mathop{\mathrm{I}m}\lambda_{5}^{\prime},g\right).\end{split}$
(108)

The r.h.s. of Eq. (108) can now be rewritten as

$\textrm{R.H.S}=\left(f-\frac{\mathop{\mathrm{R}e}\lambda_{5}}{c_{1}}\right)^{2}+\left(g+\frac{\mathop{\mathrm{I}m}\lambda_{5}}{c_{1}}\right)^{2}-\frac{1}{4c_{1}}\left((\mathop{\mathrm{R}e}\lambda_{5})^{2}+(\mathop{\mathrm{I}m}\lambda_{5})^{2}\right)+\lambda_{4}\left(d^{2}+e^{2}\right)\,,$
(109)

where

$c_{1}=\frac{\lambda_{a}}{d^{2}+e^{2}}.$ (110)

Choosing an additional direction in the field space,
$f=\frac{\mathop{\mathrm{R}e}\lambda_{5}}{c_{1}}$ and
$g=-\frac{\mathop{\mathrm{I}m}\lambda_{5}}{c_{1}}$, we can derive the
following condition,

$-\frac{1}{4}\frac{(\mathop{\mathrm{R}e}\lambda_{5})^{2}+(\mathop{\mathrm{I}m}\lambda_{5})^{2}}{\lambda_{a}}+\lambda_{4}>0\,.$
(111)

Finally, let us rewrite the r.h.s. of Eq. (108) in yet another way,

$\textrm{R.H.S}=\left(\frac{\sqrt{c_{2}}}{\sqrt{d^{2}+e^{2}}}-\sqrt{\lambda_{4}\left(d^{2}+e^{2}\right)}\right)^{2}+2\sqrt{c_{2}\lambda_{4}}-\left(\mathop{\mathrm{R}e}\lambda_{5}f-\mathop{\mathrm{I}m}\lambda_{5}g\right)\,,$
(112)

where

$c_{2}=\lambda_{a}\left(f^{2}+g^{2}\right),\quad\lambda_{b}=\sqrt{\lambda_{a}\lambda_{4}}\,.$
(113)

Analyzing the quartic potential along the direction
$\sqrt{c_{2}}=\sqrt{\lambda_{4}}(d^{2}+f^{2})$, we obtain

$V_{4}\geq\left(4\lambda_{b}^{2}-(\mathop{\mathrm{R}e}\lambda_{5})^{2}+\mathop{\mathrm{R}e}\lambda_{5}\mathop{\mathrm{I}m}\lambda_{5}\right)f^{2}+\left(4\lambda_{b}^{2}-(\mathop{\mathrm{I}m}\lambda_{5})^{2}+\mathop{\mathrm{R}e}\lambda_{5}\mathop{\mathrm{I}m}\lambda_{5}\right)g^{2}\,,$
(114)

leading straightforwardly to the last two conditions,

$\begin{split}4\lambda_{b}^{2}-(\mathop{\mathrm{R}e}\lambda_{5})^{2}+\mathop{\mathrm{R}e}\lambda_{5}\mathop{\mathrm{I}m}\lambda_{5}&>0\\\
4\lambda_{b}^{2}-(\mathop{\mathrm{I}m}\lambda_{5})^{2}+\mathop{\mathrm{R}e}\lambda_{5}\mathop{\mathrm{I}m}\lambda_{5}&>0.\end{split}$
(115)

##  Appendix D Renormalization group equations

In this Appendix, we collect the one-loop RGEs of our model computed with
SARAH Staub (2014, 2015). We denote

$\beta(X)\equiv\mu\frac{dX}{d\mu}\equiv\frac{1}{16\pi^{2}}\beta^{(1)}(X)\,.$
(116)

$\displaystyle\beta^{(1)}(g_{1})$ $\displaystyle=\frac{103g_{1}^{3}}{15}$
(117) $\displaystyle\beta^{(1)}(g_{2})$ $\displaystyle=-\frac{g_{2}^{3}}{3}$
(118) $\displaystyle\beta^{(1)}(g_{3})$ $\displaystyle=-\frac{13g_{3}^{3}}{3}$
(119) $\displaystyle\beta^{(1)}(\lambda_{1})$
$\displaystyle=-\frac{9}{5}g_{1}^{2}\lambda_{1}-9g_{2}^{2}\lambda_{1}+\frac{27g_{1}^{4}}{100}+\frac{9g_{2}^{4}}{4}+\frac{9}{10}g_{1}^{2}g_{2}^{2}+12\lambda_{1}^{2}+4\lambda_{3}^{2}+2\lambda_{4}^{2}+2\lambda_{7}^{2}+4\lambda_{3}\lambda_{4}$
$\displaystyle+12\lambda_{1}(y_{43}^{u})^{2}+12\lambda_{1}\left[(y_{24}^{u})^{2}+(y_{34}^{u})^{2}\right]-12(y_{43}^{u})^{4}-12\left[(y_{24}^{u})^{2}+(y_{34}^{u})^{2}\right]{}^{2}$
$\displaystyle+4\lambda_{1}\left[(y_{14}^{\nu})^{2}+(y_{24}^{\nu})^{2}+(y_{34}^{\nu})^{2}\right]-4\left[(y_{14}^{\nu})^{2}+(y_{24}^{\nu})^{2}+(y_{34}^{\nu})^{2}\right]^{2}$
(120) $\displaystyle\beta^{(1)}(\lambda_{2})$
$\displaystyle=-\frac{9}{5}g_{1}^{2}\lambda_{2}-9g_{2}^{2}\lambda_{2}+\frac{27g_{1}^{4}}{100}+\frac{9g_{2}^{4}}{4}+\frac{9}{10}g_{1}^{2}g_{2}^{2}+12\lambda_{2}^{2}+4\lambda_{3}^{2}+2\lambda_{4}^{2}+2\lambda_{8}^{2}+4\lambda_{3}\lambda_{4}$
$\displaystyle+12\lambda_{2}(y_{43}^{d})^{2}+12\lambda_{2}\left[(y_{14}^{d})^{2}+(y_{24}^{d})^{2}+(y_{34}^{d})^{2}\right]-12(y_{43}^{d})^{4}-12\left[(y_{14}^{d})^{2}+(y_{24}^{d})^{2}+(y_{34}^{d})^{2}\right]{}^{2}$
$\displaystyle+4\lambda_{2}(y_{43}^{e})^{2}+4\lambda_{2}\left[(y_{24}^{e})^{2}+(y_{34}^{e})^{2}\right]-4(y_{43}^{e})^{4}-4\left[(y_{24}^{e})^{2}+(y_{34}^{e})^{2}\right]{}^{2}$
$\displaystyle+4\lambda_{2}\left[(x_{14}^{\nu})^{2}+(x_{24}^{\nu})^{2}+(x_{34}^{\nu})^{2}\right]-4\left[(x_{14}^{\nu})^{2}+(x_{24}^{\nu})^{2}+(x_{34}^{\nu})^{2}\right]{}^{2}$
(121) $\displaystyle\beta^{(1)}(\lambda_{3})$
$\displaystyle=-\frac{9}{5}g_{1}^{2}\lambda_{3}-9g_{2}^{2}\lambda_{3}+\frac{27g_{1}^{4}}{100}+\frac{9g_{2}^{4}}{4}+\frac{9}{10}g_{1}^{2}g_{2}^{2}+4\lambda_{3}^{2}+2\lambda_{4}^{2}+\lambda_{5}^{2}+6\lambda_{1}\lambda_{3}+6\lambda_{2}\lambda_{3}+2\lambda_{1}\lambda_{4}$
$\displaystyle+2\lambda_{2}\lambda_{4}+2\lambda_{7}\lambda_{8}+6\lambda_{3}(y_{43}^{d})^{2}+6\lambda_{3}\left[(y_{14}^{d})^{2}+(y_{24}^{d})^{2}+(y_{34}^{d})^{2}\right]+2\lambda_{3}\left[(y_{24}^{e})^{2}+(y_{34}^{e})^{2}\right]$
$\displaystyle+2\lambda_{3}(y_{43}^{e})^{2}+6\lambda_{3}(y_{43}^{u})^{2}+6\lambda_{3}\left[(y_{24}^{u})^{2}+(y_{34}^{u})^{2}\right]+2\lambda_{3}\left[(x_{14}^{\nu})^{2}+(x_{24}^{\nu})^{2}+(x_{34}^{\nu})^{2}\right]$
$\displaystyle-4\left[x_{14}^{\nu}y_{14}^{\nu}+x_{24}^{\nu}y_{24}^{\nu}+x_{34}^{\nu}y_{34}^{\nu}\right]{}^{2}+2\lambda_{3}\left[(y_{14}^{\nu})^{2}+(y_{24}^{\nu})^{2}+(y_{34}^{\nu})^{2}\right]$
(122) $\displaystyle\beta^{(1)}(\lambda_{4})$
$\displaystyle=-\frac{9}{5}g_{1}^{2}\lambda_{4}-9g_{2}^{2}\lambda_{4}-\frac{9}{5}g_{1}^{2}g_{2}^{2}+4\lambda_{4}^{2}-\lambda_{5}^{2}+2\lambda_{1}\lambda_{4}+2\lambda_{2}\lambda_{4}+8\lambda_{3}\lambda_{4}+6\lambda_{4}(y_{43}^{u})^{2}+6\lambda_{4}(y_{43}^{d})^{2}$
$\displaystyle+6\lambda_{4}\left[(y_{24}^{u})^{2}+(y_{34}^{u})^{2}\right]-12(y_{43}^{d})^{2}(y_{43}^{u})^{2}-12\left[y_{24}^{d}y_{24}^{u}+y_{34}^{d}y_{34}^{u}\right]{}^{2}+6\lambda_{4}\left[(y_{14}^{d})^{2}+(y_{24}^{d})^{2}+(y_{34}^{d})^{2}\right]$
$\displaystyle+2\lambda_{4}(y_{43}^{e})^{2}+2\lambda_{4}\left[(y_{24}^{e})^{2}+(y_{34}^{e})^{2}\right]-4\left[y_{24}^{e}y_{24}^{\nu}+y_{34}^{e}y_{34}^{\nu}\right]{}^{2}+2\lambda_{4}\left[(x_{14}^{\nu})^{2}+(x_{24}^{\nu})^{2}+(x_{34}^{\nu})^{2}\right]$
$\displaystyle+4\left[x_{14}^{\nu}y_{14}^{\nu}+x_{24}^{\nu}y_{24}^{\nu}+x_{34}^{\nu}y_{34}^{\nu}\right]{}^{2}+2\lambda_{4}\left[(y_{14}^{\nu})^{2}+(y_{24}^{\nu})^{2}+(y_{34}^{\nu})^{2}\right]$
(123) $\displaystyle\beta^{(1)}(\lambda_{5})$
$\displaystyle=\lambda_{5}\Big{\\{}-\frac{9}{20}g_{1}^{2}-\frac{9}{4}g_{2}^{2}+\lambda_{3}-\lambda_{4}+\lambda_{6}+2\lambda_{7}+2\lambda_{8}+6(x_{34}^{Q})^{2}+3\left[(x_{42}^{u})^{2}+(x_{43}^{u})^{2}\right]$
$\displaystyle+\frac{3}{2}(y_{43}^{u})^{2}+\frac{3}{2}\left[(y_{24}^{u})^{2}+(y_{34}^{u})^{2}\right]+3\left[(x_{42}^{d})^{2}+(x_{43}^{d})^{2}\right]+\frac{3}{2}(y_{43}^{d})^{2}+\frac{3}{2}\left[(y_{14}^{d})^{2}+(y_{24}^{d})^{2}+(y_{34}^{d})^{2}\right]$
$\displaystyle+2(x_{34}^{L})^{2}+\left[(x_{42}^{e})^{2}+(x_{43}^{e})^{2}\right]+\frac{1}{2}(y_{43}^{e})^{2}+\frac{1}{2}\left[(y_{24}^{e})^{2}+(y_{34}^{e})^{2}\right]+\frac{1}{2}\left[(x_{14}^{\nu})^{2}+(x_{24}^{\nu})^{2}+(x_{34}^{\nu})^{2}\right]$
$\displaystyle+\frac{1}{2}\left[(y_{14}^{\nu})^{2}+(y_{24}^{\nu})^{2}+(y_{34}^{\nu})^{2}\right]\Big{\\}}$
(124) $\displaystyle\beta^{(1)}(\lambda_{6})$
$\displaystyle=2\lambda_{5}^{2}+10\lambda_{6}^{2}+4\lambda_{7}^{2}+4\lambda_{8}^{2}+24\lambda_{6}(x_{34}^{Q})^{2}-24(x_{34}^{Q})^{4}+12\lambda_{6}\left[(x_{42}^{u})^{2}+(x_{43}^{u})^{2}\right]$
$\displaystyle-12\left[(x_{42}^{u})^{2}+(x_{43}^{u})^{2}\right]{}^{2}+12\lambda_{6}\left[(x_{42}^{d})^{2}+(x_{43}^{d})^{2}\right]-12\left[(x_{42}^{d})^{2}+(x_{43}^{d})^{2}\right]{}^{2}$
$\displaystyle+8\lambda_{6}(x_{34}^{L})^{2}-8(x_{34}^{L})^{4}+4\lambda_{6}\left[(x_{42}^{e})^{2}+(x_{43}^{e})^{2}\right]-4\left[(x_{42}^{e})^{2}+(x_{43}^{e})^{2}\right]{}^{2}$
(125) $\displaystyle\beta^{(1)}(\lambda_{7})$
$\displaystyle=-\frac{9}{10}g_{1}^{2}\lambda_{7}-\frac{9}{2}g_{2}^{2}\lambda_{7}+2\lambda_{5}^{2}+4\lambda_{7}^{2}+6\lambda_{1}\lambda_{7}+4\lambda_{6}\lambda_{7}+4\lambda_{3}\lambda_{8}+2\lambda_{4}\lambda_{8}$
$\displaystyle-12(x_{34}^{Q})^{2}(y_{34}^{u})^{2}+12\lambda_{7}(x_{34}^{Q})^{2}+6\lambda_{7}\left[(x_{42}^{u})^{2}+(x_{43}^{u})^{2}\right]-12(x_{43}^{u})^{2}(y_{43}^{u})^{2}+6\lambda_{7}(y_{43}^{u})^{2}$
$\displaystyle+6\lambda_{7}\left[(y_{24}^{u})^{2}+(y_{34}^{u})^{2}\right]+6\lambda_{7}\left[(x_{42}^{d})^{2}+(x_{43}^{d})^{2}\right]+4\lambda_{7}(x_{34}^{L})^{2}-4(x_{34}^{L})^{2}(y_{34}^{\nu})^{2}$
$\displaystyle+2\lambda_{7}\left[(x_{42}^{e})^{2}+(x_{43}^{e})^{2}\right]+2\lambda_{7}\left[(y_{14}^{\nu})^{2}+(y_{24}^{\nu})^{2}+(y_{34}^{\nu})^{2}\right]$
(126) $\displaystyle\beta^{(1)}(\lambda_{8})$
$\displaystyle=-\frac{9}{10}g_{1}^{2}\lambda_{8}-\frac{9}{2}g_{2}^{2}\lambda_{8}+2\lambda_{5}^{2}+4\lambda_{8}^{2}+4\lambda_{3}\lambda_{7}+2\lambda_{4}\lambda_{7}+6\lambda_{2}\lambda_{8}+4\lambda_{6}\lambda_{8}-12(y_{34}^{d})^{2}(x_{34}^{Q})^{2}$