A three Higgs doublet model with symmetry-suppressed flavour changing neutral currents

We construct a three-Higgs doublet model with a flavour non-universal ${\rm U}(1)\times \mathbb{Z}_2$ symmetry. That symmetry induces suppressed flavour-changing interactions mediated by neutral scalars. New scalars with masses below the TeV scale can still successfully negotiate the constraints arising from flavour data. Such a model can thus encourage direct searches for extra Higgs bosons in the future collider experiments, and includes a non-trivial flavour structure.


Introduction
The properties of the new resonance observed at the LHC in 2012 [1,2] seem tantalizingly close to those of the Higgs boson predicted by the Standard Model (SM) (for instance, see [3,4]). The particle spectrum predicted by the SM has now been fully confirmed, but many important questions in particle physics are left unanswered: the smallness of neutrino masses, the fermion mass hierarchy, the colossal asymmetry between the quantity of matter and antimatter in the universe and the nature of dark matter are a few of such unresolved issues. They are usually taken as hints for the existence of new physics (NP) beyond the SM (BSM). Typical BSM scenarios that aim to fix one or more such shortcomings of the SM often end up extending the scalar sector of the SM. In these extensions, the 125 GeV scalar observed at the LHC is not the only scalar in the spectrum but the first one in a series of others to follow. This is an intriguing possibility which motivates us for a closer inspection of the properties of the observed scalar and inspires us to carry on our efforts to look for new resonances at collider experiments.
When it comes to extending the scalar sector of the SM, adding replicas of the SM Higgs-doublet is one of the simplest ways to do it. Such extensions do not alter the tree level value of the electroweak (EW) ρ-parameter. Two Higgs-doublet models (2HDMs), which add only one extra doublet to the SM Higgs sector, have received its fair share of attention through the years. They were proposed by T.D. Lee in 1973 [5] as a means to obtain a spontaneous breaking of the CP symmetry, and boast a rich phenomenology. Other than the possibility of spontaneous CP breaking, such models contain a richer particle spectrum, with a charged scalar and a total of three neutral ones, may feature dark matter candidates in certain scenarios, and generically give rise to the treelevel scalar-mediated flavour changing neutral currents (FCNCs). Indeed, one ominous outcome of adding extra scalar doublets is that the fermions of a particular charge will now receive their masses from multiple Yukawa matrices and consequently, diagonalization of their mass matrices will no longer guarantee the simultaneous diagonalization of the Yukawa matrices. Therefore, in general, there will exist FCNCs at tree-level mediated by neutral scalars.
Experimental data from the flavour sector -for instance, neutral meson mass differences for Kaons and Bmesons, or K data -strongly constrain such FCNCs, typically forcing the extra scalars to have masses above 1 TeV [6]. An alternative is to fine-tune the FCNC interactions so that they are very small, though for some models cancellations between CP-even and CP-odd contributions to the amplitudes off the observables mentioned allow for below-TeV scalars with a minimal fine-tuning (see, for instance, [7][8][9]). Another possibility is to assume alignment between different Yukawa matrices [10][11][12][13][14][15], though that is an ansatz which is not preserved under renormalization [16]. There is yet another possibility, however: the BGL (Branco-Grimus-Lavoura) model [17,18] is a remarkable version of the 2HDM wherein FCNC interactions are naturally smallthis results from the imposition of a flavour-violating symmetry, i.e. a symmetry which treats some generations of fermions differently from others. As a consequence of this symmetry, FCNC couplings are suppressed by off-diagonal Cabibbo-Kobayashi-Maskawa (CKM) matrix elements. The phenomenology of the model has been studied quite thoroughly (see, for instance, [19][20][21]) and it remains a valid and exciting possibility for BSM physics.
As a next-to-minimal possibility along these lines, recent years have seen a growing interest in the topic of three Higgs-doublet models (3HDMs). These models add two extra scalar doublets on top of the SM one, thereby conforming to the aesthetic appeal of having three scalar generations in harmony with the three generations of fermions. An early proposal by Weinberg [22] included two discrete symmetries which yielded flavour conservation. A recent study [23], following earlier work in [24][25][26], used a generalized CP symmetry (called CP4) to constrain the vast parameter space of 3HDM, and was found to have large regions of parameter space which conformed with experimental constraints in the flavour and scalar sectors. That model, however, used the liberty in parameters present in the Yukawa sector to fit all observables, but the FCNC couplings were not necessarily small. Furthermore, it was shown in [27] that the model contains several regions where top decays to light charged Higgs bosons is in conflict with LHC data. In [28,29] is was also studied how baryogenesis can be realized via the decay of new TeV-scale Higgs bosons in a 3HDM. Last but not least, in [30], a maximally-symmetric version of the 3HDM based on the Sp(6) symmetry was proposed.
In this article, we will attempt to implement the BGL method in a 3HDM framework. We construct a 3HDM with a U(1) × Z 2 flavor symmetry, where we also achieve the smallness of the FCNC couplings in a way similar to the BGL case. In the considered version of the model, tree-level FCNC interactions occur in the down-quark sector while those in the up-quark sector are forbidden by the flavor symmetry. In this work, we employ publicly available tools such as the generic BSM spectrum generator SARAH/SPheno interfaced with the widely-used Higgs (HiggsBounds/HiggsSignals) and flavor (Flavio) observables' analysers enabling us to thoroughly verify the model parameter space against the most relevant theoretical and experiments bounds. We find that relatively light scalars can successfully pass through the stringent experimental constraints arising from flavor data and hence may occur within the reach of the current and future collider experiments. To motivate future searches, we also outline possible signatures of nonstandard scalars present in our model. For such a purpose we use MadGraph5 aMC@NLO 2.6.2 and focus on scalar and pseudoscalar production via gluon fusion with subsequent decay into tau leptons.
Our article is organized as follows. In Sec. 2 we review the framework of BGL models. Then, in Sec. 3, we build our model, a 3HDM endowed with a U(1) × Z 2 symmetry with a non-trivial structure in the Yukawa sector. In Sec. 4 we review the constraints imposed upon the model, both theoretical -boundedness from below, unitarity, electroweak precision bounds -and experimental -LHC Higgs data and searches for heavier scalars, flavour physics data, among others. In Sec. 5 we explain in detail the procedure we followed to perform a thorough numerical scan of the model and present the results we found for the parameter space that survives all constraints imposed upon the model. We conclude in Sec. 6 with an overview of this work and a discussion of its significance.

The BGL model
The BGL model is a version of the 2HDM where the scalar interactions with fermions violate flavour -meaning, unlike the interactions of the photon and Z boson, the neutral scalars in the BGL model "jump families" like the W boson does. In 2HDM the general recipe to avoid FCNC is to only allow fermions of the same charge to couple to just one of the doublets [31]. This is usually enforced by imposing a Z 2 or U(1) [32] symmetry on the model (see also [33,34]). The reason for doing this in the first place is the fact that tree-level mediated FCNCs would make significant contributions to flavour sector observables such as the mass differences of the K 0 , B d or B s mesons, or to the K quantity, ruining an agreement found within the SM for those quantities -unless the masses of the new scalars are all of order TeV, or the FCNC Yukawa couplings are tuned to be very small.
The BGL model is remarkable since it forces the FCNC couplings to be heavily suppressed as the result of a symmetry. The model therefore provides a simple and natural explanation as to why NP contributions to flavour observables would not ruin the agreement found within the SM, without the need of any fine-tuning. To understand how this is achieved, consider the Yukawa Lagrangian for the quark sector in the 2HDM, where Q a L = (p a L n a L ) T and φ i are the SU(2) weak isospin (left-handed) quark and Higgs doublets, respectively, andφ k = iσ 2 φ * k . The p and n fields are the positively and negatively charged quark fields, respectively. Upon rotation to the mass basis they will yield the physical up and down quarks. The a and b are fermion family indices. Γ and ∆ are 3×3 Yukawa coupling matrices for the down and up sector, respectively. Upon spontaneous symmetry breaking, the scalar doublets develop neutral vacuum expectation values (VEVs), such that 1 We define tan β = v 2 /v 1 . For a CP-conserving model (both at the explicit and vacuum levels), the model will have a charged scalar H + , a pseudoscalar A and two CP-even scalars, h and H. The 2 × 2 CP-even mass matrix is diagonalized via an angle α.
The up and down quark mass matrices are then given by the eigenvalues of which will be the physical quark masses. In fact, these mass matrices will be bidiagonalized in the usual form as where m x are the physical quark masses, whereas V and U are U(3) matrices. These matrices relate the physical quark states u and d to the p and n original states in the following manner: The CKM matrix is then obtained as 1 We are assuming these VEVs are real which is the case of the BGL model but the generalization to complex ones is trivial.
We also define the following matrices, which end up being related to the Yukawa couplings between the physical scalars and quarks. In fact, with the usual conventions (see for instance [9,35]), the Yukawa Lagrangian for the physical fields may be written as where we used the notation s x ≡ sin x, c x ≡ cos x. On a side note, we can see from the above Lagrangian how in the alignment limit the lighter Higgs' Yukawa interactions are exactly like those of the SM particles: in that limit one has sin(β − α) = 1 -which forces the vertices between h and the electroweak gauge bosons to be identical to those of the SM Higgs particle -and therefore the Yukawa couplings of h to quark pairs are proportional to the quark mass, since the contribution of the N matrices vanishes.
In models with flavour conservation, each family of fermions of the same electric charge couples to a single Higgs doublet, via the imposition of discrete Z 2 [31,36] or global U(1) [32] symmetries. Then, the diagonalization of the M u and M d matrices, Eq. (3), is the same as that of matrices N u and N d and there are no flavour-violating Yukawa interactions mediated by neutral scalars. In general, though, that will not be the case and FCNCs occur at tree level.
The BGL model is based on a symmetry imposed on the whole of the Lagrangian, where some of the quark and scalar fields transform as with θ = nπ, with n an arbitrary integer. All other fields remain invariant under this symmetry. As we see, the symmetry treats differently one of the generations of quarks 2 , since only the "first family" of quarks is affected by the transformations above. In fact, there are six (not counting the leptonic sector) models of the BGL type, which depend on which generation of quarks is chosen in Eq. (8) above. For the scalar sector, the above symmetry transformation yields a Peccei-Quinn [32] scalar potential, which must be complemented with a soft breaking parameter as to yield a massive pseudoscalar particle.
In the quark sector, the symmetry transformations of Eq. (8) impose zero values in several entries of the Yukawa matrices of Eq. (1). They are found to be 2 In fact, these are unrotated quark fields, not yet corresponding to physical quarks, but the principle holds.
where, in general, "×" indicates a complex non-zero entry. Then, the form of ∆ 1 and ∆ 2 implies that the matrix M p is block diagonal, namely and it may be bi-diagonalized by unitary matrices V L and V R of the form with some phase θ R . The shape of V L is crucial for the FCNC suppression. First, though, we see that the matrices V L and V R can simultaneously bi-diagonalize ∆ 1 and ∆ 2 , and therefore M u and N u are both diagonal in the basis of physical up quarks. A simple calculation yields with m u1,2,3 the masses of the up-type quarks. However, we have not yet specified which generation is affected by the BGL transformations. Nevertheless, the above shows that for any choice of quark generation in Eq. (8) the N u matrix is diagonal in the up-type quark mass basis, and therefore -since this matrix contains the quark Yukawa couplings of the neutral scalarsno FCNC occurs in the up sector.
The down-quark sector is a different story: given the form of the V L matrix in Eq. (11) and with the CKM matrix V given by (5), we immediately obtain and the impact of this structure on the left rotation matrix for the down-quark sector is considerable. Indeed, a straightforward calculation yields, for the N d matrix (which, we remind the reader, contains the Yukawa couplings between the physical down-type quarks and scalars) for the diagonal terms, whereas the off-diagonal ones are given by Thus we see that the off-diagonal Yukawa couplings between scalars and down-type quarks -which determine the strength of the FCNC interactions -are CKM-suppressed. There is a freedom to choose the "first" family as any one of the physical quark generations, and therefore one has three BGL models with FCNC in the down-quark sector and without them in the up-sector. An analogous symmetry to that of Eq. (8) associated with a given family of up-type quarks would yield other three models, where CKM-suppressed FCNC occur in the up-quark sector and where the down-sector is free from such flavour violation interactions.
This then is how the hallmark of the BGL models is achieved: a flavour-breaking symmetry, which yields offdiagonal FCNC couplings naturally suppressed by the entries of the CKM matrix elements. In what follows we build a similar model but for the case of three Higgs doublets.

A BGL-like 3HDM
Beyond the aesthetic reason of considering three Higgs doublets in analogy with three fermion families, or the intellectual challenge of attempting to reproduce the BGL structure with a larger scalar sector (see [37] for an earlier attempt), there are other reasons to explore a 3HDM with similarly suppressed FCNCs. The BGL model is quite successful, but recent studies [19] have found that its parameter space may be quite constrained. A possible criticism one may levy at the analysis of [19] is that the latter has extended the BGL structure to the leptonic sector as well -something that is not mandatory as the model has enough freedom to accommodate a flavour-preserving leptonic sector in what concerns the Yukawa interactions -which is what we will consider here. Nonetheless, this shows that even with natural FCNC coupling suppression via off-diagonal CKM matrix elements, the BGL structure can be quite constrained from experimental data. Working within the framework of a 3HDM will in principle imply greater freedom in terms of parameters that can be adjusted to comply with experimental bounds.
There is also another reason, more theoretical and fundamental, to attempt a 3HDM study of the BGL paradigm. In many instances, comparisons of the 2HDM with 3HDMs have revealed how special a model the 2HDM is. To give only a few examples, tree-level vacuum stability against charge breaking or spontaneous CP breaking was found for charge-and-CP conserving minima within the 2HDM [38][39][40], but charge breaking minima were found to coexist with charge-preserving ones for NHDMs with N ≥ 3 [41]; a full listing of all possible symmetries of the SU(2) × U(1) invariant 2HDM was found [42,43] while for 3HDM we refer to [24,44,45]; generic boundedfrom-below [42,43] and unitarity [46] bounds were found for the 2HDM, but for the 3HDM such bounds only exist for particular versions of the model. As such, the possibility of ascertaining whether the BGL structure can be extended to a full 3HDM compels us to try to find it. And of course one can obtain an exact BGL-like 3HDM -all one needs to do is copying the procedure detailed in the previous section for the first two doublets, while keeping the third doublet VEVless. For that to happen one would consider the construction of the Inert Doublet Model (IDM) [47][48][49][50] and impose a Z 2 symmetry on the third doublet, φ 3 → −φ 3 . This model would have FCNCs in the visible sector and also a dark matter candidate stemming from the third doublet. An interesting model would be recovered, however it would not solve our aim to obtain a larger freedom to fit the flavour observables in comparison to what one has in the 2HDM BGL model. Thus we are compelled to consider a 3HDM wherein all doublets acquire VEVs.
Let us start by describing the scalar sector of the model, which contains three spin-0 SU(2) doublets, φ 1 , φ 2 and φ 3 .

The scalar sector
The scalar doublets are made to transform under the U(1) × Z 2 symmetry as follows: We further require the scalar potential to be CP-invariant, i.e to be invariant under the usual CP transformations, such that it can be written as The most general potential that softly breaks the U(1) × Z 2 flavor symmetry is given by Here, due to the CP symmetry all the parameters in V = V (φ 1 , φ 2 , φ 3 ) are real. We have introduced real soft breaking terms, in V soft , to avoid the appearance of a massless axion in the physical spectrum. Recall that the same procedure was necessary for the 2HDM BGL [51,52], due to the analogous breaking of the U(1) symmetry given by Eq. (8).
After spontaneous symmetry breaking, all doublets acquire real VEVs 3 and are expanded as where v k represent the VEVs of each doublet which satisfy v 2 The minimization of the potential yields three equations that can be conveniently resolved by expressing the quadratic mass parameters µ 2 1 , µ 2 2 and µ 2 3 in terms of the three VEVs and other couplings as follows: For latter use, we parameterize the VEVs as, and setting v 13 = v 2 1 + v 2 3 , define the following orthogonal matrix which rotates the gauge eigenstates into the so-called Higgs basis, greatly simplifying the analysis of the scalar sector, We now turn our attention to the physical scalar spectrum of the model. Since we are considering a potential with explicit CP conservation and a vacuum which does not spontaneously break CP, the neutral scalars have definite CP quantum numbers. The scalar spectrum of the model is composed of a pair of pseudoscalars, a trio of CP-even scalars and a pair of charged scalars, to be discussed in what follows.
In this work, we have studied the properties of the Higgs sector in the so-called Higgs alignment limit such that one of the physical scalars coincides with the SM Higgs boson (i.e. features its mass and interactions). In order to ensure this in the input data prepared for our parameter scans we would like to utilise an inversion procedure and require the alignment limit at the level of input parameters. Such an inversion procedure would enable us to express the parameters of the scalar potential in terms of physical masses, VEVs and mixing angles.
The mass terms for the pseudoscalar sector can be straightforwardly extracted from the scalar potential -they will correspond to the terms quadratic in the z k (k = 1, 2, 3) fields, after one has replaced the expression for the doublets of Eq. (20) into the potential (18) and (19). One obtains where M 2 P is the 3 × 3 pseudoscalar mass matrix which takes a block-diagonal form in the Higgs basis with the help of the orthogonal matrix O β defined in (23), where Thus, apart from the three VEVs, only λ 10 and the soft mass parameters, µ 2 12 , µ 2 13 , µ 2 23 , enter the pseudoscalar mass matrix. The line and column of zeroes in this matrix obviously tells us that it has a zero eigenvaluewhich of course is the neutral Goldstone boson responsible for the longitudinal polarisation of the massive Z boson.
The pseudoscalar mass matrix is further diagonalized to the mass basis via an additional rotation, where the matrix O γ1 is defined as The full diagonalization of M 2 P from Eq. (24) is therefore accomplished with the matrix product O β · O γ1 , and this is important when one wishes to write the Yukawa interactions between the two physical pseudoscalars and the physical quarks.
An analogous procedure can be performed in other sectors. For instance, the 3 × 3 charged Higgs sector mass matrix, M 2 C , can also be block diagonalized by the same orthogonal matrix as: where again the single line and column of zeros yields a zero eigenvalue, the massless charged Goldstone boson providing the longitudinal polarisation of the massive W boson. In the above matrix, we find Then, one switches to the mass basis in the charged scalar mass matrix as follows, with the charged mixing matrix and where m C1 and m C2 denote the masses of the two physical charged scalars, H ± 1 and H ± 2 , respectively.
Repeating the procedure above also for the CP-even states, we obtain where M 2 S is a 3 × 3 symmetric mass matrix. In explicit form, Switching to the Higgs basis,  we notice that the state H 0 has the same gauge and Yukawa couplings at tree level as those of the SM Higgs boson.
The physical CP even scalars, h, H 1 and H 2 , are obtained via a different orthogonal rotation: where O α is a 3 × 3 orthogonal matrix which can be conveniently parameterized as with, Therefore, M 2 S should be diagonalized via the following orthogonal transformation: where the first eigenvalue corresponds to the SM-like Higgs boson state h, with m h 125 GeV.
Last but not least, it is instructive to discuss the alignment condition in our model. While a rotation to the Higgs basis is performed with the O β matrix such that Eq. (34) holds, the mass eigenstates can be written in terms of the Higgs basis ones as  where one defines the matrix The alignment limit is achieved once the SM-like Higgs boson completely overlaps with H 0 which, in practice, results in the condition With our VEV parametrization in Eq. (22) the alignment condition can be cast as cos α 1 cos β 2 sin(α 2 + β 1 ) + sin α 1 sin β 2 = 1 (41) with reduces to the result in [53] if one identifies This means that putting m h = 125 GeV, α 1 = β 2 and α 2 = −β 1 + π/2 will ensure the presence of a 125 GeV SM Higgs boson in the spectrum -that is the exact alignment limit of this model, forcing the interactions between h and the electroweak gauge bosons Z, W (as well as with SM fermions, see below) to be exactly identical to those of the SM.
In practice, the exact alignment implies that (M 2 S ) 11 = m 2 h and (M 2 S ) 12 = (M 2 S ) 13 = 0 where we define the Higgs basis mass matrixM This can be further solved with respect to λ 1 , λ 2 and λ 10 such that one can write At this point, it is instructive to summarise the above steps. First of all, in order to make our numerical calculations technically feasible and time efficient, in this work the analysis of the scalar spectrum (couplings, mixing and masses) is performed entirely at tree level. We note that the scalar potential in Eqs. (18) and (19) contains sixteen real parameters. Among them, the quadratic parameters µ 2 1 , µ 2 2 and µ 2 3 can be traded in favor of the three VEVs, v 1 , v 2 and v 3 or equivalently v, tan β 1 and tan β 2 . In our numerical studies we take advantage of the exact alignment limit in order to randomly sample tan β 1 , tan β 2 , λ 3,...,9 as well the soft parameters µ 2 13 , µ 2 21 and µ 2 23 such that, using Eq. (44), one obtains the correct λ 1 , λ 2 and λ 10 quartic couplings compatible with an exact alignment of the SM-like Higgs boson. While off-alignment deviations are beyond the scope of this article, we provide in Appendix A generic formulas to obtain the gauge eigenbasis parameters if the physical masses and mixing angles are provided as inputs.

The Yukawa sector
Alongside the scalar field transformations of Eq. (16) the following quark fields are assumed to transform nontrivially under the U(1) × Z 2 flavor symmetry: with α the same arbitrary phase of Eq. (16), and the rest of the quark fields remain unaffected under said symmetry transformations. In Eq. (45), as before, Q La = (p La , n La ) T denotes the left-handed quark doublet of the a-th generation whereas p Ra and n Ra denote the a-th generation (unrotated) up (positively charged) and down (negatively charged) type quark singlets respectively. Notice the similarity between these transformation laws and those of the BGL model, Eq. (8).
The quark Yukawa Lagrangian for a 3HDM will then have the general form where as before Γ k and ∆ k stand for the Yukawa matrices in the down and up quark sectors respectively. Due to specific charge assignments given for the Higgs doublets and quark fields under U(1) × Z 2 these Yukawa matrices will have the following textures: Therefore, the quark mass matrices that emerge from these Yukawa matrices have the following structure: We then rotate from the p and n fields to the physical quark states u and d via rotation matrices V L , V R , U L and U R identical to those of Sec. 2. We thus obtain diagonal mass matrices as in Eq. (3), and the CKM matrix is, as before, given by V = V † L U L . Let us now analyse carefully the Yukawa couplings between the neutral scalar eigenstates and the physical quarks, with particular attention to any FCNC couplings which may arise.
In the alignment limit, with α 1 = β 1 and α 2 = β 2 , the physical scalar h completely overlaps with H 0 . In that limit, the other physical scalars, H 1 and H 2 , will, in general, be an orthogonal mixture of the intermediate states defined above, H 1 and H 2 . Now, the terms in the Yukawa Lagrangian pertaining to the interactions between CP-even scalars and quarks are from which, using the rotation matrix of Eq. (34) to express the h k in terms of H 0 , we can obtain In writing the last step, we have made use of Eqs. (3) and (48). Thus we see that H 0 possesses SM like Yukawa coupling at tree level. This is a close analogy to the BGL model, where we explained how, in the exact 2HDM alignment limit, the h state had identical Yukawa interactions to those of the SM Higgs boson.
In a similar manner, we can write down the Yukawa couplings of H 1 and H 2 with the down type quarks as follows: where the matrices N d1 and N d2 are given by To simplify further the expressions for N d1 and N d2 , we go back to the textures of the mass matrices in Eq. (48). From the block diagonal structure of M u , one can conclude that the corresponding bidiagonalizing matrices, V L and V R , should have block diagonal structures too. In fact, we can choose with the understanding that the phase of (M u ) 33 can always be absorbed into (V R ) 33 . Here, unlike the BGL example of section 2, we are choosing to single out the third family. Then, from Eq. (5) we obtain which means that the third row of U L is identical to that of the CKM matrix, as occurred in the 2HDM BGL case. To proceed further, it is useful to define the following projection matrix Thus, in view of the structures of the Yukawa matrices, we obtain the following relations in the down quark sector: Using Eqs. (54) and (56), the expressions for N d1 and N d2 can now be simplified so that: These equations tell us that the FCNC interactions of H 2 are exactly BGL-like -all off-diagonal elements in N d2 are CKM-suppressed. That however is not the case for H 1 -the first term in the right-hand side of Eq. (57a) is a matrix whose off-diagonal entries are suppressed by two CKM matrix elements, as in the BGL model, but the second term's FCNC couplings are suppressed by only one CKM matrix element. To estimate the size of (U R ) 3B , we note that if (Γ 1 ) 31 and (Γ 1 ) 32 were zero in Eq. (47), then we could choose U R to be block diagonal as well. However, in view of the smallness of the off-diagonal elements in the CKM matrix, the elements of Γ 1 should also be quite small. Therefore it is reasonable to assume that the elements (U R ) 3B (for B = 3) are also small. Considering all these facts, one can expect that the FCNC couplings in the down quark sector controlled by N d1 and N d2 will be under control.
In this work, as was stated above, we work in the exact alignment limit implying that the physical state h couples to the SM fermions and gauge bosons with the same strength as the SM Higgs boson. Ultimately, of course, we will be close to the alignment limit -the current LHC data requires it, see Sec. 4.2 -but not necessarily exactly in it, so the physical CP-even/odd scalar states would have FCNC interactions given by linear combinations of the N d1 and N d2 matrices via the rotation matrices defined in Sec. 3.1. As such, one obtains a model that is not as "clean" as the 2HDM BGL, but where one still sees how FCNC interactions arise which are CKM-suppressed due to the symmetries imposed upon the potential -the suppression is therefore natural and not the result of a fine tuning. A quantitative study of the possible impact arising from off-alignment corrections on FCNC observables in the current model is a subject of a dedicated work in the future.
A similar exercise in the up quark sector would reveal that there are no scalar mediated FCNC at tree level in the up sector. This is due to the special structures of the up type Yukawa matrices, which, in turn, are dictated by the fermionic charges given in Eq. (45). But it should be noted that, just like in the usual BGL models, the charges in Eq. (45) can be appropriately shuffled so that the tree level FCNC couplings reside entirely in the up quark sector instead of the down quark sector. And within each sector one still has the possibility of choosing FCNC associated with a given family. However, we choose the current variant -FCNC in the down sector, associated with the third family -because it will be the most restrictive one. There is a wealth of experimental data limiting such FCNCs, thus this version of the model will be the most restricted one. Therefore, if this version can survive all constraints we will throw at it, other versions would more easily pass the same constraints and may be the subject of future studies.
In passing, it should be mentioned that the leptonic fields are assumed to couple only to φ 1 in the Yukawa sector. This can be achieved very simply by assigning the following transformations to the leptonic fields where L La = (ν La , La ) T and Ra denote, respectively, the left-handed lepton doublet and the right handed charged lepton singlet of a-th generation. Since the charged leptons receive their masses from a single scalar doublet, there will be no FCNC couplings at tree level in the leptonic sector. Thus, all constraints from observables such as µ → eγ are automatically satisfied, since for simplicity in our model the lepton number is assumed to be exactly conserved. Additionally, it is also worth mentioning that in this first study of a 3HDM BGL-like scenario we do not introduce right handed neutrinos, i.e., neutrinos are assumed to be massless in our model.
One major challenge in producing a BSM theory with a non-trivial Yukawa sector, i.e. with FCNC interactions, resides in being able to successfully fit the quark mass spectrum and CKM mixing, simultaneously. In fact, it is a highly non-trivial -and time-consuming -task to find values for Yukawa couplings and scalar VEVs capable of fitting quark masses which differ by many orders of magnitude. Add to that the difficulty in having to simultaneously being capable of fitting the CKM matrix entries and a simultaneous fit to the quark and scalar sector becomes a very difficult achievement.
In this work, we follow an inversion procedure in some ways similar to that implemented in Appendix A for the scalar sector. In essence, the inversion here means expressing the Lagrangian parameters of a given BSM scenario (partially or completely) in terms of physical observables and mixing angles connecting the gauge basis with the physical one. In the model under consideration, such inversion can be unambiguously realized since the corresponding system of equations is linear and non-singular and thus yields a unique non-trivial solution, upon an appropriate choice of input parameters.
The inversion procedure in the Yukawa sector consists in, literally, inverting the usual fitting logic of the Yukawa sector: instead of scanning over Yukawa couplings and VEVs defining some sort of χ 2 function whose minimization would yield an acceptable quark mass spectrum and CKM matrix, we do the opposite. Namely, quark masses and CKM matrix elements are our initial inputs, and we scan over the bidiagonalization matrices which pass from the interaction basis to the mass eigenstate basis. To and the CKM matrix V defined above. Using the unitarity of the rotation matrices we can invert the above equations, and since the definition of the CKM matrix implies that U L = V L V , we can write where U L has been replaced by V L and the CKM matrix. In Eq. (60) the quark masses, CKM matrix and VEVs will be the inputs (the VEVs obtained from a previous partial scan of the scalar sector, already ensuring that the alignment limit is satisfied). The unknowns are the rotation matrices V L , V R and U R . Since any U (3) matrix can be parameterized by three angles and six phases we have a total of 27 free parameters with which one would attempt a fit to Eq. (60). However, using our inversion method there is no fit involved and, due to the Yukawa textures, V L,R become 2 × 2 block diagonal matrices. This means that, in the up sector, we can parameterize our rotation to the mass basis with two angles, α L and α R , defined as Therefore, the five real non-zero elements of ∆ in M p are traded for α L , α R and the physical up-type quark masses m u , m c and m t . In the down sector, on the other hand, both U L and U R are generic 3 × 3 matrices. While the former is fixed as U L = V L V , with the CKM matrix expressed in the Wolfenstein form ( see PDG for details [54]), the latter can be parameterized with three angles β 1,2,3 R if we assume, as we do in our numerical implementation, that the only complex CP-phase comes from the CKM matrix. This means that, for our scenario, the seven non-zero elements of Γ in M n can be consistently described with eight real couplings, which are replaced by the four Wolfenstein parameters λ, A,ρ andη, three quark masses, m d , m s and m b and a randomly generated angle which we call β R . In particular, denoting the solutions of Eq. (60) for the angles as , which are numerically computed, one can express U R as with β R a free parameter.
All in all, given the physical quark masses and CKM mixing (with a complex CP-phase), our numerical sampling of the Yukawa sector is consistently achieved solely with three angles, α L , α R and β R , upon inversion of Eq. (60).
Therefore with the VEVs and the values obtained for the nonzero entries, it is a simple matter to reconstruct the elements of the Γ and ∆ matrices -notice that each non-zero entry of M p (M n ) has the contribution of a single ∆ (Γ) matrix element such that the reconstruction of the Yukawa matrices is unequivocal.
With this simple procedure we ensure that the very different quark masses are automatically and exactly reproduced for every single sampled point, and so does the non-trivial structure of the CKM matrix. The computational challenge then is to scan over three rotation angles in order to reproduce the correct flavour observables in consistency with strict experimental constraints.

Constraints on the model
Any realistic BSM theory needs to do, at least, as good a job as the SM in describing well measured particle physics properties. For multi-scalar scenarios, there is a wealth of theoretical and experimental bounds that are imposed upon the model's parameter space. The main challenge then is to find the domains of validity for a given model for which the parameter space points pass all the relevant restrictions. In this section, we summarize the main theoretical and experimental constraints that need to be satisfied in order to validate our BGL-like 3HDM framework.

Theoretical constraints
For models with a scalar content larger than that of the SM, special attention needs to be focused on the possibility of the scalar potential becoming unbounded from below i.e. tending to minus infinity for some direction along which the fields are assuming arbitrarily large values. This imposes constraints on the model's scalar quartic couplings, as the quartic part of the potential clearly dominates over the quadratic (or even an eventual cubic) one when the scalar fields tend to infinity. This is already a concern for the SM -it is the reason why the SM Higgs quartic coupling λ is taken positive. For the 2HDM, generic conditions were found in Refs. [42,43] but for the 3HDM, on the other hand, there is no such generic boundedness-from-below prescription that can be straightforwardly implemented in a parameter scan (see [55] for a U(1) × U(1) version of the 3HDM). Still, some necessary conditions are easy to find. Analysing the shape of the scalar potential of Eq. (18), we see that if one takes each of the doublets φ i to infinity separately, the potential will tend to −∞ unless Likewise, following a procedure similar to the one used in the 2HDM [18], if one takes two doublets (i, j) to infinity but such that φ † i φ j = 0 (which is realised, for instance, when the upper component of one of the doublets and the lower component of the other doublet are zero) one obtains a positive value of the potential for any value of the fields if We can also adapt the boundedness-from-below necessary conditions from Ref. [56] (see Eqs. (21)-(24) there), being careful with the fact that the potential of that work is different from the one considered here (our potential has a more restrictive symmetry, and therefore has fewer quartic couplings). This translates into a generalisation of the above conditions, which become Ultimately, these necessary conditions eliminate a great deal of parameter space. Even though they are not the sufficient ones, they are still expected to cover most of the parameter space regions that lead to a boundedfrom-below potential.
Another strong constraint on the quartic couplings of potential is the requirement for the theory to be unitary. For the SM, this implies an upper bound on the mass of the Higgs boson [57,58], and similar studies have been applied to the 2HDM (general unitarity conditions are presented in Ref. [46]) and other models with a larger scalar content. Essentially, the method consists in computing all scalar-scalar J = 0 scattering amplitudes (usually denoted a 0 ) and requiring that they respect unitarity in the high energy regime. This translates into an upper bound on those amplitudes, |a 0 | < 1/2. Theories with many scalars complicate the calculation somewhat due to a growing multiplicity of such scattering amplitudes and bounds must be imposed upon eigenvalues of the S-matrix. General unitarity bounds for a 3HDM with a Z 2 × Z 2 -symmetric potential are presented in Ref. [56], of which our U(1) × Z 2 symmetry is a special case. Since our model has a larger symmetry, it has less parameters than that of Ref. [56] and one could read off from Eqs. (91) to (102) of that article the unitarity constraints imposed upon the quartic couplings. However, as explained below, we will use the SARAH/SPheno machinery instead to take into account these bounds.
Finally, a "standard" constraint on multiscalar models is to verify their compliance with electroweak precision bounds i.e. constraints on the oblique S, T and U parameters. Models with N Higgs doublets automatically satisfy ρ = 1 at tree-level, which means that bounds on the oblique parameter S will be easily satisfied. However, that is no longer true for the T parameter, which typically needs to be computed for any given model. Constraints on T typically enforce, for very high masses, degeneracies between the extra scalars in the mass spectrum of the model. The results from Ref. [56] are of no help to us in this case, as the expressions for the oblique parameters given there are only valid for a 3HDM version of the inert type (where one of the doublets is VEVless and naturally decouples from the gauge sector).
Instead of computing the unitarity bounds and oblique corrections analytically, in this work, we have adopted another strategy and used the publicly available SARAH/SPheno framework [59][60][61] which enables one to compute them numerically in a given particle physics model for a particular parameter space point (for an earlier implementation of this procedure, see e.g. Ref. [62]). In our numerical analysis substantiated in detail below, we have implemented the 3HDM model under consideration into this framework where both the unitarity constraints and the electroweak precision bounds on S, T, U parameters have been consistently incorporated in the parameter scans designed to search for the valid physical regions of the model. In particular, the oblique parameters are computed by SPheno-4.0.4 at one loop-order using tree-level masses and then are verified against the experimental bounds for every parameter space point.

Experimental constraints
Since the discovery of the Higgs-like state in 2012 the LHC collaborations have been able to verify that its properties correspond to those expected for the SM Higgs boson, within rather small experimental error bars. In practical terms, this means that the couplings of the 125 GeV h state to electroweak gauge bosons and fermions in a BSM model cannot deviate too much from the corresponding SM couplings. A convenient way of parameterizing the Higgs interaction strengths is by introducing the κ-formalism, defining the dimensionless quantities, through the decay widths, Γ, of the Higgs boson to a certain final state X (typically, ZZ, W W , ττ and bb), computed both in a given BSM scenario as well as in the SM. This definition implies that the exact SM behaviour would correspond to κ = 1. Notice that current LHC measurements [3,4] only allow deviations from unity roughly up to 20% for the several final states.
Requiring h to be SM-like naturally suppresses the couplings to gauge bosons of the heavier CP-even states H 1 and H 2 (these three states' gauge couplings obey a sum rule, due to gauge invariance). Since one of the most constraining channels in the search for heavier resonances at the LHC is precisely via di-Z boson production, most of the constraints coming from those searches are automatically satisfied. However, there is still an ample parameter space for which the heavier scalars also have suppressed production cross sections and could thus have avoided detection so far.
The basic properties of the Higgs boson in the BGL-like 3HDM under consideration, such as its physical couplings, decay rates and branching ratios, have been computed numerically in SPheno-4.0.4 for each input parameter-space point in a dedicated numerical scan. A set of input points has been prepared in a separate special routine using an inversion procedure in the scalar sector in order to reproduce the exact Higgs alignment limit, with fixed VEV v 246 GeV and Higgs mass m h 125 GeV values, as detailed above. Then, a candidate point has been chosen to satisfy the electroweak precision bounds, the boundedness-from-below and unitarity constraints discussed in the previous subsection. On top of that, in order to take into account the latest data from the LHC on the properties of the observed 125 GeV state we used HiggsSignals-2.2.3beta [63], whereas limits from the LEP, Tevatron and LHC direct searches for new CP-even and CP-odd scalars were taken into account through the use of HiggsBounds-5.3.2beta [64][65][66]. The latter two packages are linked to SPheno-4.0.4 such that they get all the necessary inputs containing masses and decay widths of the SM and BSM Higgs states for each parameter-space point. A point that passes all these constraints is considered to be valid, at least, from the point of view of the scalar sector constraints.
As for the flavour sector, there are numerous observables that have to be computed and checked for each parameter space point that satisfies the scalar sector constraints. The inverted procedure outlined in Sec. 3.2 already ensures that our choice of the parameter-space points, even before a numerical scan, automatically reproduces the correct measured values of the quark masses and of the CKM matrix elements. But the presence of charged scalars, as well as neutral ones with tree-level CKM-suppressed FCNC interactions, means that we must verify particularly sensitive quantities, such as the b → sγ width, the neutral K-and B-meson mass differences and the CP-violating K phase, among others. To this end, we employ the FlavourKit package [67] as part of the SPheno tool in order to compute the full list of Wilson coefficients for a given parameter-space point. The Wilson coefficients are then passed over to the Flavio python package [68] enabling to compute all the relevant flavour physics observables for a given point and to compare them to the corresponding SM values.
In particular, out of an extensive list, we have found that the most relevant 4 quark flavour violation (QFV) observables to consider are B → χ s γ, B s → µµ, ∆M s , ∆M d and ε K , which we collectively denote as O 3HDM . In our analysis we quantify deviations from the SM prediction as a ratio O 3HDM /O SM . Defining the QFV experimentally measured value, its experimental error and the SM prediction theoretical uncertainty as O Exp , σ Exp and σ SM respectively, one can use the error propagation formalism to obtain the error associated to the ratio O Exp /O SM that reads as In Tab. 1 we numerically determine σ for each of the five QFV observables mentioned above which will later be used to define the error bars in our plots.  Given the BGL-like nature of the FCNC interactions in the considered 3HDM, the QFV observables appear to be already in the right ballpark of typical values being not too far from the SM results (see below). Thus, there is no need to include the flavour observables in the numerical fit -we have just computed those observables for each point that has passed the theoretical, electroweak precision tests and Higgs sector constraints. In other words, when scanning over the BGL-3HDM parameter space requiring exact alignment and using the inversion procedure described for the fermion sector, the parameters found already give QFV observable values not too far from the expected values. As such there is no practical need to include fitting those variables along with the scalar sector (a task which would be highly work-intensive), since the overall efficiency of the scan is not compromised by points not complying with QFV observables. Let us now turn to a discussion of our numerical results.

Results
As it was discussed above, we use an inversion procedure in order to automatically obtain the correct quark masses, CKM matrix elements, the measured Higgs boson mass as well as an exact alignment of the latter to the SM, for all sampled parameter space points. For such a purpose we have built an internal routine, that we shall denote as pre-SPheno, where Eq. (44) and an expanded form of Eq. (60) are implemented. With this we compute the gauge eigenbasis Lagrangian parameters and then proceed with the calculation passing them to SPheno. The input free parameters were then randomly sampled in the ranges given in Tab. 2.  [54].
In essence, for for all that follows, we have imposed a priori certain basic constraints implemented analytically: the correct electroweak symmetry breaking must occur, thus the doublets have VEVs satisfying v 2 1 + v 2 2 + v 2 3 = (246 GeV) 2 ; the lightest CP-even mass is fixed to be 125 GeV; the λ 1 , λ 2 and λ 10 quartic scalar couplings where chosen to provide the exact Higgs alignment limit; and the remaining angles in the fermion sector besides α L , α R and β R were calculated to comply with the correct quark masses and CKM mixing. All the other constraints, such as the necessary boundedness-from-below conditions described in Sec. 4.1, the unitarity, electroweak precision and Higgs-sector phenomenological constraints are checked numerically using the coupled chain of computer packages -SPheno, HiggsSignals (HS) and HiggsBounds (HB) -in a dedicated numerical scan, while the flavour physics constraints are verified a posteriori using FlavourKit and Flavio tools.
Let us now investigate the effect that the various constraints, both on the scalar and flavour sectors, have on the allowed parameter space of the model.

Allowed parameter space
We show in Fig. 1 the EW precision observables in the S-U , S-T and U -T projections. Current precision measurements [54] provide the allowed regions, where S-T are 92% correlated, while S-U and T -U are −80% and −93% anti-correlated, respectively. We compare our results with the EW fit in Eq. (68) and require consistency with the best fit point at 95% C.L. taking into account the correlation between the oblique parameters.
We see that even with our generation of parameters satisfying exact alignment compliance with electroweak precision constraints is not guaranteed, and plenty of points are rejected. However, as expected in a model with multiple doublets, it is not difficult to find regions of parameter space for which all constraints on the oblique parameters are satisfied. Points outside 95% CL for the oblique parameter STU analysis Points within the 95% CL for the oblique parameter STU analysis In Fig. 2 we show the effect of non-flavour constraints on the allowed parameter space. Here, the impact of restrictions from the LHC experiments, both in terms of measurement of the Higgs bosons' properties or in the searches for extra scalars, incorporated in the HS/HB framework have been analysed. Unitarity bounds on the model's quartic couplings are also imposed, as well as precision electroweak constraints via the S, T and U parameters, each leading to a considerable reduction of the allowed parameter space. We see a close correlation between m A1 and m H1 for large values, stemming mostly from unitarity and electroweak precision constraints. Note, the same tendency of near-degeneracy is observed in the mass spectrum of the 2HDM. Furthermore our scan generates very low masses for the scalars, which are excluded by various direct collider searches and implemented in HS/HB. It is important to mention that the size of the input soft masses, together with that of the quartic couplings in Tab. 2 set, approximately, the scale of the physical BSM scalars.
In Fig. 3 we show how some of the QFV observables might further constrain the model's parameter space that survives the Higgs physics, unitarity and electroweak precision constraints. For instance, the dependency of the ratio of the b → sγ width computed in the BGL-like 3HDM to the expected SM value as a function of m H1 is shown in Fig. 3(a). Here, we observe a dispersion around the SM value such that some of the points deviate by more than 2σ. The 1σ band is defined in the first line of Tab. 1. In analogy to many known versions of the 2HDM, the b → sγ constraint is a very important one, excluding a number of parameter points which otherwise could be perfectly acceptable. Not all flavour variables yield strong constraints, though -in Fig. 3(b) we show the values obtained within our parameter scan for the Kaon system CP-violating K phase. One notices a rather minuscule variation around the SM value after all other QFV observables have been constrained to lie within a 2σ interval of their respective SM-expected values. This is clearly an indication that there are no substantial FCNC contributions to this observable in the considered BGL-like 3HDM.
Note that in Figs. 2 and 3 we have shown only the regions of the parameter space where all constraints from boundedness from below, unitarity, electroweak precision variables and direct searches are obeyed. Having ascertained the relevance of the constraints imposed on the scalar sector, we wish to analyse in detail the impact of the QFV observables, leaving temporarily aside the remaining phenomenological constraints. Since our model has tree-level FCNCs, the inverted procedure described in Sec. 3.2 does not immediately guarantee a good fit to QFV observables such as the mass differences of the neutral K, B d and B s mesons, the already mentioned K CP phase, branching ratios such as B s → µ + µ − etc. Since the SM already does a good job at describing the quark and lepton sector behaviour, we need to verify that NP contributions to those quantities do not ruin the existing agreement between theory and experiment.   As we have already mentioned, our model, being BGL-like, ought to allow for an easy fit in the flavour sector, an assumption we now put to the test. In Fig. 4 we see how an agreement with b → sγ constraints is not automatic -in fact, we observe a number of (grey) points with larger than 2σ deviation from the SM prediction, and quite a few disagreeing with the SM b → sγ values between 1 and 2σ. These larger than 1σ deviations are seen to occur mostly (though not exclusively) for lower masses of the lightest charged scalar, m ± H1 ≤ 700 GeV. This is similar to what occurs for the Type II 2HDM [70]. Since in our model the down-type quark masses do not arise from a single Γ matrix, it is in fact natural that we find regions of parameter space for which observables such as b → sγ behave in a similar manner to a Type II 2HDM. However, we see that our 3HDM can fit this observable for much lower charged Higgs masses than 700 GeV, as occurs, for instance in a Type I 2HDMagain to be expected, certain regions of our parameter space should mimic well the Type-I behaviour. A similar phenomenon was observed for a 2HDM with tree-level FCNCs, see [9].
We further observe that the values of the B → X s γ width in our model approach the corresponding SM value for very large values of the lightest charged Higgs boson mass 5 . This is not surprising since NP contributions to this observable depend on the inverse of the square of the extra scalars' masses and are thus expected to approach zero as those masses tend to infinity. In Fig. 5 on the other hand, considering again the full set of phenomenological constraints, we observe how the inverted procedure we are using to constrain the Yukawa sector yields an excellent agreement with other QFV observables -there we plot the values of K as a function of the lightest charged Higgs boson mass, and see how close it gets to the SM value for all the generated points. We see that this observable attains, in this model, values extremely close to the SM prediction, with deviations of the order of ∼ 0.01%. To put these results in context, the current experimental uncertainty on K stands at less than 0.5% of its central value. The minimal value of the charged Higgs boson mass that still reproduces the experimental value of K and satisfies all constraints is found to be ∼ 150 GeV.  In the left panel, (a), we show the results for b → sγ, namely, the ratio of 3HDM-to-SM branching fractions for B → X s γ reaction while in the right panel, (b), we plot an analogous ratio for K , both in terms of the H 1 mass. The colour code is as in Fig. 2 and grey points are excluded, at 2σ level, by at least one QFV observable.
For completeness, let us also consider the B-meson mass differences. These are the observables which in the SM are generated by one-loop box diagrams but also receive tree-level contributions in theories with scalar mediated FCNC interactions in the down-quarks sector. Again, and as expected, we see in Figs. 6 and 7 that the values obtained in the BGL-like 3HDM for ∆M Bs and ∆M B d approach their SM values for large enough masses of the extra scalars. We also see that our scanning procedure produces values of ∆M B d extremely close to that of the SM (even for lower masses), with a larger dispersion found in ∆M Bs , still within a 2σ variation. This is clearly due to the fact that we chose a specific structure for the Yukawa matrices in Eq. (57) in order to single out the third generation. Furthermore, for a BGL-like model, the FCNC interactions are expected to be suppressed by the CKM matrix elements, which, for the B-meson oscillation observables under consideration, explains how contributions to ∆M B d , which involve a "jump" across two generations, are more suppressed than those contributions to ∆M Bs , for which scalars only "jump" one generation in their QFV interactions.
While we do not show all the numerical results explicitly, we have analysed a wealth of other flavour physics observables, encountering 1σ agreement with current experimental bounds for all of them. These included the remaining QFV observables such as neutral Kaon mass differences, neutral B mesons decays to muon and electron pairs and other leptonic sector measurements, Z → bb observables etc.

A degree of fine-tuning
The remarkable agreement found in the previous section for QFV observables, even when including the FCNC contributions from the extra scalars present in our model, needs to be analysed in depth. Such an agreement can arise in several ways. For instance, if all extra scalar boson masses are very high, then the NP contributions take very small values and an agreement with the SM value is easily found. Another possibility is that the FCNC Points consistent with all QFV observables within 1-2σ uncertainty Points ruled out by at least one QFV observable Figure 4: The ratio of 3HDM-to-SM branching fractions for B → X s γ decay in the considered BGL-like 3HDM as a function of the H ± 1 mass. Red-shaded points are those for which all QFV constraints are satisfied at most at 2σ (many even at 1σ) level, grey points are those for which at least one QFV observable is in disagreement with the current measurements by more than 2σ. The horizontal lines account for the 1 and 2σ uncertainties as in the first line of Tab. 1. The "temperature" gradient of colour shows the lightest pseudoscalar mass. Yukawa interactions are naturally small, something which occurs in the 2HDM BGL model. Their strong suppression follows from the small off-diagonal CKM matrix elements as a consequence of the symmetry of the model, which we claim to also happen in the current 3HDM version. And finally, there is also the possibility of a fine tuning having occurred in the scanning procedure, "unnaturally" causing cancellations between different NP contributions. Let us now demonstrate that this last possibility does not occur in our model.   Figure 7: B d mass difference as a function of heavier CP-even Higgs boson masses. The colour code is the same as in Fig. 4.
In order to see an example of such possible fine tuning, consider the Kaon system and observables such as the K 0 −K 0 mass difference, or the CP-violating phase K . The matrix element describing the transitionK 0 → K 0 , M 21 , receives contributions from the SM, via box diagrams, and from NP, through tree-level FCNCs in the scalar sector: The NP terms arise in our model from tree-level Feynman diagrams and thus can in principle overwhelm the SM result. These diagrams represent the tree-level exchanges of CP-even and CP-odd scalars with FCNC interactions which, using the vacuum-insertion approximation (see Refs. [6,7,9]), are found to be where f K and m K are the K-meson decay constant and mass, respectively, and the following combinations of FCNC couplings were defined, for each scalar/pseudoscalar X = {h, H 1 , H 2 , A 1 , A 2 }, as follows Here, the matrices N dX are the Yukawa matrices for down-type quark interactions with each scalar X. We observe in Eq. (69) that CP-even contributions tend to cancel CP-odd ones. This opens up the possibility to fit the Kaon observables if the CP-even and CP-odd terms, though potentially large by themselves, cancel up to the n th decimal place. Therefore, such a cancellation would produce a result that, while in nominal agreement with experiments, is "unnatural", and would be potentially challenged by higher order corrections. Similar fine-tunings should be investigated in other observables, such as the B-meson mass differences, or semileptonic quantities such as the branching ratio of B s → µ + µ − etc.
In order to investigate whether our results were fine-tuned or not, we have undertaken the following procedure: • Chose a combination of parameters for which all theoretical and experimental constraints are satisfied.
Here, we were focused on the parameter space points that give rise to the predicted values for QFV observables that are at most within 2σ of their SM counterparts.
• Fixed all the model parameters except for the soft breaking mass terms, which were then smeared with a random variation of no more than 1% about their initial values. The scalar mass spectrum was then re-calculated and we chose those situations for which there were variations of less than 5% on all scalar boson masses.
• With the new values of the scalar boson masses (and all remaining mixing angles and Yukawa couplings being fixed to the values prior to the smearing procedure described above) we recalculate the QFV observables and compare them with their initial values.
If there is a fine tuning in the calculation of K , for instance, a small variation in the masses, such as m A1 or m H1 , should lead to a much larger variation of the observable. We see the results of this procedure in Fig. 8. Here we show variations observed for the branching ratios of B → X s γ and B s → µ + µ − decays. Similar variations were found for other observables, such as the meson mass differences ∆M d and ∆M s or the branching ratio for B s → µ + µ − decay. Other observables had by far smaller variations -for comparison, the absolute value of variations on K , for instance, were always found to be less than 0.01%.
Thus, we observe that, while keeping all mixing angles and Yukawa couplings the same, smearing the physical scalar masses by up to 4% around their starting values for phenomenologicaly viable points, induces variations of less than 2.5 % on the values of typical QFV observables. We therefore conclude that in our BGL-like 3HDM the best parameter space points that pass all the considered constraints are not the result of an accidental fine-tuning, rather the FCNC interactions of the model have been rendered naturally small by the symmetries of the model, just as it occurs in the standard 2HDM BGL.

LHC predictions
As we saw in the previous sections, the BGL-like 3HDM under discussion can fit, despite the presence of tree-level FCNC interactions in the down sector, all current experimental constraints without any fine tuning. Furthermore, it can do so even with extra scalar masses below ∼ 500 GeV raising a tantalizing question: could such scalars have already been observed at the LHC, or can the current LHC data be used to exclude portions of the model's parameter space? Besides, what can one expect vis-à-vis future LHC sensitivity to potentially discover the extra scalars predicted in our model?
To begin with, we must recall that we have chosen to work in the exact alignment limit in this model. Thus, the CP-even scalars H 1 and H 2 have couplings to Z (and W ) boson pairs exactly equal to zero. As such, experimental searches for extra scalar resonances decaying into Z or W boson pairs [71][72][73] are automatically satisfied in our parameter scan. Deviations from the alignment limit to be considered in a more general analysis might somewhat change that state of affairs. Given how SM-like the 125 GeV scalar appears to be in the LHC measurements, the couplings of H 1 and H 2 to electroweak gauge bosons should always be heavily suppressed, so we do not perform such more generic off-alignment analysis in the current work. However, the sum rule for scalar-gauge boson couplings -which in the exact alignment limit makes H 1 and H 2 gauge-phobic -does not apply to Yukawa interactions. For instance, in the 2HDM or in SUSY models, interactions of the pseudoscalar A to fermion pairs may be enhanced (or suppressed) by a factor of tan β. A promising avenue in searches for additional scalars is therefore the di-tau channel, where the current and future LHC sensitivities may well reveal their presence.
In Fig. 9 we show the cross section for the production of the lightest pseudoscalar A 1 in a gluon-gluon fusion process multiplied by its branching ratio into tau pairs, for a center-of-mass energy of 13 TeV. The CMS 65% and 95% exclusion bounds were extracted from [74]. The pseudoscalar production cross section via gluongluon fusion was obtained using MadGraph5 aMC@NLO 2.6.2 [75]. We have used an interface with SPheno where masses, decay widths and branching fractions are calculated and linked to MadGraph5 for every single generated point. The relevance of requiring an agreement with QFV observables is emphasized in this plot, where grey points represent scenarios where at least one of such observables was found to deviate from its SM value by more than 2σ. The blue points denote scenarios for which a full 1σ agreement was found for all QFV observables while the darker shades of blue also survive all remaining theoretical and experimental constraints. As was to be expected, the total signal strength for this channel diminishes as the pseudoscalar mass increases -see in particular the sharp drop-off around m A1 ∼ 375 GeV, that is twice the top mass. Indeed, for larger pseudoscalar masses, a new decay channel A 1 → tt becomes kinematically allowed and tends to reduce the branching ratio for A 1 → ττ . Notice, however, that there is a number of points for lower m A1 which can almost be probed by the current CMS bounds.
Until the end of LHC operation we can expect an increase in accumulated luminosity by at least a factor of 100, which would roughly lower the exclusion lines shown in Fig. 9 by an order of magnitude. As such, we can expect the searches in this channel to at least exclude parts of the parameter space for m A1 < 400 GeV. In fact, we see in Fig. 9 that the maximum of the signal strength occurs for m A1 350 GeV, which is unsurprising, given that this value roughly corresponds to twice the top mass. In fact, it is well known that the gluon-gluon fusion CMS 95% C.L. Bound CMS 65% C.L. Bound Excluded Points Points consistent with QFV observables Points that survive all constraints Figure 9: The pseudoscalar A 1 production cross section in gluon-gluon fusion at the LHC center of mass energy of 13 TeV, times its decay branching ratio to ττ as a function of m A1 . Grey points represent the sets of scenarios excluded by not obeying at least one QFV observable at 2σ. Blue points, both dark and light shades, correspond to an agreement with all QFV observables at least at 2σ level. Those points that are further allowed by all imposed constraints are represented by the dark blue points. The 1 and 2σ observation limits available from the CMS Collaboration for searches in this channel are taken from [74].
cross section has a local maximum for a c.o.m. energy equal to twice the top mass, both for the production of a CP-even or a CP-odd scalar.
The di-tau channel is also appropriate in searches for a heavier CP-even state, as we see in Fig. 10. As before, take notice of the expected sharp drop in the value of the signal rate for masses m H1 > 2m t . Both in direct production via gluon-gluon fusion into H 1 , or in its associated production with a bottom quark pair, the obtained signal strength including the branching ratio for H 1 → ττ is very close to the current CMS sensitivity for the lower mass region. Thus, we see that our BGL-like 3HDM is close to being probed by the current LHC data, and before the end of the next LHC run certain parts of its parameter space can also be tested in direct searches for BSM scalars. We therefore provide five representative benchmark points in Tab. 3 to be searched for in the LHC run-III. These were chosen such that they obey all theoretical and experimental contraints on the scalar, gauge and fermion sectors, and further satisfying the following criteria: • BP1 corresponds to the lightest CP-even BSM Higgs boson found in our scans with mass m H1 = 249 GeV. This point also corresponds the lightest charged Higgs scenario with m H ± 1 = 101 GeV; • BP2 represents the second-to-lightest CP-even and charged BSM Higgs bosons found in our anlysis with masses m H1 = 285 GeV and m H ± 1 = 146 GeV; • BP3 and BP4 correspond to the lithest and next-to-lightest CP-odd Higgs boson found in our scan with masses m A1 = 161 GeV and 206 GeV respectively; • BP5 offers an early discovery or early exclusion scenario in the gg → A 1 → τ τ channel, m A1 = 338 GeV, where the signal strength was found to be the closest one to the CMS bound.
• BP6 corresponds to an early discovery/exclusion scenario in the gg → H 1 bb → τ τ bb channel, m H1 = 313 GeV, where the signal strength was found to be the closest one to the CMS bound.
• Last but not least, BP7 represents an early discovery/exclusion scenario in the gg → H 1 → τ τ channel, m H1 = 353 GeV, where the signal strength was found to be the closest one to the CMS bound.   Table 3: A selection of seven benchmark points. All masses are given in GeV and cross sections in pb. These correspond to the lightest scalars found that respect all QFV, electroweak, Higgs and theoretical constraints (BP1 -BP4), and to three early discovery/exclusion points that mostly approach the experimental bounds in Fig. 9 (BP5) as well as in Fig. 10-right (BP6) and Fig. 10-left (BP7). [pb] (a) (b) Figure 10: The signal strength for production of a CP-even scalar via the gluon-gluon fusion mechanism (a) times its branching ratio to ττ , and (b) with associated production of a bb-pair times its branching fraction to ττ , as a function of the lightest CP-even mass, m H1 . The colour code is the same as in Fig. 9 and the exclusion bounds in both panels were also taken from [74].
Note that the entire scalar spectrum in BP3, BP4 and BP5 is lighter than 1 TeV and potentially at the reach of the LHC run-III. Furthermore, it is remarkable to note that the lightest charged Higgs in BP1 is allowed to be lighter than the SM Higgs boson while conforming with all experimental constraints. On the other hand, in BP1 and BP2 the heavy scalar masses m H2 , m A2 and m H ± 2 are larger than 4 TeV while in BP6 and BP6 their masses are approximately 1.5 TeV and 1.1 TeV. We also provide in Tab. 3 both the production cross sections and the branching fractions calculated for each of the studied channels as well as the 3HDM-to-SM ratio of each of the five QFV observables in Tab. 1. While the former are relevant for direct searches for new scalars at the LHC, the latter may be probed in flavour experiments.
An extended parameter space domain which relaxes the condition of exact Higgs alignment -while maintaining the full agreement with the current LHC measurements regarding the properties of the lightest Higgs bosonwould no doubt show a larger excluded region by means of the measurements in the gg → H 1 bb → ττ bb channel.
Other search channels might also be considered (such as searches in decays to top pairs) but we relegate those analyses for a future work where one would need to perform a more thorough scan of the model's parameter space. Our main goal here is to prove, with a simple example, that the model is of interest for the ongoing LHC searches.

Conclusions
We have presented a new BGL-like 3HDM model with symmetry-suppressed FCNC interactions mediated by neutral scalars, in a close similarity the well-known example of BGL 2HDM. We have applied the basic theoretical (unitarity, boundedness from below) and experimental (electroweak, Higgs and flavour) constraints and identified the domains of validity of the model and main phenomenological implications. Our analysis has enabled us to narrow down the allowed parameter space regions which simultaneously fit all theoretical and experimental constraints. We have also discussed the possibility of probing the model at the LHC via gluon fusion production of new CP-even and CP-odd Higgs bosons and subsequent decay into ττ pairs as well as via an associated production of the new CP-even Higgs state and bb pair. In particular, we have observed that the BGL-like 3HDM offers a possibility for lighter than conventionally allowed non-standard scalars, at the reach of the LHC III. We have identified and described seven benchmark scenarios that can be used in experimental searches for Higgs partners at forthcoming LHC runs.
Our analysis determined the most sensitive flavour violation channels, and has revealed that the BGL-like mechanism induced by the U(1) × Z 2 flavour symmetry is indeed responsible for the suppression of FCNCs rather than any accidental cancellation or any fine-tuning. Indeed, it results from the BGL nature of the 3HDM under consideration that the most stringent constraints on the model's parameter space are not the QFV ones but rather the Higgs physics observables from direct searches. As one of the possible future avenues, theoretical investigation of the proposed BGL-like mechanism can be continued towards generalisation of the lepton and neutrino sectors by implementing the U(1) × Z 2 symmetry there and studying its consequences on LFV and neutrino mass generation.

A Generic off-alignment conditions
An alternative inversion procedure to the one used in our numerical calculations can be given. In particular, we provide in this appendix closed expressions for the quartic couplings λ 1,...,10 and for the soft breaking mass terms µ 2 12 and µ 2 23 in terms of all physical scalar masses, mixing angles and VEVs.
Equating these to the corresponding quantities in Eq. (26), one obtains three equations which can be resolved with respect to the potential parameters, for example, λ 10 , µ 2 12 and µ 2 23 in terms of two physical masses, m A1 and m A2 , one mixing angle γ 1 , the third soft mass parameter, µ 2 13 , and the Higgs VEVs (or, equivalently, mixing angles β 1,2 and v). The results reads, For the charged scalar masses one can use the expressions for λ 10 , µ 2 12 and µ 2 23 found above and, solving the eigenvalue problem for Eq. (31), extract the relations for other three parameters of the potential -λ 7 , λ 8 , and λ 9 couplings -in terms of the physical masses of the H ± 1,2 and A 1,2 states, the µ 2 13 parameter, the mixing angles γ 1,2 and the Higgs VEVs (or, equivalently, mixing angles β 1,2 and v). The result reads as Finally, for the CP-even scalar masses, inverting Eq. (37) we get, which enables us to solve for the remaining six quartic couplings λ 1,2,3,4,5,6 in terms of the physical CP-even scalar masses, the Higgs VEVs and the mixing angles α 1,2,3 , as well as the quartic couplings λ 7,8,9,10 , and the soft parameters µ 2 12 and µ 2 23 , obtained above for the CP-odd and charged scalar sectors. The final expressions can be expressed as − v 2 c α1 m 2 H1 c α3 (s α1 c α2 c α3 + s α2 s α3 ) + m 2 H2 s α3 (s α1 c α2 s α3 − s α2 c α3 ) In essence, a generic, off-alignment scan, can be performed by exchanging µ 2 12 , µ 2 23 and λ 10 for m A1 , m A2 and the mixing angle γ 1 using Eqs. (72) to (74). The remaining nine quartic couplings, as discussed above, can be expressed in terms of five physical masses (three CP-even scalars, and two charged scalars), the remaining soft-breaking parameter, µ 2 13 , and four mixing angles (three in the CP-even sector and one in the charged scalar sector).