Anomaly-free 2HDMs with a gauged abelian symmetry and two generations of right-handed neutrinos

Despite the popularity of two-Higgs-doublet models (2HDMs) with a gauged abelian flavor symmetry, the allowed Yukawa textures have only been partly explored so far. In this work, we classify and compare every anomaly-free instance, in the case of having two-generations of right-handed neutrinos and a type-I seesaw mechanism. We found 16 valid implementations in total, out of which 11 agree well with current experimental bounds. To our knowledge, neither of these models have been considered previously.


I. INTRODUCTION
Despite there being relatively strict bounds on Higgs couplings to the SM gauge bosons and heavy fermions [1], the Higgs sector remains among the least constrained in the SM. The simplest example of a non-minimal Higgs sector is built by adding one extra Higgs doublet and is known as the 2HDM (for a detailed review, see e.g. Refs. [2][3][4]).
The 2HDM was first proposed by T. D. Lee, over 40 years ago, to motivate spontaneous CP violation [5] and has evolved into one of the most well-explored extensions of the SM today. A generic 2HDM (i.e. one without any imposed symmetries) is, however, not very enlightening, as its scalar potential and Yukawa sector contains a large number of free parameters. Besides, the generic 2HDM features flavor-changing neutral currents (FCNCs) at tree-level, which are heavily constrained by experiments [6].
By imposing a continuous or discrete abelian symmetry, one can convert a 2HDM into a more predictive framework, with the possible realizations for the quark sector presented in Ref. [7] and confirmed in Ref. [8,9]. One particularly well-explored path is to look for models that avoid tree-level FCNCs by construction. The simplest example of this is when all right-handed fermions of the same charge couple to only one of the Higgs doublets, by the implementation of a discrete Z 2 symmetry, and are referred to as being naturally flavor conserving (NFC) [10,11]. Another, softer, method of avoiding the tree-level FCNCs, is by imposing the Yukawa couplings of the two scalar doublets to align in flavor space, as done in the so-called flavor-aligned 2HDM [12][13][14].
Rather than entirely forbidding tree-level FCNCs, an alternative procedure is to suppress them sufficiently to comply with experimental data. A particularly simple and important scenario of such a minimally flavor violating (MFV) scheme is realized in Branco, Grimus and Lavoura's (BGL) model [15,16]. In Ref. [17] such a BGL * astrid.ordell@thep.lu.se † roman.pasechnik@thep.lu.se ‡ hugo.serodio@thep.lu.se formulation was extended to incorporate a mechanism for generating neutrino masses, while a gauged version of the BGL model was implemented for the first time in Ref. [18]. Due to their simplicity and elegance, most of the available models on the market today realize either NFC or MFV scenarios. While this is a reasonable outcome, their dominance has led to a significant research gap in terms of possible Yukawa textures. For example, as noted in our earlier work [19], more than half of all anomaly-free implementations of a 2HDM with a gauged U (1) flavor symmetry had never been previously considered.
In this paper, we again classify and compare all anomaly-free implementations of a 2HDM with a gauged abelian flavor symmetry, but this time around also incorporating a mechanism for generation of neutrino masses. In particular, we consider a type-I seesaw mechanism [20][21][22][23] with two generations of right-handed neutrinos [24][25][26][27][28][29]. 1 This corresponds to having two massive and one massless generation of left-handed neutrinos, which is the minimal possible setting in agreement with current experimental bounds [1]. From this, we found that 11 models, out of the total 16, were in good agreement with data, with three of them being particularly promising.
The paper is organized as follows: In Sec. II and III, we introduce our model and the method of finding anomalyfree solutions, respectively. Next, in Sec. IV, we present the explicit form of all anomaly-free implementations, and, finally, Sec. V, VI and VII involves the scan, its results and our conclusions.

II. MODEL SPECIFICS
In this paper, we consider models where the SM particle content is extended by two generations of righthanded neutrinos, a Higgs doublet Φ, a Z gauge boson and a scalar singlet S. In general, all matter-sectors (SM and non-SM) are charged under the new U (1) symmetry and transform as where X χ is the corresponding U (1) charge, which is in general flavor dependent, where j is a flavor index, a = 1, 2, and where the zero superscript denotes the flavor basis.
While a subset of the charges can be set to zero, we always require that the singlet, and at least one of the Higgs doublets, are charged under U (1) . As such, the scalar masses are not solely determined by the EW scale and by the dimensionless couplings in the scalar potential (bounded by unitarity), but can be adjusted via the VEV of the scalar singlet.
In the flavor basis, the Yukawa interactions are then given by with the scalar fields parameterized as where v S and v a , a = 1, 2, are the VEVs of the scalar singlet and the Higgs doublets, respectively. With the additional abelian symmetry, we can set two of the VEVs real, e.g. α 1 = α 2 = 0. After EW symmetry breaking, the mass matrices for the quarks and charged leptons are given by and equivalently for M d and M e , but with ∆ exchanged for Γ and Π, respectively. For the neutrino sector, on the other hand, we have with In the limit of M R m D , the light eigenstates are the left-handed Majorana particles, such that the effective mass term is given by As the heavy eigenstates decouple, the scalar potential and the new gauge sector have the same form as in Ref. [19].

III. PROCEDURE FOR FINDING ANOMALY-FREE SOLUTIONS
For two generations of massive neutrinos, we have the following physical requirements in the leptonic sector (i) No massless charged leptons, i.e. det M e = 0; (ii) Two generations of massive neutrinos, i.e. rank M ν = 2; as apparent from Eq. (7). With the neutrinos gaining their mass via a type-I seesaw mechanism, condition number (ii) corresponds to M R being a 2 × 2 matrix with a non-zero determinant and m D being a 3 × 2 matrix with rank 2.
The key is to find an efficient method for identifying all textures that fulfill these constraints, as naively considering every possible combination would require an immense amount of computational power. To this end, we will use the procedure introduced in Ref. [19], i.e. using a set of minimal constraints, in addition to avoiding textures that are equivalent up to permutations of flavor indices.

A. The charged lepton sector
Let us begin with the charged lepton sector, which is the only sector with the same minimal textures as in Ref. [19]. Here, the requirement of a non-zero determinant forces the combined texture of Π 1 and Π 2 to have at least one non-zero entry in each row and each column. In other words, there are six possible textures for M eall being permutations of the diagonal texture As we have not yet introduced the sectors involving Dirac neutrinos, Majorana neutrinos and quarks, permutations of row and columns simply amount to a relabelling of flavor indices. Hence, we only need to consider one of the possible six textures, e.g. the diagonal one presented in Eq. (8). This leaves four possibilities for the minimal textures of Π 1 and Π 2 , namely Note that, as these are so-called minimal textures, additional non-zero entries are allowed. The textures are then translated into a set of linear equations for the U (1) charges, since, in order for the abelian flavor symmetry in Eq. (1) to be a symmetry of the Lagrangian, it is required that which corresponds to where i, j are the flavor indices running from one to three, and where there is no summation over the repeated indices. When converting the textures in Eq. (9), we exclude constraints on the form X i −X ej = X Φa , as we are want to allow for additional non-zero entries. As such, the textures correspond to Note that, in excluding constraints on the form X i − X ej = X Φa , the conditions in Eq. (12) are sufficient for generating all viable lepton textures.

B. The Majorana neutrinos
For the Majorana neutrinos, there are nine possible minimal textures for A, B and C that fulfill the constraints of M R being symmetric and having a non-zero determinant. Once again, we do not need to consider any permutations of rows and columns, as these permutations are separate from the ones in the charged lepton sector, and hence simply amount to a relabelling of indices. The first three minimal textures are given by which corresponds to the constraints X ν1 + X ν2 = 0, X ν1 + X ν2 = −X S and X ν1 + X ν2 = X S , respectively, since invariance under the flavor symmetry requires that Again, there are no constraints in the form of 2X ν 1(2) = 0 or 2X ν 1(2) = ±X S and, as such, the final textures are not limited to the ones in Eq. 13. Next, we have the remaining six minimal textures for the Majorana neutrinos, namely 1 : where {1, 2, 3} can be any permutation of {A, B, C}. The corresponding constraints are given by Furthermore, with the phase sensitive part of the scalar potential given by one of the following four terms, invariance requires for one the following four constraints to be met For more details, see Ref. [19]. Note that, in cases where only A has a non-zero texture, all four solutions are allowed. C. The Dirac neutrinos Next, we have the possible textures for Dirac neutrinos, where the permutations of rows are not independent of those in the charged lepton sector, and where permutation of columns are not separate from those in the Majorana neutrino sector. As such, we must consider every possible permutation explicitly, which for a rank 2 m D corresponds to a total of 6 possible textures 1 : where each matrix comes in four versions, depending on whether the non-zero texture sits in Σ 1 or Σ 2 . For example, for texture number 1 we have the four possibilities with the corresponding constraints given by as invariance under the flavor symmetry requires that

D. The quark sector
For the quark sector, on the other hand, we can use the textures from Ref. [7] that corresponds to continuous symmetries, rather than repeating the procedure with minimal constraints. The quark textures in Ref. [7] are then translated into a set of linear constraints using the same method as above, i.e. by demanding invariance under the flavor symmetry However, this time around, we do make use of the non-equalities, i.e. replacing the zero textures with X qi − X dj = X Φa and X qi − X uj = −X Φa , respectively, since the result in Ref. [7] corresponds to complete textures and not to minimal ones.

E. Combining all sectors
When combining the constraints from all sectors above, in addition to including the anomaly cancellation conditions involving U (1) , 2 we end up with one large system of equations for the U (1) charges. In total, there are linear equations coming from the charged lepton textures, Majorana neutrino textures, Dirac neutrino textures, quark textures and from the invariance of the scalar potential, in addition to a mix of linear and non-linear constraints coming from the anomaly constraints. Any rational solution to this system is then labeled as a valid model. Note that, as we loop over every possible combination of constraints, we are guaranteed to find all valid solutions, in comparison to e.g. a scanning procedure, in which one might only find a subset of them.
To avoid degenerate solutions once everything is combined, we make sure that neither of the models can be reached from one another via the transformations where P is the three-dimensional representation of the permutation group S 3 , P is the two-dimensional representation of S 2 and where G = A, B, C. That is, i, j, k, l, m ranges from one to six, while n ranges from one to two.

IV. ANOMALY-FREE MODELS
Up to the implementation of Eq. (18) there are a total of 16 valid models. Below they are gathered into five separate categories, Category A to E, based on whether they share the same textures. The corresponding charges for each model are presented in Tab. II.
Note that three of the Majorana-and Dirac neutrino textures below (in E3, E4 and E5) have been studied previously in Ref. [24], but in the context of having only the SM particle content.
A. Category A Category A consists of two lepton textures and four quark textures, where all eight combinations correspond to valid models, model A1 to A8. The two lepton textures are given by and the four quark textures by (Q1) : Note that lepton texture E1 comes in one additional version, with the textures of B and C exchanged, and X S replaced with −X S , while lepton texture E6 allows for any of the four solutions in Eq. (18).
The eight possible combinations of the lepton and quark textures are then referred to as Model A1 to A4 : E1 with Q1 to Q4, Model A5 to A8 : E6 with Q1 to Q4.
Here, each lepton texture has a maximum number of parameters that can be made real from rephasings of the charged leptons, Majorana neutrinos and Dirac neutrinos. One possible choice, which allows for the maximum number of real parameters, is shown directly in the textures above, where asterisks corresponds to complex entries and crosses to real entries. For example, lepton texture E1 has a maximum of nine parameters that can be made real simultaneously. Note that this notation applies for the entirety of Sec. IV, but not in any of the remaining sections.
Similarly, each quark texture has a maximum number of Yukawa couplings that can made real by rephasings of the right-and left-handed quark fields.

B. Category B
For Category B, there are two lepton textures and two quark textures, with all four combinations corresponding to valid models, model B1 to B4. Here, the lepton textures are given by (E2) : and the two quark textures by (Q5) : Note also that both lepton textures come in an additional version, with X S replaced for −X S and with the textures of B and C exchanged. The same thing goes for the lepton texture in Category C and D.
C. Category C For Category C, there is one lepton texture and two quark textures, with both combinations corresponding to valid models, model C1 and C2. Here, the one lepton texture is given by : : with the two possible combinations defined as Model D1 (D2) : E5 with Q9 (Q10).

V. SCANNING PROCEDURE AND PHENOMENOLOGICAL TESTS
The model scan is divided into four subsequent parts: the scalar potential scan, the leptonic scan, the quark scan and the full scan. This separation is done to increase its efficiency -for example, as shown in the result section, half of the models can be excluded already after the leptonic scan. Below follows a summary of each stage.
In this section, we evaluate the performance of each parameter point from the scan based on its largest individual pull. In short, each parameter point comes with a set of corresponding theoretical predictions for each observable, where the pull is defined as the difference between the experimentally measured value and the theoretical prediction, divided by the one sigma deviation of the measurement. As such, the pull is given in units of sigma, with each observable having its separate pull. When evaluating the performance of a parameter point, we use only the (magnitude of the) pull of the observable that deviates the most from its measured value, i.e. the largest individual pull.

A. Scalar potential scan
With the couplings of the scalar potential defined as in Ref. [19], we allow for the dimensionless, quartic couplings to vary over the range λ ∈ [10 −10 , 1] and the dimensionful, trilinear parameter to range over a 1 ∈ [−5, 0] TeV. The fixed input for each minimization, on the other hand, are the scalar VEVs, β and v S . These parameters are picked to lie on a 100 × 100 grid, with β-values chosen such that tan β ∈ [10 −3 , 10 3 ] in log-scale and with v S ∈ [1, 10 4 ] TeV.
For each of the 10 000 possible values for (β, v S ), we then find values for the scalar potential parameters that optimize: • The tree-level tadpole equations being fulfilled; • The lightest massive neutral scalar having the mass of the SM Higgs boson, m H = 125.10 ± 0.14 GeV [34]; • One of the Goldstone bosons and the Higgs field being aligned with the SM limit; • The scalar contributions 3 to the oblique S, T and 3 The remaining contributions to the oblique parameters are incorporated at a later stage (in the full scan).
By then removing any parameter points with a largest individual pull above 2σ, 97% of the initial 10 000 points survive. That is, the output of this part of the scan is 9 700 points in the (β, v S ) space.
Note that we allow for a (small) misalignment, as the alignment condition is a part of the optimization rather than being implemented as a strict constraint. The tadpole equations are also a part of the optimization procedure, as our code is constructed to handle any number of Higgs doublets and any number of complex scalar singlets, where general analytical solutions are not available. In simpler systems, where the analytical solutions are known, we carefully monitor any deviation on this part.

B. Leptonic scan
The leptonic sector contains the dimensionless Yukawa couplings |y| ∈ [10 −10 , 1] with arg(y) ∈ [10 −10 , 6.28], the dimensionful Majorana parameters in A ∈ [10, 10 5 ] TeV and the dimensionless Majorana parameters in B, C ∈ [10 −3 , 1]. For each (β, v S ) grid point that survives the minimization of the scalar sector, we then scan, for both inverted ordering (IO) and normal ordering (NO), over these ranges to fit: • The running charged lepton masses in [35]; • The two squared mass differences for neutrinos [34], Here, any points with an individual pull above 4σ are disregarded, which further reduces the (β, v S ) parameter space. The percentage of (β, v S ) values that survive this cut, in addition to the number of β values this corresponds to, are shown in Tab. I. For example, for model A1 with IO, 63% of the (β, v S ) values survive, i.e. 6 300, which in this specific case corresponded to 93% of the β values remaining, i.e. 93 values.
In the analysis carried out in Sec. VI, we see that IO outperforms NO for all models. 4 Hence, the remainder of the scan treats solely IO.  ) values out of the initial 10 000 after the leptonic scan (2nd column), percentage of surviving tan β values relative to the initial 100 before and after the quark scan (3rd column), and the number of input points and best-fit points after the full scan (4th column).

C. Quark scan
For the quark sector, only β enter as an input, while the magnitude and argument of the Yukawa couplings are allowed to range over the same values as in the leptonic sector. The scan procedure is then finding the set of Yukawa parameters that optimize the fit of: • The running quark masses in [35]; • The angles and CP phase of the CKM mixing matrix, parametrized in terms of the Wolfenstein parameters [34] λ = 0.22453 ± 0.00044, A = 0.836 ± 0.015, The percentages of surviving points, after using an individual pull of 2σ as the cut, are again shown in Tab. I. For model A1 with IO, this is 91 %, i.e. 91 out of the 93 β values survive the quark scan. As the two β values that were killed off can at most correspond to 100 v S values each, model A1 with IO has a minimum of 6 100 (β, v S ) values after the quark scan.
Unlike in the previous scans, we save all parameter points below the 2σ limit. There can hence be several points with the same β value, but with different values for the magnitude and argument of the Yukawa couplings. As a result, the number of (β, v S ) parameter points can here exceed the initial 10 000. For example, model A1 with IO has 78 215 such points, used as input for the full scan.

D. Full scan
In the full scan, we then combine the output parameters from all previous minimizations of the largest individual pulls and use them as fixed input values. The only free parameter left to adjust is hence the gauge coupling of U (1) , which we allow for to vary in the range g ∈ [5×10 −4 , 1]. This scan contains, on top of all sectors previously described, phenomenological constraints for: • Electroweak observables; • Meson sector observables; • Collider constraints.
Here, the electroweak observables include Z-pole pseudo observables, oblique parameters, off-pole crosssections, rare top decays, atomic parity violation, electric dipole moments and muon magnetic moments, while the meson observables involve mass splittings, kaon sector CP asymmetry, B-sector CP-violating observables, leptonic decay and radiative decay. For the collider constraints, we consider only the ones coming from direct searches of the Z boson. For more details, see Ref. [19].
Besides the observables considered in Ref. [19], we include two additional lepton flavor violating (LFV) observables, namely two kinds of charged lepton decay -→ γ and → 3 . To evaluate the new physics (NP) contribution to → γ, we begin with defining the effective Hamiltonian Leading order NP contribution to → γ. Note that the both dashed and wiggled line is used to indicate that the propagator can be either a Z boson or a neutral NP scalar.
Leading order NP contribution to j → i l k .
with the operators, for on-shell matching, defined as From matching this (at the NP scale) to the leading order NP contributions in Fig. 1, the Wilson coefficients are given by where F 2 and G 2 are the so-called Pauli-and EDM form factors, calculated with (and defined as in) Package X [36], for each parameter point. To verify the result, we compared its analytic form in the limit of massless initialand final states with the formulae presented in Ref. [37]. Note that the evaluation of the form factors in their exact form, i.e. with no such limit taken, requires high precision for numerical stability, and also that the only contributions to F 2 and G 2 in Fig. 1 come from diagrams where the detached photon is attached to the leptonic propagator.
For the 3-body lepton decay, j → i l k , the leading order NP contribution is instead a tree-level diagram, as shown in Fig. 2. Here, the effective Hamiltonian (in the massless final state approximation) is given by with an implicit summation over repeated indices and with X, Y = L, R. From matching the effective scattering amplitude to the amplitude shown in Fig. 2, we then obtain the Wilson coefficients with ∆ and Π defined as where S 0 is any neutral NP scalar and where all particles are defined to be in their mass basis. Note that, in Feynman-'t Hooft gauge, the Goldstone contribution is zero in the case of neglecting the masses of final state leptons.
The branching ratio is then given by with N s = 2 in the case of having two identical particles in the final state (e.g. τ − → µ − µ + µ − ), and N s = 1 otherwise, in agreement with Ref. [38]. Here, the overall factor of 1/2 is the phase space reduction in the case of having two identical particles in the final state, while the prefactors of N s comes from there being two possible contractions when X = Y and i = k. For pioneering work on this subject, see Ref. [39].
Note that there are also contributions coming from the SM fields that mix with either Z or S 0 . This contribution is incorporated into the Wilson coefficients at the scale where the SM fields are integrated out, and have the same form as in the equations presented above. However, as this contribution is suppressed in the case of a small mixing angle, it is often negligible.
The number of best-fit points after the full scan, i.e. number of points with the lowest individual pull for each minimization, is shown in Tab. I. At this stage, any points with an individual pull above 6σ are disregarded. The reason for allowing for a larger deviation here than in the previous scans, is to properly display the shape of the distribution when plotted in Sec. VI. Also, as there are large hadronic uncertainties in the meson sector, a model with e.g. a 4σ deviation should not necessarily be disregarded, provided that the largest individual pull comes from an observable in the meson sector and not the EW sector.
Note also that the range for g is divided into four equal parts, with each parameter point scanned for all ranges. This is why the amount of best-fit points can exceed the number of input points for some of the models in Tab. I, but can never go higher than four times its value. As a rough measure of the overall performance of a model, we can hence compare the number best-fit points with four times the number of input points. For example, model A1 with IO has a ratio of 40850/(4 · 78215) ∼ 0.13, to be compared with e.g. ∼ 0.35 for A3 and ∼ 0.18 for A4.
The four models with the largest ratio are A3, A6, B1 and B2, while no parameter points survive the full scan for models A2, B3, B4, D1 and D2. Note however that the ratio has a slight bias due to the various models not having exactly the same amount of statistics (the number of parameter points having survived the earlier steps in the scan of course varies in between models).

VI. COMPARING THE MODELS
In this section, distribution of the largest individual pulls are depicted in the form of box plots. In a box plot, the distribution is ordered by magnitude and then split into four equally sized parts, commonly referred to as quartiles. The box itself extends from the first to the third quartile, while the so-called whiskers extend out to the value with the largest magnitude within 1.5 times the length of the box in either direction. If there are any points small or large enough to be outside the extent of the whiskers, these are classified as outliers and shown as isolated black dots. The middle line corresponds to the median of the set.
Note that the largest individual pull is not necessarily a fair measure. For example, if we have a parameter point with only one deviating observable, it would be considered to perform equally to that of a parameter point where e.g. all observables deviate with that same amount. However, it is still a useful first indicator to whether a model is promising or not.

A. The leptonic sector
With each of the 16 models in Sec. IV coming in two varieties -one with NO and one with IO -there are a total of 32 models for us to compare. The result from the second stage of the scan, i.e. from comparing the six lepton textures for IO and NO, is shown in Fig. 3. Here, we see that all textures tend to perform better for IO, with the lowest individual pull ranging down to values below 10 −2 σ for all textures but E5, rather than down to around 1.5σ, as is the case for NO. In fact, for lepton texture E5, there are no valid parameter points at all for NO below the cut-off.
The source of this behavior can be identified by studying the pulls of individual observables, as shown for lepton texture E6 in Fig. 4. Here, we see that the limiting observable for NO is the CP phase, with a pull that never reaches below 2σ. The same type of behavior is exercised by all of the lepton textures, which is why the quark scan and the full scan was carried out solely for lepton textures with IO, reducing the total number of models from 32 to 16.
FIG. 4. Distribution of pulls for model E6, with inverted ordering (left) and normal ordering (right).

B. The full scan
The distributions of best-fit points after the full scan are shown in Fig. 5. Note that five models are absent, namely model A2, B3, B4, D1 and D2, for which there were no parameter points with an individual pull below 7σ. Overall, there are five models (A1, A3, A4, B1 and B2) with parameter points for which the largest deviation of any observable is below 3σ, and neither of these points are outliers of the distribution. Among these, model A3, B1 and B2 have the largest ratio of best-fit points, as shown in Tab. I.

VII. CONCLUSIONS
Since the birth of 2HDMs, naturally flavor conserving implementations have been dominating the market. While they deserve the attention, their dominance has resulted in a sizable research gap. To start filling this gap, we previously classified all viable 2HDM-U(1) extensions with quarks and charged leptons in Ref. [19], and, in the present work, extended the analysis to include the neutrino sector in the instance of a type-I seesaw mechanism.
The three most promising models, based on the percentage of surviving best-fit points and the distribution of the largest individual pulls, are models A3, B1 and B2, all with an inverted-hierarchical neutrino mass spectrum.
While a scan is not a complete exploration of the parameter space, we find the results presented in this work promising. In particular, with best fit points being in the main bulk of the distribution, and not in some highly finetuned region of parameter space (they are not outliers), it should be relatively easy for future collaborations to find similar minima. As such, we hope that this work, and in particular model A3, B1 and B2, can prove useful for future studies.

Model
Xq