Classification of anomaly-free 2HDMs with a gauged U(1)′ symmetry

Two-Higgs-doublet models (2HDMs) with a flavored U(1)′ gauge symmetry are a popular extension to the Standard Model (SM), yet they currently lack a complete survey. In this paper, we present a full classification of anomaly-free 2HDMs within the SM fermion content, resulting in a total of eleven distinct models. Four of these are relatively well-studied, while the rest either partially, or entirely, lack previous treatment. We study these textures under a variety of experimental bounds, focusing mainly on the previously unexplored models. This work is intended to act as a catalog to models worth considering in greater detail.


I. INTRODUCTION
With the discovery of the Higgs boson [1,2], the Standard Model (SM) has reinforced itself as the most accurate description of High Energy Physics. While the SM has endured all kinds of experimental scrutiny, it still leaves room for sizable New Physics (NP) effects. For the last decades, many extensions have been proposed. One of the most popular SM extensions is the Two-Higgsdoublet model (2HDM) [3,4], where a single additional copy of the Higgs doublet is added to the SM. Such a simple scenario introduces, nevertheless, new types of interactions; a general 2HDM implementation is plagued by tree-level flavor-changing neutral currents (FCNCs) mediated by the new scalars [5]. To circumvent such a problem, the Yukawa textures are constrained with various additional symmetries. In the so-called natural flavor conserving 2HDM implementations, an additional Z 2 is added, leading to the complete absence of tree-level FC-NCs [6,7]. Larger symmetries can also be introduced, and a popular scenario is the one known as the Branco-Grimus-Lavoura (BGL) model [8], where a Z 3 symmetry is used. In the BGL, there are FCNCs present at treelevel; however, due to the specific symmetry implementation, they turn out to be highly suppressedà la Minimal Flavor Violation [9,10].
A complete survey of 2HDMs with an Abelian flavor symmetry, either discrete or continuous, has been done by Ferreira & Silva [11] (later confirmed through different methods in [12] and [13]). From the initially available 3 18 model implementations, and under the imposition of physical constraints in the quark sector, the authors were able to reduce this number down to 246 model candidates. In the vast majority of these implementations, the symmetry group is either the U (1) or a Z n≥3 , which from the scalar potential point of view are not distinguishable and will always lead to the presence of an accidental Goldstone boson [14]. Since many of these implementations use the global version of the symmetry, a common * astrid.ordell@thep.lu.se † roman.pasechnik@thep.lu.se ‡ hugo.serodio@thep.lu.se practice is to softly break it in the scalar potential. Another way is to promote the symmetry to a local one, turning the extra Goldstone boson into the longitudinal component of the new gauge field, i.e. the Z boson.
In the present work, we look into the general classification and phenomenology of flavored U (1) gauge symmetries in the context of 2HDMs. This was initiated as a master thesis project, and some of the details on the procedure can be found in [15]. We find that, out of the 246, only 11 distinct models survive anomaly cancellation. Out of this small subset, five of the models can be found in the literature with partial or complete phenomenological studies. The remaining scenarios were virtually unexplored until now. All these findings assumed only the SM fermionic content. Neutrinos can get their mass by a large number of different mechanisms and, therefore, deserve a detailed study on their own. The paper is organized as follows: in Sec. II, we introduce the framework, notation and the general procedure used for the classification. Section III contains a detailed implementation of all the 11 viable models, while, in Sec. IV, we present an analysis of phenomenological constraints of the new models. The results in this section serve as a first step in the phenomenological validation of such models with relatively light NP fields. We conclude and summarize our findings in Sec. VI. Additional details on the scalar potential and the anomaly conditions can be found in the Appendices.

II. ABELIAN FLAVOR SYMMETRIES IN 2HDMS
A. 2HDM Yukawa sector In this work, we consider a 2HDM with an additional gauged U (1) symmetry. There is always at least one Higgs doublet charged under U (1) ; however, if the only scale entering into its spontaneous breaking is the EW scale, it becomes extensively challenging to accommodate a relatively heavy Z gauge boson without deviating from arXiv:1909.05548v1 [hep-ph] 12 Sep 2019 the SM Z currents 1 . The particle content, in the scalar sector, is thus enlarged by an additional scalar singlet S, charged under the new Abelian symmetry.
We parameterize the scalar fields as with a = 1, 2, v a,S the vacuum expectation values (vevs) of the scalar fields and v 2 1 + v 2 2 = v 246 GeV. The presence of two Abelian symmetries allow us to rephase the scalar fields such that α 1 = α 2 = 0 and α S = 0 can be chosen without loss of generality. This choice for the phases will be used throughout this paper. Note that the inclusion of the scalar singlet, even though essential for phenomenological tests, does not effect the model classification. For details on the scalar potential, see Appendix A.
The Yukawa interactions for quarks and charged leptons are, in the flavor basis, given by withΦ a = iσ 2 Φ * a , where σ 2 is the Pauli matrix. As stated in the introduction, the fermionic content is the SM one with the fields defined in the generic flavor basis. After EW symmetry breaking, the mass matrices take the form These matrices are then diagonalized via the unitary field where D f is a diagonal matrix with the ordered masses. In the 2HDM there is a unique orthogonal combination, per fermionic sector, than can be defined and that carries all the information about flavor changing itera-tions in the mass eigenbasis, namely Any off-diagonal component in equation (5) can be further emphasized by presenting it on the following form with t β ≡ tan β = v 2 /v 1 , and where the latter term can be expressed in terms of projectors on the quark mass matrices. For further convenience, we can single out the sources of tree-level FCNCs mediated by scalar fields into the dimensionless quantities The specific form for these quantities is texture dependent and will be specified for all models in section III.

B. The new gauge sector
In this section, we introduce the relevant interactions of the new gauge field, in addition to discussing the role played by the scalar singlet. Under the corresponding Abelian symmetry, the field transformations are given by The charges of the fields are in general flavor dependent and denoted by the label of the corresponding field, except for the charge of the scalars, which are denoted by x a and x S for the doublets and singlet, respectively. We will, throughout the paper, use a compact matrix notation for the charges. For example, for q L , and similarly for all the other fields. A convenient basis to write the relevant Lagrangian for the gauge sector is the would-be SM basis, i.e. the mass eigenbasis for the Z gauge boson in the absence of mixing with the new Z boson. In such a basis, the fields and free parameters carry a hat. The neutral gauge interactions after EW symmetry breaking are then given by where ψ = {u 0 L/R , d 0 L/R , e 0 L/R , ν 0 L } and φ = {R a , I a , ρ, η}, and where we assume the gauge kinetic mixing to be zero.
Here, the second line concerns mass terms for theẐ − Z system, which in the would-be SM basis are explicitly given bŷ withf =ê/s θ c θ , whereê is the electromagnetic charge in the would-be SM basis. In the absence of mass mixing, i.e. δM 2 = 0, the gauge bosonẐ is identified with the SM Z boson and we can lose the hat notation. We would then have the MS weak mixing angle s 2 θ 0.231157(23), while the on-shell weak mixing angle, extracted by νN experiments [16], is given by 0.2237 (9). For a non-zero mass mixing, the value of s θ gets modified [17] (see Sec. IV for more details).
In the most general scenario, the mixing betweenẐ andẐ is removed through the orthogonal field transformation with such that, in the mass eigenbasis for the gauge fields, the masses of Z and Z are given by The final two lines of Eq. 10 define the interactions between the massive neutral gauge bosons and the fermionic and bosonic matter, respectively. The couplings are matrices in flavor space, explicitly given bŷ with i = u L,R , d L,R , L , e L . In the mass eigenbasis for the gauge fields, we then have that where T 3 is the isospin, Q the electric charge, and with Q ≡ T 3 + Y , where Y is the hypercharge. Note that, in the mass eigenbasis for the matter fields, Q i Z retain its form while Q i Z acquire a generic flavor structure, i.e.
In accordance with the scalar sector, all sources of treelevel FCNCs mediated by the neutral gauge bosons, are encoded in the dimensionless quantity Ξ i . Further details on the specific form of these couplings are given for all models in Sec. III.
C. General procedure for anomaly-free solutions In 2011, Ferreira & Silva classified all possible implementations of a global Abelian symmetry in the quark sector of a 2HDM [11]. In short, they used that a flavored symmetry transformation, as the one defined in Eq. 8, imposes constraints on the Yukawa couplings. In order for the Yukawa Lagrangian, in Eq. 2, to remain invariant under such field transformations, the following constraints have to be fulfilled (for the quark sector) such that and similarly for the up sector, with a = 1, 2, j = 1, 2, 3, and with no summation over repeated indices 2 . On top of these constraints, one can add three extra physical requirements: (i) No massless up-type quarks, i.e. det M u = 0 ; (ii) No massless down-type quarks, i.e. det M d = 0 ; (iii) Possibility for Dirac-type CP-violation at tree level, [11], a weaker formulation of (iii) was used, demanding the CKM mixing matrix to be non-block diagonal 3 . From these simple requirements, in addition to excluding any models equivalent up to permutations 4 , the number of 2 The case of a possible discrete Abelian symmetry was also tackled in [11]. However, those cases are of no relevance in our study. 3 For the models obtained, there is no real difference between the two formulations, but the latter has a more straightforward implementation when dealing with generic textures. 4 Any transformation that preserves the flavor symmetry, i.e. any simultaneous permutation of rows or independent permutation of columns, will result in an equivalent model. The row (column) permutations simply alter the constraints in such a way that the corresponding left (right) charges are exchanged, resulting in the procedure merely amounting to a relabelling of flavor indices.
viable texture combinations were reduced down from 3 18 to 246. In contrast to Ferreira & Silva, we now wish to classify all possible implementations of a gauged Abelian symmetry in the 2HDM, in addition to including not only the SM quarks but also the SM charged leptons. With the Abelian symmetry now being gauged, we must also make sure that the solutions are anomaly-free. The six anomaly conditions which do not cancel trivially are presented in Appendix B and involve There has been plenty of activity over the last decades on finding efficient ways of extracting generic solutions of such a system of equations [18][19][20][21] 5 .
To find all valid models, we use, as a starting point, the subset of textures in [11] corresponding to a continuous symmetry. The textures are then, in correspondence with the anomaly conditions, converted into constraints for the charges just like in Eq. 19. When including the charged leptons, there is no need to perform the whole procedure in [11] once more, as the lepton textures will get highly constrained by the anomaly equations. Hence, we only need to impose the additional physical requirement (iv) No massless charged leptons, i.e. det M e = 0.
For a 2HDM, this additional requirement can only be fulfilled when the combined texture of Π 1 and Π 2 have at least one non-zero entry in each row and each column. There are hence six possible minimal combined textures, but since they are all equivalent up to permutations, it is only of interest to consider one of them, e.g. the diagonal one. For this choice, the models have one of the following four textures (I) Π 1 : 5 There are different approaches where the anomaly constraints can be relaxed. The so-called Green-Schwarz mechanism [22][23][24] or effective anomalous U (1) scenarios [25][26][27][28], are some of the popular ones. We shall not pursue this line further and consider only models that cancel all the anomaly equations.
which can once more be turned into constraints, namely If a solution exists where all constraints are met, while simultaneously all charges are rational numbers and the quarks textures are not destroyed, the model is labeled as valid 6 . Out of the 246 models in [11], there are 116 of them corresponding to continuous symmetries, out of which only 11 non-degenerate models survive these constraints.
To identify degenerate solutions, we consider all possible permutations of the textures, i.e. every possible combination of i, j, k for the transformation where P is the 3-dimensional representation of the permutation group S 3 , such that i, j, k ranges from one to six. The permutations on the left are shared by both quark sectors, while the permutations on the right are independent. With this corresponding to an ordered sampling with replacement, there is a total of 6 3 = 216 possible permutations for each model. We should stress, once more, that besides the constraints (I)-(IV) being the minimum requirement for all leptons to gain a mass, they are also sufficient. The system of anomaly equations is so constraining for these models, that no additional textures for the charged leptons are possible 7 . The only exception to this is when all leptons, and the Higgs doublet that they couple to, have zero U (1) charge. However, as this outcome is already included as a specific charge assignment for the models in Sec. III A, the constraints (I)-(IV) are indeed sufficient. This means that we can directly identify the Yukawa entries with the charged lepton masses, and in doing so, we get six distinct implementations for each of the four cases, one for every possible assignment of the flavor indices.

III. ANOMALY-FREE MODELS
Up to permutations, there are a total of seven valid models. Four of these come in two editions, a and b, where b has the flipped texture with respect to a, i.e.

Model
Global U (1) For each model, we present the allowed textures and some analytical predictions, with the corresponding charges specified in Tab. II. An overview on which of the models that have been studied previously, in either its global or gauged form, can be found in Tab. I. As a first validity check for the less studied models, this section includes a parameter scan for fitting the quark masses and mixings (MMs). The best-fit points then serve as input values when accounting for all the remaining observables in Sec. IV.
For all models, except M5 and M6, the free parameters in the Yukawa sector are: • Yukawa couplings modulus ∈ [10 −5 , 5]; • Yukawa couplings argument ∈ [10 −10 , 2π[; The scan was performed by giving, as input, 300 different tan β values evenly distributed in logarithmic scale. For the models M5 and M6, on the other hand, it is more convenient to parameterize the Yukawa couplings with the physical quark masses and CKM mixing angles (see respective model subsection for the explicit parameterization). For the fermions masses and CKM mixing matrix data, we have used [16]. The charged lepton mass matrix is diagonal in all models and, therefore, the extraction of their Yukawa couplings is straightforward.

A. Leptons that couple exclusively to one Higgs
In this section, we present all models fulfilling the constraints specified in section II C, while simultaneously having the leptons coupling to only one of the Higgs doublets. In all these models, there exists a particular charge assignment for which one of the scalar doublets, as well as all charged leptons, are uncharged under U (1) . In this special case, all models can fulfill the anomaly constraints even without the introduction of charged leptons.

Model M1
Let us being with the one model that is naturally flavor conserving, namely with the corresponding charges specified in Tab. II. As familiar, this model has no tree-level FCNCs, since K d,u = 0 and Ξ i are diagonal.

Model M2
Next, we have the textures By rephasings of the left and right-handed quarks, we can select, at most, real values for eight of the ten Yukawa couplings. One viable choice is to use real values for all couplings apart from the ones in Γ 2 . With this, we have a reduction from 20 real parameters, down to 12, in the quark sector. For the lepton sector, there are, for all models in this paper, enough degrees of freedom to remove the complex phase for all couplings.
From Eq. 24, we see that there are FCNCs in all of the quark sectors, as neither of the sets are Abelian (proof given in [43]). Nevertheless, the small number of free parameters allows us to extract some additional information. Using the principal invariants for H u = M u M † u , we can reduce the number of free parameters to two. The relatively compact solutions for the up-quark sector read with a and b defined as that remain free.
The sources of FCNCs in the gauge sector are given by ), with x = y and with the projectors in the new basis defined as in [39], namely where V is the CKM matrix, V = U † uL U dL and with the projectors defined as (P a ) ij ≡ δ ij δ ja , with no summation over repeated indices and with i, j, a taking on values from one to three. For the scalar sector, the relevant sources of FCNCs can be expressed as In Fig. 1, we plot the preferred magnitudes for elements (M u ) 22 and (M u ) 21 when fitting MM. The red parameter points correspond to a region of parameter space with all deviations below 1σ, while blue points have a deviation below 3σ. In the plot, we see a strong preference towards either (M u ) 22 or (M u ) 21 being of the order of the top mass. Similar features can be found in the down sector when plotting (M d ) 23 as a function of (M d ) 22 . Red crosses correspond to regions with all deviations below 1σ, and blue to deviations below 3σ.

Model M3a and M3b
For model M3a, the textures are given by For this model, there are again FCNCs in all quark sectors, but this time around, there are further restrictions. With the charges of the first two generations being degenerate for all quarks, i.e. q 1 = q 2 , u 1 = u 2 and d 1 = d 2 , we have that with x = y. While, for the scalar-mediated FCNCs, Model M3b, on the other hand, has the same expression for Ξ as model M3a, while the scalar-mediated FC-NCs are now given by with q = u, d. Figure 2 shows the preferred values for elements of the mass matrices. Here we see that for model M3b, (M u ) 33 prefers having values around the up mass, while the remaining elements have a larger flexibility. The equivalent plot for model M3a, shows a strong preference for (M d ) 33 to be of the order of the down quark mass.
The figure is a so-called boxplot, where the array of data is sorted by magnitude and then split into four equally sized parts, each referred to as a quartile. The black line within the box is the median of the data set, while the box itself extends from the first to the third quartile. The so-called whiskers span from the smallest to the largest values in the set, out of the points not classified as outliers. An outlier is here defined as any point further away from the first or third quartile than one and a half times the length of the box. For a normal distribution, this would correspond to 0.7 percent of the data. The outliers are marked as black dots in the figure and, due to the low number of good points in their vicinity, they indicate the fine-tuned regions of the parameter space. For both models M3a and M3b, at most eight out of the fourteen Yukawa couplings can be set real from rephasings of the left and right-handed quarks. There are hence 20 parameters in the quark sector, in addition to the three parameters in the lepton sector. One possible choice is for all couplings in Γ 2 and ∆ 2 to be real, in addition to the three diagonal couplings in ∆ 1 and (Γ 1 ) 33 .

B. Leptons that couple to both Higgs doublets
This section contains all models fulfilling the constraints specified in section II C, while simultaneously having leptons that couple to both Higgs doublets. Contrary to the models discussed in the previous subsection, the inclusion of charged leptons carrying U (1) charge is now crucial for all charge assignments.

Model M4a and M4b
For model M4a, the textures are given by with the corresponding charges specified in Tab. II. Model M4a can have at most eight simultaneously real Yukawa couplings from rephasings of the quark fields. There are hence 12 parameters in the quark sector, in addition to the three parameters in the lepton sector.
One allowed scenario is for all couplings in Γ 1 , ∆ 1 and ∆ 2 to be real, in addition to (Γ 2 ) 11 (similar for M4b).
Since both combinations M u M † u and M † u M u are real and block diagonal, the diagonalization matrices U uL and U uR are simply given by a real orthogonal rotation. Using the same notation as in Eq. 12, the rotation angle for U uL will be given by and the same for tan 2θ uR , but with (M u ) 23 and (M u ) 32 exchanged in all places. This can be further simplified by once again using the principal invariants, which has three possible solutions in the up sector, namely and two equivalent expressions but with the masses exchanged as m c ↔ m t and m c ↔ m u . Comparing this with the result of the scan in Fig. 3, we see, as expected, that (M u ) 11 always take on one of the three possible uptype quark masses. The gauge mediated FCNC couplings are then given by with x = y and the scalar mediated FCNCs by For model M4b, we have a similar scenario, but with the roles of the up and down sectors exchanged.

Model M5
Next, we have a so-called "right model of type E" [39], with the textures Similarly, for the scalar sector, we have As stated at the beginning of this section, for model M5 we use the physical quark masses and mixing matrix as inputs. Since the left-handed quark doublet charges are fully degenerate under U (1) , any weak basis transformation (WBT) in this sector will never affect any observable. Therefore, applying the WBT q 0 L → U uL q 0 L , we are instead in the basis The only free quantities are the right-handed rotation matrices. The Yukawa couplings get, therefore, fully determined by the above mass matrices.

Model M6a and M6b
Model M6a is the familiar BGL model [8], with the corresponding charges specified in Tab. II, and with FCNCs only in the left-handed down sector. As is well-known for BGL-type models, U uL and U uR are block diagonal. In the right-handed down sector, all charges are degenerate, leading to a diagonal Ξ dR . For the other three quark sectors, on the other hand, there are two degenerate charges yielding with the projector index a = 1, 2, 3 for when (M u ) 33 = m u , m c , m t , respectively. For the scalar mediated FCNCs we have with a defined as before. Model M6b obeys similar relations, but with the role of the up and down sectors exchanged. Here the FCNCs are instead limited to the left-handed up sector. For models M6a and M6b, there are hence no new parameters.
Here, we can again choose a basis in which the expressions for the mass matrices are simplified. Besides the WBT used for model M5, we can now transform also the right-handed quarks such that in the new basis without influencing the conserved currents, due to the way the degeneracy of charges coincides with the block diagonality of the right-handed mixing matrices.

Model M7a and M7b
Finally, model M7a has the textures with the corresponding charges specified in Tab. II and with FCNCs in all quark sectors. Here, there are two degenerate charges in all quark sectors except for the right-handed down quarks, such that the FCNC sources in the gauge sector are given by with x = y. While for the scalar-sector we get v 2 ).
Model M7b obeys similar relations, but with the down and up sector exchanged. In Fig. 4, we show the preferred values for entries of the up quark mass matrices, where red crosses correspond to regions with all deviations below 1σ, and blue points to a deviation below 3σ. Similar features can be found in the down sector of M7a. For model M7a and M7b, there are, at most, nine out of fourteen Yukawa couplings that can be set real from rephasings of the quark fields. As such, there are 19 parameters in the quark sector. One allowed selection of real entries correspond to all couplings in Γ 2 and ∆ 2 , in addition to (∆ 1 ) 12 and (∆ 1 ) 31 . 33 and (M u ) 31 in model M7b. Here, red crosses correspond to parameter points with all deviations smaller than 1σ, and blue ones to all deviations below 3σ.

C. Models tendencies
The overall tendencies for fitting the quark MMs for the models are gathered in Tab. III. Here, the second column displays the sum of all χ 2 for the very best point for TABLE III. Display of the effort required for fitting MMs for the various models. The second column presents the total χ 2 of the best-fit point, the third column shows the percentage of points with a χ 2 below one, and the fourth column gives the observables with the most significant deviations from data, for the points with a deviation in the range 1σ − 5σ.
each model, where we see that all models contain regions in parameter space with an excellent agreement with the quarks MMs. The third column instead presents the percentage of points with an individual χ 2 below one, from which we can tell that model M3a has the most preferable landscape, while model M4a has the least preferable one. And, finally, the fourth column presents the two most difficult observables to fit. For most models, this is the strange mass. Note that the largest individual deviation is still always below 5σ.

IV. SCAN AND PHENOMENOLOGICAL TESTS
In this section, we present the details on the scan performed for each of the models found in Sec. III. For all models except M5 and M6, we have given as input the values of the best-fit points form the masses and mixing scan. We are then left to range over: • Scalar singlet vev v S ∈ [10 2 , 10 6 ] GeV; • Dimensionless scalar parameters λ 12,S,Si ∈ [−1, 1]; • Dimensionful scalar parameter a 1 ∈ [ −5, 5] TeV.
The mass parameters m 2 i and m 2 S are fixed by the tadpole conditions, while s θ and f are left as free parameters (for more details, see Sec. IV A). In the scan, we compute the masses for the neutral gauge bosons, neutral and charged scalars and their respective couplings to fermions and bosons. With that information, we then compute their corrections to a set of physical observables detailed below.
For models M5 and M6, the approach is slightly different. As explained in the previous section, there is no need to perform a masses and mixing scan since these parameters can be given as input. In model M6, the Yukawa sector is fully determined by the known quark masses and CKM mixing, such that the only new free parameter is tan β, for which we choose β to be in the interval [10 −3 , 1.56]. In model M5, we can once more set the quarks' masses and mixing to the physical ones, but we are still left with two generic unitary matrices in the up and down right-handed quark sectors. These unitary matrices are parameterized as V R = KR 12 R 13 R 23 , with K a diagonal matrix with 3 phases bounded by [10 −1 , 2π[ and R ij a complex rotation in the ij plane with both the angle and phase constrained to the interval [10 −5 , 2π[. Next, we give a more detailed description of the observables used in the scan.

A. Electroweak and low energy constraints
As already discussed in Sec. II, the presence of mixing between the "would be" SMẐ and the new gauge boson Z will in turn lead to modifications on the currents of the physical Z boson. The electroweak sector of the SM has been probed with a very high level of accuracy, such that only a very small mixing is in general allowed. In our study of the EW constraints, we look at the Z-pole pseudo-observables To compute these observables, we follow closely the LEP-TOP [44,45] approach, which includes one-loop EW corrections. In this approachᾱ, m Z and G F are given as inputs, where G F is the Fermi coupling constant. However, since the presence of the new Z boson alters the Z mass, we include m Z andᾱ in the fit. We then take as the defining relation for s θ , whereᾱ includes the running due to the three lepton and five "light" quark loops leading toᾱ −1 ≡ α(m 2 Z ) −1 = 128.878 ± 0.090 [45]. The analytical definitions for these pseudo-observables are the standard ones, and their one-loop corrections can be found in [44,45]. The effects of the scalars were taken into account through the computation of the oblique parameters [46], along the same lines as Burgess [47].
We also compute the off-pole cross sections with both Z and Z , for the final states µ + µ − , τ + τ − and qq, within the energy range √ s = [130, 207] GeV given by LEP-II [48].
The previous observables do not really constrain the interactions of the Z boson with the top quark. We have therefore included also the top width in the fit, in addition to the rare top decays t → Zq and t → hq.
Low energy constraints such as atomic parity violation, electric dipole moments (EDMs) and the muon magnetic moment, can place significant bounds on the parameter space of our scenarios. Atomic parity violation for 133 55 Cs and 205 81 Tl have been included in the fit [49,50]. For the neutron EDM, we have considered one-loop corrections from neutral and charged scalars [51] and Barr-Zee diagrams [52,53], while for the muon magnetic moment both the gauge and scalar one-loop contributions are taken into account [54]. The EDM of the electron, on the other hand, is virtually zero since there are no FCNCs in the lepton sector.

B. Meson constraints
A flavor changing mediator, scalar or vector, is strongly constrained from meson observables. In the fit we have included: The most dominant bounds usually come from the ∆F = 2 transitions observables, ∆M M 0 and K . The SM contribution for the K 0 and B 0 q systems is dominated by the top-loop box diagrams, while the SM prediction for the D 0 system is less trivial. The short-distance effects are virtually zero, but long-distance contributions can be sizable (at the level of the current experimental values). However, long-distance effects are plagued with significant hadronic uncertainties, making it difficult to estimate the SM contribution. For this, we follow the approach of [58,59] and attribute the experimental values to short-distance NP. This will serve as an upperbound, such that NP contributions do not greatly exceed the current experimental bound. Box diagrams for the charged Higgs were also included in the numerical scan, see [5] for the general expressions. The hadronic matrix elements for ∆F = 2 processes are very important and were included together with the QCD running of the Wilson coefficients from the NP scale down to the meson scale, this was done in line with refs. [55][56][57]. We took a conservative approach for the meson oscillation observables and for the Kaon sector, with the error for ∆M K taken as 20% of the central value [16] while for | K | we take the error to be 10%. We also use a 20% error for the D-meson sector. For the B-meson sector, since it is less sensitive to long-distance contributions, we took the error for the ∆M Bq to be on the 5% level.
The b → s transitions are very relevant for these classes of models and have gotten the attention of a large community since the 2014 LHCb results [78] seemed to indicate a breaking of lepton universality. There have been several global fit analyses [79][80][81][82], pointing the preferred region for some of the four-Fermi effective operators. The experimental status is very likely to change in the next couple of years. Currently LHCb reports [83] suggesting a 2σ − 3σ deviation from the SM, given [84]. The Belle collaboration also measured this observable and gives values more compatible with the SM [85], i.e.
To compute the contributions of NP to these observables, we start by defining the effective Hamiltonian with Q ( ) For the SM Wilson coefficients, we have C SM 9 −C SM 10 4.2 and zero for their primed counterparts. The R K and R K * ratios are then given by [87] R K/K * 1 + ∆ +/− + Σ +/− , with ∆ +/− the SM-NP interference term and Σ +/− the pure NP contribution. Their explicit form is given by ∆ ± 0.24 Re C µµ,NP 9−10 ± C µµ,NP and where we defined C ( ) 9−10 ≡ C ( ) 9 − C ( ) 10 . These Wilson coefficients, in our framework, take the form with λ Bs t ≡ V tb V * ts and Z 0 = {Z, Z }, and equivalent for C 10 , but with the sign in front Q e L Z 0 flipped. The corresponding results are discussed in Sec. V C.
The rare decay modes B → {K, K * , X s }νν can be important tests for NP. The SM short-distance predictions were computed in [88] and are still an order of magnitude below the experimental bounds [89][90][91], leaving plenty of room for NP. We can parameterize the NP effects through the ratios with the SM predictions, see ref. [88].

C. Collider constraints
At LHC several direct searches for new neutral gauge and scalar bosons have been performed. As we shall see in the next section, due to the necessity of a large vev for the scalar singlet S, the new scalar fields tend to be heavier than the Z and, in most of the models, out of reach for the current colliders. We have therefore only looked at the collider constraints coming from direct searches of Z gauge boson.
We follow a simplified approach by working in the narrow width approximation (NWA). For all models, this turns out to be a good approximation in the low m Z range, i.e. below a few tens of TeV. Under the NWA we split the cross section under production and decay, i.e.
where the production cross section is given by [92,93] σ(pp → Z ) = 4 3 with q = u, c, d, s, b. The integration over the parton distribution functions is included in the definition of ω q (s, m 2 Z ) given in [94], these were computed using the PDF set NNPDF23 nnlo as 0119 qed [95].
The experimental studies used in this analysis are presented in Tab. IV.

V. MODEL CONSTRAINTS AND PREDICTIONS
In this section, we study all models except for M1, as it has no FCNCs and has already been extensively studied with the corresponding charges given in Tab. II. Some general results for these charges are gathered in Fig. 5, Tab. V and Tab. VI. Figure 5 shows the ability of the various models to fit all observables. Here we see that, for the explored regions of parameter space, Model M2, M3a, M3b, M4a, M5, M7a and M7b can all reach low deviations with data. In particular, model M2, M4a and M7b do so without the corresponding parameter points being outliers of the distributions. Surprisingly enough, the BGL model lacks parameter points below ∼ 3.9σ.
This could be understood from the fact that, while the FCNCs are naturally suppressed in the BGL model, they are also fixed, which reduces the flexibility when fitting parameters. Table V displays again the effort required for fitting the various models. Here, we show points with an individual deviation below 5σ and m Z < 30 TeV. The second column contains the number of parameter points, out of the initial 300, that survived this requirement when fitting the MMs. The MM points were then used as input to the minimization of the remaining physical observables, with the number of surviving parameter points shown in the third column. Many of these points share the same MM input but differ in values for the remaining free parameters of the model. Next, the fourth column displays the most constraining collider channels, out of the ones in Tab. IV. The corresponding channel is only displayed whenever it actually excludes any of the parameter points from the scan, for that corresponding model. Overall, collider constraints are dominated by the µ + µ − channel. We recall that for all model implementations chosen, the electron carries no charge under U (1) . Finally, the fifth column identifies the observables with the largest individual χ 2 , with the percentages reflecting their relative dominance. For example, for model M4b, 78% of the FIT points had their largest deviation coming from the meson observable ∆M K .  While a more complete scanning of the parameters space would be desirable, the results in Tab. V already show some of the model tendencies. For almost all models, ∆M K is the observable with the highest deviation. However, since the meson sector is plagued with uncertainties on hadronic matrix elements and long-distance contributions, models with deviations around 4σ should not be completely discarded. Such sizable deviations in the EW sector are more difficult to accommodate making, for example, model M4a less appealing.
In general, all the models considered require a vast hierarchy between the EW scale and the vev of the scalar singlet, resulting in the decoupling of the neutral scalar singlet (H 0 3 ), and a full doublet of scalars (H 0 1 , H 0 2 and H + ). In Tab. VI, we have presented the lowest possible scalar masses obtained in the scan for the various models. Here, we see that, even for the parameter points with the lowest scalar masses, the scale hierarchies are evident. Only models M3b, M5, M6a, M6b and M7a have scalar masses below 4 TeV. These are, however, quite isolated  points, as for all models except M6, this amounts to less than 1% of the FIT points in Tab V. For model M6a and M6b, however, v S can be of the order of a few TeV, allowing the scalar masses to come down to a few hundred GeV. For M6a these amount to around 3% of the FIT points and 10% for model M6b.
As already mentioned, we have only included the direct searches on Z , such that similar searches for scalars could easily rule out these low mass points. Therefore, Tab. VI should only serve as an indication for which models we might expect relatively light scalars. Also, the results presented here correspond to the charge assignments given in Eq. 57. A different choice could help lower the bound for the scalar masses. For example, in model M5, there were many allowed points with scalar masses in the range [500, 10 3 ] GeV. All these were only excluded once the collider searches were included. A different coupling to light quarks could then easily have revived these points.

A. Collider bounds
In Fig. 6, we plot g as a function of the Z mass, for four of the models. Every parameter point excluded by collider constraints have been removed from the plot, and instead replaced by a red, shaded area, for all larger values of g with the same Z mass. As such, the red, shaded area, corresponds roughly to the area excluded by collider constraints. This area should not be seen as a forbidden region of the parameter space, but more as a less favorable one, due to collider constraints. For example, in model M5, there are still several points inside this shaded red area that survived the collider bounds. These correspond to cases where the Z coupling strength to light quarks is small, leading to a reduction in the production cross section. The blue, shaded area, with a larger grid, indicates any regions of parameter space not probed by our scans. From Eqs. 11 and 14 we see that, in the case of a small mixing angle, θ M , and a large hierarchy between the vev of the scalar singlet and EW scale, the Z mass is approximately given by m Z ∼ g v S x S , explaining the mostly linear behavior of the plots. This also means that the tilt of the plot corresponds roughly to the inverse of the product of the U (1) charge of the scalar singlet, and its vev, i.e. (v S x S ) −1 .
Given the charge assignments in Eq. 57, the ratio Γ Z /m Z was always below the 5% level for all models, except for M5 (7%) and M3b (13%). When g is large and m Z is small, we can have a Landau pole appearing at lower scales. A simple estimation, using the 1-loop renormalization group equation for the evolution of g , gives, for the point (g , m Z ) = (1, 10 TeV), a Landau pole at 10 9 GeV for M3b and 10 17 GeV for M5. Once g is lower than 0.5, all models have a Landau pole well above the GUT scale.
Overall, models M2 and M5 show the most promising features. For model M2, the collider constraints do not even exclude any of the points. In comparison, the BGL model (not included in the plot), looks similar to model M5, but without reaching equally low χ 2 values.

B. Ternary plots
In Fig. 7-9, we plot the preferred values for the U (1) charges in the mass eigenbasis of the gauge fields, Q i Z /g . The colored regions have a deviation below 5σ for all observables 8 , with the red points corresponding to a Z − Z mixing below one percent and the blue points to a mixing above one percent. As such, for the red points, the charge matrix in the mass and flavor eigenbases are approximately equal, i.e.
In each plot, the black stars mark the U (1) charges in the flavor eigenbasis of the quarks. Hence, for any red parameter point positioned in the vicinity of a star, we have that Ξ = U † XU ∼ X, i.e. that the corresponding base changing matrix is close to unity. With this, the ternary plots can be used to visualize the amount of FCNCs that occur in a particular sector.
With the CKM matrix given by U † uL U dL , the mixing matrices for the left-handed sectors are always rather close to unity. Hence, we will only include plots for the right-handed sectors. Take for example model M2, with x = 1/3, y = 0 and z = 1. For these values, the charges for all quarks in the flavor eigenbasis are given by such that there are six black stars in Fig. 7, namely all permutations of {1/3, 0, −1/3}. With the red points deviating by a small amount from the black stars, some off-diagonal contributions are allowed. For model M3a, there are two degenerate charges in all quark sectors. As such, there are three black stars in Fig. 8, rather than six. Here, the right-handed up sector shows a significant deviation, while for the down sector, the parameter points are practically on top of one of the black stars. Similar features can be found in model M3b, but with the FCNCs being slightly subdued in the up sector, and with some off-diagonal elements also in the down sector.
In Fig. 9, we have the preferred U (1) charges for model M7a. Here, there are sizable off-diagonal elements in both the right-handed sectors. For model M7b, the base changing matrices are slightly closer to unity.
Overall, we see that the ternary plots are dominated by red parameter points, i.e. there is a strong preference toward a small Z − Z mixing.

C. Flavor anomalies
The values for the ratios of C NP( ) 9,10 , with the charge assignments in Eq. 57, are shown in Tab. VII. Given these charges, and the identification of the electron as the charged lepton with no U (1) charge, the Wilson coefficients C NPe( ) 9,10 are automatically zero for all models. It turns out that, for this choice, there is always a scenario where the muon couples vectorially to the Z , leaving C NPµ( ) 9 as the only non-zero Wilson coefficient. In Tab. VII, these cases are indicated by gray, shaded cells, which are also the scenarios used in the scan. Note that, for models M5 and M6b (and M1), the unprimed Wilson coefficients are zero, as there are no FCNCs in the left-handed down sector. Similarly, the primed Wilson coefficients are zero for models M6a and M6b (and M1), as they have no FCNCs in the right-handed down sector.
For most models, the predictions for R K and R K * where just the SM ones. However, for model M5, the presence of C NPµ 9 was enough to have sizable deviations from the SM prediction. In Fig. 10, we plot the ratio for R * K , R K for model M5, together with the 1σ bounds for the LCHb and Belle measurements, as specified in Sec. IV B. Here we see that there are parameter points that can account for the current anomaly. Note that, in the plot, we are showing the values only in the region of interest for the anomaly -there exist plenty of points in agreement with the SM prediction and also with a larger  . Here, the asterisk indicates entries for which C NPµ( ) 9 is zero, while shaded cells corresponds to the flavor configurations used in the scans.
value for R K than 1. We have also looked into the predictions of these models for the ratios R ν K,K * as defined in [88]. Most model predictions were compatible with the SM results. However, models M4a and M5 allowed some significant deviations. For model M4a (R ν K , R ν K * ) varied from (0.8, 1.1) to (1.3, 0.85) linearly. In model M5, since there are only right-handed FCNCs, the two ratios are equal and were in the range [0.95, 1.45].

VI. CONCLUSIONS
This work has two central results. One is the classification of all anomaly-free 2HDMs with a gauged flavor symmetry, and the other is the analysis of their phenomenological consistency. Together, the two components can be used to identify the current research gap for this type of 2HDMs and work as a starting point for further studies.
Out of the 116 model candidates corresponding to continuous symmetries in [11], only 11 models survived after imposing anomaly cancellation. Six of these models were studied here for the first time. In the phenomenological study, each model was subjected to a range of observational constraints, and while there were good parameter points for all 11 models, we wish to give extra attention to models M2 and M5. For these two models, one could easily find extended regions of parameter space with a deviation below 3σ, in addition to model M5 being able to accommodate the R K , R * K anomaly with only the presence of C NPµ 9 . As such, they deserve further consideration.