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

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.


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 [1] and extra dimensions [2][3][4] 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 [5][6][7][8][9][10][11], 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. [12] 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. [12] 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.[12,13], Zmediated Flavour Changing Neutral Currents in Ref. [14], and a deviation from the SM prediction in the measured value of the anomalous magnetic moment of muon in Refs.[13,15,16].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. [16] were characterized by relatively low (∼ 200 GeV) masses of the VL lepton doublets and large (∼ 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. [16] 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.[13][14][15][16].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. [16] it was required that all the dimensionless parameters of the lagrangian remain perturbative (in a loose sense of being smaller than √ 4π for the gauge/Yukawa and smaller than 4π 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. [16] were disfavored.
Last but not least, we refine the derivation of the stability conditions for the scalar potential which in Refs.[13,16] 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.[13,16], 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 τ 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. [12].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.[12][13][14]16].
The particle content of the model is summarized in Table I.The SM fermion sector, collectively denoted as ψ i (ψ 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 (ψ 4 , ψ 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 ϕ.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.
All the renormalizable Yukawa interactions between the SM and NP fermions which are allowed by the extended gauge symmetry can be schematically written as: where H is either H u or H d and M ψ 4 (M ψ 4 ) denotes the VL doublet (singlet) mass parameter.Note that with the U(1) X charges given in Table I 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.

A. Hierarchy of masses
Once the neutral components of the scalar fields develop their vevs, the 5 × 5 fermions mass matrices are generated.Since their upper 3 × 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 [12] that this allows one to choose a flavor basis in which the fermion mass matrices acquire the following form: In the above, the term in parentheses assumes a non-zero value in the mass matrix of the downtype quarks, while it is zero for the up-type quarks and charged leptons.The exact forms of the matrices M u , M d and M e are presented in Appendix A.
In order to calculate the masses of the physical quarks and leptons, the 5 × 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 M ψ on the couplings y ψ i4 , y ψ 43 , x ψ 4i , x ψ 34 and the masses M ψ( ψ) 4 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 and defining tan β = v u /v d , the masses of the third and second family quarks and leptons are approximately given by (see also Ref. [12] for a related derivation) 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 ≪ 1 or x u 42 ≪ 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 ϕ should be of the same order as the VL mass parameter M Q 4 .Note, however, that in our phenomenological analysis we always perform the numerical diagonalization of the mass matrices (A2), (A5) and (A8).
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 ≈ 34, puts relevant constraints on the allowed parameter space of the model.In fact, we have The relation (7) leads to two distinct classes of solutions.In the first one, with both the Yukawa couplings of order one, tan β ∼ O(10) is required.In the other one, with tan β ∼ O(1), a large hierarchy between the up and down sector couplings, y d 43 ≪ 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, In the neutrino sector, the corresponding mass matrix is 7 × 7 and its explicit form can be found in Eq. (A12).The resulting masses of the heavy neutrinos read 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. [16]).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.

B. CKM mixing matrix
The full 5 × 5 mixing matrix takes the following form [14]: where V u L and V d L are the left-handed mixing matrices of Eqs.(A6) and (A9) which diagonalize the up-and down-type quark mass matrices M u and 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 ± .Following the strategy of Ref. [14] and working under the assumption that v u,d /M Q,u,d 4 ≪ 1, we can approximate the 3 × 3 CKM matrix as where Based on the conclusions from Sec. II A one expects x u , x d ≪ 1.Note also that: • The element V us of the CKM matrix is given by y d 14 /y d 24 in our model.The presence of a non-zero coupling y d 14 is thus crucial to generate the Cabibbo angle of the right size.We also expect y d 14 ≈ 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 ≈ 0.017.It then follows that x u ≈ −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 [17].Defining, for example, one obtains Plugging the Wolfenstein parameters extracted from the global fit [18] into Eq.( 17) and comparing it with the experimental determination of the CKM matrix elements reported in Ref. [18], one can estimate to what extent the measured structure of the CKM matrix can be reproduced in our model.One obtains 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). 1 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.
A. Scalar masses in the alignment limit In the interaction basis, the most generic renormalizable scalar potential of the model defined in Table I takes the form [16]: where µ 2 u,d,ϕ are dimensionful mass parameters, λ 1,2,••• ,8 denote dimensionless quartic coupling constants, and µ 2 sb 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 ϕ around their vacuum states, 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, One must have 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, Re H 0 u , Re H 0 d and Re ϕ, account for three CP-even Higgs bosons.The corresponding mass matrix M 2 CP−even (see Eq. (B1)) can be diagonalized with a mixing matrix R h defined in Eq. (B2).The masses of three physical neutral scalars, h 1 , h 2 and h 3 , correspond to the eigenvalues of M 2 CP−even , 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 λ 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 where the equality imposes a perfect alignment condition.The masses of the CP-even scalars in the alignment limit read with A 23 and B 23 defined in Eq. (B12) and Eq.(B13), respectively, and The CP-odd scalar mass matrix in the basis Im H 0 u , Im H 0 d , Im ϕ , M 2 CP−odd , is defined in Eq. (B14).After the diagonalization, the physical CP-odd spectrum consists of one massless Goldstone boson and two massive pseudoscalars, a 1 and a with the masses given by Note that λ 5 < 0 and µ 2 sb > 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 Charged , is defined in Eq. (B15).After the diagonalization with a mixing matrix R β , one is left with a massless charged Goldstone boson and a charged Higgs boson, The corresponding mass reads in this case As a closing remark, let us notice that the alignment condition (25) indicates In order to preserve the perturbativity of λ 2 (more on this in Sec.V), the term in parentheses needs to be fine-tuned with a precision O(1/ tan 2 β) or better, effectively fixing λ 3 ≈ λ 1 with the same accuracy.On the other hand, Eq. ( 26) implies that we can identify λ 1 with the quartic coupling of the SM, λ 1 = 0.258, as long as tan β ≳ 3. Similarly, the alignment condition (24) gives Perturbativity of λ 8 then requires λ 7 ∼ O(1/ tan 2 β) and λ 5 ∼ O(1/ tan β).

B. 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 λ i .To derive the 'bounded-from-below' constraints, one should analyze all possible directions along which the scalar fields H u , H d and ϕ 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): where Since in this study we do not investigate the CP violation, we assume that all the parameters of the lagrangian are real, indicating Im λ 5 = 0. Note also that several novel conditions w.r.t. the findings of Refs.[13,16] are identified in Eq. (36).

C. 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.[19]).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++ [20] (the C++ version of [21]) to find all treeand 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.

A. Muon anomalous magnetic moment
The discrepancy between the SM prediction  and the experimental measurement of the anomalous magnetic moment of the muon has been confirmed separately by the Brookhaven National Laboratory [44] and the Fermilab experimental groups [45,46], giving rise to the combined 5.1 σ anomaly: In a generic NP model which features heavy scalars ϕ i and fermions ψ j coupled to the SM muons via the Yukawa-type interactions y ij L ϕ i ψj P L µ and y ij R ϕ i ψj P R µ (where P L,R = (1 ∓ γ 5 )/2 are the usual projection operators), a well-known one-loop contribution to the muon anomalous magnetic moment reads where M ϕ i is the physical mass of a heavy scalar, M ψ j is the physical mass of a heavy fermion, , and the electric charges of ϕ i and ψ j are related as The loop functions are defined in the following way: The first addend in Eq. ( 38) captures the loop chirality-conserving contributions to ∆a µ .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 [47,48].We will thus focus on the second addend in Eq. (38), which corresponds to the loop chirality-flipping contributions to ∆a µ .
In the framework of the model defined in Table I, 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 ∆a µ read for the CP-even scalars and for the CP-odd scalars.The one-loop contributions to ∆a µ from the neutral leptons and charged scalars are given by Figure 1: The one-loop chirality-flipping contributions to ∆a µ mediated by (a) a neutral (pseudo)scalar/charged lepton exchange, and (b) a charged scalar/neutral lepton exchange.
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 (A11), as well as on the elements of the mixing matrices R h (Eq.( 23)), R a (Eq.( 29)), and V e L/R (Eq.(A3)).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 ∆a µ with the numerical package SPheno [49,50].

B. 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 τ → µγ decay receives contributions from the one-loop diagrams analogous to those of ∆a µ .The corresponding branching ratio (BR) is given by [51] where Γ τ = 2.3 × 10 −12 [18] indicates the total decay width of the tau, α em is the fine structure constant, and the decay amplitude A ij L reads The corresponding amplitude A ij R is obtained from Eq. ( 44) by replacing L ↔ R. Just like it was in the ∆a µ case, the main contribution to BR(τ → µγ) originates from the second addend in Eq. (44).The current experimental 90% confidence level (C.L.) upper bound on BR (τ → µγ) from the Belle collaboration reads [52]: The τ → 3 µ 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. [51].The 90% C.L. upper bound on BR (τ → 3 µ) by the Belle collaboration reads [53]: C. 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 β decay, from the experimental data on kaon decays, and from the lattice calculations [54][55][56][57].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 V ud and V us resulting from the β decay and from the kaon decays.The second anomaly originates from two different measurements of V us : from the semileptonic K → πlν and the leptonic K → µν decay, respectively.
The experimental upper bound on the CKM deviation from the unitarity reads [18] ∆ 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 [55,[58][59][60].In such a setting deviations from the unitarity of the three-dimensional CKM matrix can arise quite naturally.Since the model defined in Table I 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 I 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 Here Λ 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 = C m µ /Λ.
One can now derive from Eq. ( 48) rough estimates of the energy scale associated with a hypothetical NP contributing to ∆a µ at different loop orders, tree level : C ≈ 1, Λ ≈ 3000 GeV (49) and so on.
Going beyond the EFT approximation, let us investigate a one-loop chirality flipping contribution to ∆a µ like the one in Eq. (38).Assuming that it arises from an unspecified UV completion of our model above the scale Λ, it can be estimated by the corresponding UV mass m Λ and the UV Yukawa couplings y L/R (Λ) as By demanding that the new contribution (52) does not shift our phenomenological predictions for ∆a µ by more than 3 σ, we can derive a lower bound on the UV mass, For the Yukawa couplings at the upper edge of perturbativity, y L (Λ) = y R (Λ) = √ 4π, Eq. ( 53) translates into a conservative estimation of the scale of validity of our phenomenological model, In other words, the model can not be UV completed below m Λ .An immediate consequence of Eq. ( 54) is that the commonly employed perturbativity bounds, which read ≲ √ 4π for the gauge/Yukawa and ≲ 4π for the quartic couplings, need to be imposed on the running parameters of the model evaluated at the scale Λ rather than on the bare couplings of the lagrangian (as it was done, for example, in Ref. [16]).
To implement the RG-based perturbativity constraints, we follow the RG flow of all the coupling constants from the scale µ 0 = 1.5 TeV, which is a proxy for the NP scale in our model, to Λ = 50 TeV.The one-loop RG equations (RGEs) were computed using the publicly available numerical code SARAH [61,62] 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 µ 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 µ 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 µ 0 are set to 0.7.
In the scalar sector, the perturbativity bounds are presumably most relevant for the couplings λ 1 , λ 6 and λ 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. (D5), Eq. (D10) and Eq.(D11), respectively).In Fig. 2 we illustrate the RG running of λ 1 , λ 6 and λ 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. [16] are disfavored.
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 → 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 (χ 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 χ 2 function.Finally, we provide a discussion of some experimental signatures that these benchmark scenarios could produce.

A. Scanning methodology
In Table II 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 β.
In the scalar sector the alignment conditions ( 24) and ( 25) are imposed, leading to the limited scanning ranges for λ 3 , λ 5 and λ 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.Table II: Scanning ranges for the input parameters of the model defined in Table I.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 GeV 2 .
A tentative lower bound of 1200 GeV 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 A. The scanning range for v ϕ then follows from the requirement of reproducing the correct mass of the top quark, as discussed in Sec.II A. 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 B. Finally, the range for µ 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 [18].
The experimental constraints employed in our numerical scan are listed in Table III.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 [18].Since the uncertainties for m µ and m τ 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 χ 2 -statistic function as where O model i indicates the value of an observable calculated in our model, O cen i is the central value of its experimental measurement, O err i is the corresponding experimental error, and the sum runs over all the measured observables listed in Table III.The upper bounds, corresponding to the last three rows of Table III, are not included in the χ 2 function, but applied as hard-cuts instead (a point in the parameter space is rejected if such a condition is not satisfied).
To minimize the χ 2 function, we adopt the following strategy.First, we perform an initial scan of the parameter space consistent with Table II.As a result, we obtain a seed which is then used to minimize the χ 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 κ%, and κ 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 ∆ CKM , BR(τ → µγ) and BR(τ → 3µ), 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++ [20] and we keep only those points whose vacuum is identified as "stable".

B. Benchmark scenarios
In Table IV we present input parameters for three best-fit benchmark scenarios identified by performing the numerical scan discussed in Sec.VI A. The corresponding mass spectra are summarized in Table V while the breakdown of individual contributions to the χ 2 function is shown in Table VI.
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 χ 2 function is smaller than 0.7, with an exception of m s in BP1, in which case we have χ 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 χ 2 stems from this very sector.As anticipated in Sec.II B, the main contribution to the χ 2 function is given by the element |V td |, with the corresponding χ 2 V td ranging from 16 for BP1 to 25 The two remaining best-fit points follow the same pattern.Incidentally, note that the CKM anomaly is 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π ( √ 4π) up to 1000 TeV.We can therefore conclude that the validity range of our model extends well beyond the putative scale of 50 TeV.We also checked that the maximal eigenvalue of the scattering matrix computed by SPheno is O(10 −2 ) for all the benchmark scenarios, indicating that the perturbative unitarity bound is satisfied as well.
Masses of the NP leptons are determined, to a large extent, by correctly fitting the experimental value of ∆a µ (an overall contribution from this observable to the total χ 2 function does not exceed 0.7 in all the benchmark scenarios).Contributions to ∆a µ from the individual one-loop diagrams of Fig. (1) are summarized in Table VII.We present separately fractions of ∆a µ generated by the charged scalars h ± 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".TOT stands for the total χ 2 function of each best-fit scenario.
We observe that the largest contributions to ∆a µ arise from the charged scalar/heavy neutrino loops.We thus disprove the conclusions of Refs.[13,16] 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 ∆a µ 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 Contributions to ∆a µ × 10 9 Charged scalars CP-even scalars  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 ′ν 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 ∆a µ are suppressed.Additionally, one also observes cancellations between the individual contributions to ∆a µ 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., [63,64] for a discussion).Finally, we should mention the size of the BRs for the lepton flavor violating decays τ → µγ and τ → 3µ in our model, which are of the order (3 − 4) × 10 −8 for the former and (6 − 9) × 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 [65,66], 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.

A. 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. [67,68] 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.
The most recent analysis from ATLAS, based on the data from proton-proton collisions at a centre-of-mass energy of √ s = 13 TeV, corresponding to an integrated luminosity of 139 fb −1 [69], considered the pair production of VL top partners T and VL bottom partners B with the decay channels T → Zt, ht, W b and B → Zb, hb, W t and with large missing transverse momentum.The for the EW singlets. 3The analogous results from CMS can be found in Ref. [70].Similarly, by assuming that at least one of the VLQs decays into a Z boson with the BR=100%, the 13 TeV ATLAS search [71] obtained even stronger bounds, 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 [72] the cross sections for the pair production of U 1 , U 2 , D 1 and D 2 .The results are presented in Table VIII.By comparing these numbers with the observed experimental 95% C.L. upper bounds on the signal cross section from Ref. [69] (to give an example, σ exp 95% (p p → T T ) = 4 × 10 −3 pb for M 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.
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 2 For the VLQ mass larger than 800 GeV this indicates BR(T → Zt)=BR(T → ht) = 50% and BR(B → W t) = 100% [69]. 3For the VLQ mass larger than 800 GeV this indicates BR(T : 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.[73][74][75] for T and in Ref. [76] for B. where , the cross section (60) reduces to the signal cross section constrained by the experimental collaborations, σ exp 95% ≈ σ pp → QQ .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: 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.[73][74][75], while the single VL B quark production in Ref. [76].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 O(10 −5 ) pb, while for D 1 it amounts to 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.[77,78].

B. 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 IX: Cross sections (in pb) for the pair production of the VLLs for our three benchmark scenarios.Masses are in GeV.
Table IX.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: BP3 : 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, 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: The dominant branching ratios for the subsequent decays of E 1 and N 3,4 , evaluated with SPheno, are collected in Table X.In all three cases the VLLs decay predominantly to the SM muons, which is a direct consequence of the fact that we impose ∆a µ as a constraint in our likelihood function and the largish muon-lepton-scalar Yukawa couplings are preferred.
A closer look at Table X 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 C).As an example, let us consider a process p p → E 1 Ē1 → µ μ a 1 a 1 for the benchmark scenario BP1.The lightest pseudoscalar can decay in this case either to a b b pair (with the BR of 28%) or to ν N 1,2 (with the BR of 69%).The decay of the heavy neutrinos then proceeds as N 1,2 → e ± /µ ± W ± with the BR of 56%, or N 1,2 → ν Z with the BR of 28%.We can thus expect the following distinctive experimental signatures emerging from the Table X: Dominant BRs for the decays of E 1 and N 3,4 .and CMS [80] collaborations.
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 IX 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 fb −1 of data from the 13 TeV run [79], looks for the VLLs coupled predominantly to taus.The analogous CMS analysis based on the 77.4 fb −1 of data can be found in Ref. [80].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).
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.[79] and [80] can be written as follows: Combining the cross sections from Table IX with the relevant BRs calculated with SPheno, we obtain BP3 : σ (p p → τ − τ + l − l + q q) = 5.3 × 10 −7 pb .
If we now compare the predictions of Eq. ( 65) with the corresponding experimental 95% C.L. exclusion cross sections from Ref. [79], we can conclude that the benchmark scenarios identified in Sec.VI are not excluded by the dedicated LHC searches for the VLLs. 4 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.
Table XI: Dominant BRs (> 5%) for the decays of the exotic scalars.

C. 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. [82] 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 XI).
To facilitate the analysis, we use the publicly available code HiggsTools [83], 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 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 obs > 1.We observe that the highest values of R obs (0.6 for BP1 and BP3, 0.14 for BP2) are reached for the h 2 → τ + τ − and a 1 → τ + τ − decays constrained by the ATLAS 139 fb −1 analysis [84], despite very low decay BRs in this channel.
To investigate it in more details, we calculated the a 1 /h 2 → τ + τ − cross sections with MadGraph5 MC@NLO.The results are reported in the last three columns of Table XII.These are to be compared with the 95% C.L. experimental lower bounds on the cross section reported in the third column of Table XII.We find a very good agreement with the output of HiggsTools in terms of the parameter R 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 XII 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 τ − τ + 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, γγ, ZZ Table XII: 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.
and W W , are O(10 −8 ) and 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 → t + t − decay through the measurement of an effective coupling g a 1 /h 2 tt .This scenario was investigated by CMS in Ref. [85].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) 5 we conclude that no additional constraint on our model arises from this particular search.

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. [12] 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.[13,16].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 ∆a µ .This finding is qualitatively different from the conclusions obtained in Refs.[13,16] where only the charged lepton loops were considered.We would like to emphasize that the dominance of the heavy neutrino contribution to ∆a µ 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 → τ + τ − , 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.The mass matrix for the charged leptons, M e , can be derived from Eq. ( 1) after identifying the generic fermions ψ with the corresponding lepton fields from Table I and the generic scalar H with H d .One thus has with the following components of the SU(2) L doublets: L iL = (ν iL , e iL ) T , L 4L = (ν 4L , e 4L ) T and L 4R = ( ν 4R , e 4R ) T .As a result, the mass matrix reads e 4R e 1L 0 0 0 0 0 where M L 4 (M e 4 ) denotes the mass of the VL lepton doublet (singlet) and x L 34 ≡ 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 × 5 charged lepton mass matrix M e can be diagonalized by means of two unitary transformations V e L and V e R , In Sec.II the approximate expressions for the eigenvalues m µ and m τ 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.

Up-type quarks
In analogy to the charged lepton sector, the mass matrix for the up-type quarks, M u , can be derived from Eq. (1) after taking H = H u and making the following identification: In Eq. (A4) the SU(2) L doublets have the following components: The corresponding mass matrix with the SARAH sign con-vention reads The up-type quark mass matrix M u can be diagonalized via the mixing matrices V u L and V u R as 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.

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, The corresponding down-type quark mass matrix with the SARAH sign convention reads 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 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 M u [12].As a result, the Yukawa coupling y d 14 is present in M d .The down-type quark mass matrix can be diagonalized by the unitary matrices V d L and V d R , 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 d 14 v d has important consequences for the phenomenology of the model defined in Table I.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 d 14 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.(B1) The matrix (B1) can be diagonalized by an orthogonal matrix R h parameterized with three mixing angles.We will denote them as α 12 for the (H u , H d ) mixing, α 13 for the (H u , ϕ) mixing, and α 23 for the (H d , ϕ) mixing.In this parametrization, the mixing matrix R h is given by 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, where X stands for the SM fermions and gauge bosons.For the model defined in Table I, the reduced couplings to quarks and charged leptons are given by while the reduced couplings to the EW gauge bosons read In this study we choose to work in the alignment limit, which is defined as a set of constraints on the quartic couplings λ 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, cos α 12 cos α 13 sin β ≈ 1, sin α 12 cos α 13 cos β ≈ 1 .

(B6)
It can be easily verify that Eq. (B6) leads to the following conditions on the CP-even scalars mixing angles, 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. [93] 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, (C1) It is convenient to parameterize each quartic term in the following way, (C2) To make our results more general, we allow λ 5 to be complex.Note that a, b, c ≥ 0 by definition, and In terms of the new parameters, the scalar potential (C1) can be rewritten as where one defines

Figure 2 :
Figure 2: The RG running of the quartic couplings λ 1 , λ 6 and λ 7 for a randomly chosen benchmark point which satisfies all the constraints discussed in Secs.II and III.The renormalization scale µ ranges from 1.5 TeV to 1000 TeV.µ 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.

Figure 3 :
Figure 3: Pair production of the VLQs T and B via the gluon fusion at the LHC considered by the ATLAS collaboration in Ref. [69].

Figure 5 :
Figure5: Pair production of the VLLs E 1 and N 3,4 at the LHC, as considered by the ATLAS[79] and CMS[80] collaborations.

Table I :
The particle content of the NP model considered in this study.
Measurement Central Value Exp.Error Measurement Central Value Exp.Error Table III:The experimental measurements which we employ in our numerical scan.Masses are in GeV.

Table V :
for BP3.Smaller yet still relevant contributions come from the entries |V ub | and |V ts |.Finally, an order 10 contribution to the χ 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 Mass spectra for three best-fit benchmark scenarios.All masses are in GeV.

Table VI :
Breakdown of the χ 2 contributions from various observables implemented in the χ 2 function of Eq. (55).The CKM contributions which are smaller than 3 are not shown.χ 2 Q , χ 2 L and χ 2 V indicate total χ 2 contributions from the quark masses, lepton masses, and the CKM matrix elements, respectively.χ 2

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