Axion-Neutrino Interplay in a Gauged Two-Higgs-Doublet Model

We propose a gauged two-Higgs-doublet model (2HDM) featuring an anomalous Peccei-Quinn symmetry, $U(1)_{PQ}$. Dangerous tree-level flavour-changing neutral currents, common in 2HDMs, are forbidden by the extra gauge symmetry, $U(1)_X$. In our construction, the solutions to the important issues of neutrino masses, dark matter and the strong CP problem are interrelated. Neutrino masses are generated via a Dirac seesaw mechanism and are suppressed by the ratio of the $U(1)_X$ and the $U(1)_{PQ}$ breaking scales. Naturally small neutrino masses suggest that the breaking of $U(1)_X$ occurs at a relatively low scale, which may lead to observable signals in near-future experiments. Interestingly, spontaneous symmetry breaking does not lead to mixing between the $U(1)_X$ gauge boson, $Z^\prime$, and the standard $Z$. For the expected large values of the $U(1)_{PQ}$ scale, the associated axion becomes"invisible", with DFSZ-like couplings, and may account for the observed abundance of cold dark matter. Moreover, a viable parameter space region, which falls within the expected sensitivities of forthcoming axion searches, is identified. We also observe that the flavour-violating process of kaon decaying into pion plus axion, $K^+ \to \pi^+ a$, is further suppressed by the $U(1)_X$ scale, providing a rather weak lower bound for the axion decay constant $f_a$.


I. INTRODUCTION
The origin of small neutrino masses and the nature of dark matter (DM) are two of the most pressing issues with no answers within the Standard Model (SM) of particle physics. Nevertheless, there exist plenty of other open questions suggesting the need of physics beyond the SM, e.g. the non-observation of a CP-violating phase in the strong interaction sector.
The observation of neutrino oscillations [1][2][3] has led to an understanding that, contrary to the SM picture, neutrinos are massive, albeit extremely light. A plethora of new physics proposals has been put forward to explain the smallness of neutrino masses: from the seesaw mechanism and its various realisations, see e.g. [4] -relying on new physics at very large scales -to radiative mechanisms [5][6][7][8] -taking place at much lower scales, possibly within experimental reach. On the experimental side, many efforts have been helping us to determine not only neutrino masses per se but also other intrinsically related properties, such as neutrino mass ordering and absolute scale, CP phase and whether neutrinos are their own anti-particles [9].
Another major drawback of the SM is the absence of a suitable candidate to account for the observed dark matter relic abundance [10], whose evidence arises from many sources [11]: from studies of galaxy rotation curves to cosmic microwave background data. Among the most appealing DM candidates, there are the weakly interacting massive particles or WIMPs, which, despite various experimental searches, have not yet been observed [12]. On the other hand, axions -originally proposed as a key ingredient of the Peccei-Quinn (PQ) solution to the strong CP problem [13][14][15] (for reviews, see [16,17]) -define another well-known class of DM candidates [18][19][20], having the advantage of being capable to evade strong constraints coming from WIMP searches.
The exciting possibility of linking neutrino masses to dark matter and the strong CP problem via axions has been investigated in different scenarios. For Majorana neutrinos, the implementation of seesaw mechanisms, where the large seesaw scale is identified with the PQ scale -already considered many decades ago [21,22] -has been explored in various proposals more recently, see e.g. [23][24][25][26][27][28].
The latter case has been attracting more attention over the last few years since experiments, such as searches for neutrinoless double beta decays [34], have not so far found any evidence for lepton number violation, which could confirm the Majorana nature of neutrinos.
In this work, we propose a model in which the issues of neutrino masses, dark matter and strong CP problem are addressed simultaneously. The SM group is enlarged by an extra gauge symmetry, U (1) X , as well as the global Peccei-Quinn symmetry, U (1) P Q . The model contains two Higgs doublets, as in two-Higgs-doublet models (2HDMs) [35], plus two singlets. As a result of the U (1) X charge distribution, the Higgs doublets couple to different fermions, preventing the emergence of dangerous flavour-changing neutral currents (FCNCs) at tree level. In the fermion sector, a minimal field content, including extra quarks and neutral leptons, is added to ensure gauge anomaly cancellation as well as a consistent generation of neutrino masses. Different constructions of 2HDMs with a U (1) X symmetry have already been proposed to explain the suppression of FCNCs [36], along with the implementation of WIMP [37] and axion [38] dark matter candidates, seesaw mechanism for the neutrino masses [39][40][41] and other phenomenological issues [42].
Neutrino masses are generated via a Dirac seesaw mechanism, coming with the suppression factor v ϕ /v σ , where v ϕ and v σ are the U (1) X and U (1) P Q scales, respectively. Naturally small neutrino masses suggest that the breaking of U (1) X occurs at a relatively low scale, v ϕ v σ , which may lead to observable signals in near-future experiments. Interestingly, spontaneous symmetry breaking does not lead to mixing between the U (1) X gauge boson, Z , and the SM Z. For the expected large values of the U (1) P Q scale, v σ , the associated axion becomes "invisible", with Dine-Fischler-Srednicki-Zhitnitsky (DFSZ)-like couplings [43,44], and may account for the observed abundance of cold dark matter. Moreover, a viable parameter space region, which falls within the expected sensitivities of forthcoming axion searches, is identified. Therefore, the alluring axionneutrino interplay renders the model theoretically consistent and phenomenologically rich.
The remaining of the paper is organised as follows. In Sec. II, we discuss the model building rationale and present the field content and symmetry properties. In Sec. III, the scalar spectrum is derived, and the orthonormalisation of the Goldstone bosons is thoroughly discussed. Next, in Sec. IV, we obtain the gauge spectrum, augmented by the new boson Z , and elucidate the non-mixing property between Z and Z . The fermion sector is explored in Section V, where the mixing patterns between the extra and the standard fermions are obtained via the diagonalisation of their mass matrices. Naturally small neutrino masses are shown to be generated via a Dirac seesaw mechanism. We turn our attention to the axion physics in Sec. VI, where we show the main axion properties, derive the model-dependent axion couplings to photons and fermions as well as investigate relevant phenomenological consequences. Our final remarks are made in Sec. VII.

II. MODEL BUILDING
We start by considering a 2HDM with a U (1) X gauge symmetry under which the Higgs doublets, Φ u and Φ d , carry different charges, forbidding the appearance of dangerous Higgs-mediated FCNCs.
In addition to the Higgs doublets, we introduce a SU (2) L Higgs singlet, ϕ, also charged under the new local group. We assume that ϕ acquires a vacuum expectation value (vev) above the electroweak scale, spontaneously breaking U (1) X and generating a mass to the associated gauge boson Z . Taking advantage of the presence of two Higgs doublets, fundamental ingredients of DFSZ-type axion models, an anomalous Peccei-Quinn (PQ) symmetry U (1) P Q is implemented.
This is made possible with the introduction of a second singlet, σ, whose large vev breaks U (1) P Q , triggering the PQ mechanism that solves the strong CP problem. The (pseudo-) Goldstone boson of the U (1) P Q spontaneous breaking is identified with the axion field, and it can play the role of cold dark matter.
In order for the Peccei-Quinn symmetry to be realisedà la DFSZ, the Higgs doublets must also be charged under U (1) P Q , and each of them, namely Φ u and Φ d , should couple to either the righthanded up-type quarks, u aR , or down-type quarks, d aR , respectively. However, if only standard quarks are present, U (1) X anomaly cancellation -in particular, for the [SU (3) C ] 2 ×U (1) X anomaly -is achieved once the scalar doublets are identically charged under U (1) X : to have X Φu = X Φ d so that U (1) X is responsible for the absence of Higgs-mediated FCNCs, we need to extend the quark sector of our model in such a way that all U (1) X anomalies vanish.
In the quest for minimal solutions, we add n pairs of quarks (vector-like under the SM group), k nL,R , carrying the same electric charge q k , and try to find minimal sets of (n, q k ) for which all U (1) X anomalies are cancelled. This is obviously only possible if the new quarks are chirally charged under U (1) X , and, for simplicity, we assume that they get their U (1) X charges, as well as masses, via tree-level couplings to ϕ. We can divide our search into two major cases depending on whether the right-handed charged leptons, e aR , couple to Φ d , as in the type-I DFSZ model or type-II 2HDM, or Φ u , as in the type-II DFSZ model or flipped (type-Y) 2HDM. In the former case, one of the simplest solutions is (n, q k ) = (3, 2/3), whereas for the latter case, we find (n, q k ) = (3, −1/3).
In this work, we focus on the second case, which requires n = 3 pairs of extra quarks carrying the same electric charge as the down-type quarks: At last, we include three right-handed neutrinos, ν aR , and three pairs of neutral leptons n aL,R , which are vector-like under the gauge symmetries. The presence of such fields allows for a consistent generation of small neutrino masses via a Dirac seesaw mechanism, taking place via an interplay among all scales in the model.
In Table I, we present the fermion and scalar contents together with their symmetry transformations. In the U (1) global column, we have five independent charges which can be linked to the five Abelian symmetries in the model: the five global symmetries of the model, including U (1) P Q for which q σ ≡ P Q σ = 0. Amongst them, three satisfy anomaly-free conditions, displayed in the U (1) af ree column, including U (1) X for which l ϕ ≡ X ϕ = 0. and U (1) L , where the last two are the baryon and lepton number symmetries. For instance, let us choose the generator basis to be (q n L , q Q L , q Φu , q ϕ , q σ ). In this case, the symmetries U (1) L and U (1) B are generated by (1, 0, 0, 0, 0) and (0, 1/3, 0, 0, 0), respectively. For U (1) Y , the generator can be identified as (0, 1/6, 1/2, 0, 0), which is clearly linearly independent, but not necessarily orthogonal, with respect to the previous two generators. The generators of the last two symmetries, U (1) X and U (1) P Q , are also linearly independent among themselves and with respect to the other three. By construction, the generator of U (1) X has a non-zero fourth and a zero fifth entry: (X n L , X Q L , X Φu , X ϕ = 0, X σ = 0), whereas U (1) P Q is the only symmetry for which the last entry must be different from zero: (P Q n L , P Q Q L , P Q Φu , P Q ϕ , P Q σ = 0). The exact charges that define these generators will be properly derived in the next section once the scalar spectrum is obtained, in particular, when the orthogonal Goldstone bosons are identified. This procedure allows for the unambiguous identification of the physical charges, preventing any misleading choice [45]. Finally, the U (1) af ree column represents a subgroup of U (1) global , containing only anomaly-free solutions.
To find U (1) af ree , we impose that all the coefficients arising from the anomalies below must vanish For instance, the vanishing of the anomaly coefficient in I is achieved for q σ ≡ l σ = 0. As for the coefficient II, in addition to the previous constraint, its vanishing requires that q n L ≡ l n L = −3l Q L + l Φu + l ϕ . Once these two constraints are imposed, all anomaly coefficients vanish identically, as shown in Appendix A. Consequently, after the imposition of two constraints, the number of independent charges goes from five, in the U (1) global column, to only three, in the U (1) af ree column. Such charges can be grouped in the basis (l Q L , l Φu , l ϕ ) and be identified as the generators of U (1) B−L : (1/3, 0, 0), U (1) Y : (1/6, 1/2, 0) and U (1) X : (X Q L , X Φu , X ϕ = 0).
Although the symmetries U (1) Y , U (1) X and U (1) P Q are all broken spontaneously, the U (1) B−L symmetry will remain intact, ensuring the Dirac nature of neutrinos, whose masses are generated via a (Dirac) seesaw mechanism.

III. SCALAR SECTOR
The scalar potential can be written as In the limit λ 4 → 0, the scalar potential has four independent U (1) global symmetries related to the phase redefinitions allowed for each scalar field. When the λ 4 term is introduced, one of the four possible linear combinations of these symmetries is explicitly broken so that V is left invariant under only three Abelian groups. Two of them can be identified with the gauged U (1) Y and U (1) X symmetries, while the remaining one is the global U (1) P Q symmetry. As discussed in Sec. II, the U (1) P Q (U (1) X ) charges can be obtained from the U (1) global (U (1) af ree ) column in Table I when taking q σ ≡ P Q σ = 0 (l ϕ ≡ X ϕ = 0).
In order to derive the scalar spectrum, we decompose the scalar doublets as with v 2 u + v 2 d ≡ v = 246 GeV, whilst for the singlets, we have Once all scalars acquire vevs, the following spontaneous symmetry breaking pattern takes place where G SM stands for the SM gauge group and the global symmetries are shown in parentheses.
Notice that, as discussed in Sec. II, no scalar field is charged under the accidental U (1) B and U (1) L symmetries so that they remain fully conserved. The breaking of five generators leads to four would-be Goldstone bosons, absorbed by the gauge sector, plus a pseudo-Goldstone boson, the axion, as we derive in what follows.
Substituting Eqs. (3) and (4) into Eq. (2), we extract the following tadpole equations To find the physical spectrum, we solve the equations above for the dimensionful parameters µ u , µ d , µ ϕ and µ σ , and plug them back into the potential.
The scalar spectrum contains two charged fields, which are defined in terms of (φ ± u , φ ± d ) as The first field, φ ± , is a physical charged scalar, whose mass can be around the electroweak scale, as in 2HDMs, as long as λ 4 remains very small. The smallness of such a parameter is naturally protected since in its absence the potential exhibits an enhanced symmetry. The second scalar, G ± , which remains massless, is the Goldstone boson absorbed by the gauge sector, making the SM vector boson W ± massive.
For the neutral fields, we divide them into the CP-even and the CP-odd components. Starting with the CP-even scalars, in the basis (s u , s d , s ϕ , s σ ), we can write the squared mass matrix below The matrix in Eq. (8) contains the three energy scales present in the model, and it is expected that two scalars will be heavy with masses proportional to the vevs v ϕ and v σ . It is a typical feature of the axion models that the mass matrix of the CP-even scalars contains hierarchical vacuum expectation values. To have a Higgs boson consistent with the observed one, an adjustment of the parameters is required. We will not develop it further once this is not the focus of the present work.
A. CP-odd sector: identifying Goldstone bosons and abelian charges As a result of the polar parametrisation in Eqs. (3) and (4), the terms in the potential involving only the CP-odd states can be succinctly written as Upon expanding the cosine function, we find that only one state becomes massive at this point.
The massive state is proportional to the argument of the cosine function, which, when normalised, with a squared mass given by The pseudoscalar field A becomes massless in the limit λ 4 → 0. In fact, as mentioned below Eq.
implies the existence of an extra global symmetry in the scalar potential whose spontaneous breaking would identify the field A as the associated Goldstone boson. Thus, under the assumption that λ 4 can be naturally small (λ 4 1), since its vanishing increases the symmetries of the system, m A could also be around the electroweak scale, for example. Moreover, according to the vev hierarchy, A couples predominantly to the Standard Model fields once its components are mainly along the a u and a d field space directions.
As for the remaining fields, they are Goldstone bosons associated with the spontaneous breaking of three abelian symmetries: U (1) Y , U (1) X and U (1) P Q , defining a degenerate 3-d space. In order to identify the three linearly independent CP-odd scalars, it is instructive to write down the conserved current associated with each U (1) symmetry in the model along the CP-odd scalars.
As usual, we assume that under a given global U (1) c symmetry, a scalar field φ transforms as and ω c is the infinitesimal continuous parameter of U (1) c . Noether's theorem tells us that the presence of a U (1) c symmetry -in our case c = Y, X, P Q -implies the conservation of the following current where the ellipsis corresponds to the contributions from all non-scalar fields charged under the symmetry. Using the polar decomposition for the scalars, as in Eq. (3) and (4), we find the conserved current along the CP-odd scalars to be where, in the last step, we have defined the linear combination with The field G c is precisely the massless field associated with the spontaneously broken U (1) c generator as predicted by Goldstone's theorem ( 0|J c µ |G c = ip µ G c ). We are now well equipped to determine the Goldstone bosons by applying the expression in Eq. (14) to our model. Moreover, by imposing the physical condition of orthogonality among the CP-odd states, we are able to fix the U(1) charges of the scalars in terms of the vevs, ensuring that the charges in Table I are unambiguously chosen.
The first (would-be) Goldstone boson, G Z ≡ G Y , comes from the breaking of U (1) Y and is absorbed by the massive vector boson Z via the Higgs mechanism. As only the SU (2) L doublets carry hypercharge, we can use Eq. (14) to obtain, as expected, Notice that G Z is automatically orthogonal to A, as it should be. Had we not known beforehand the hypercharges of the doublets, we could also have identified Eq. (15) from its orthogonality to A, which would in turn provide us with the relation Y Φu = Y Φ d .
A second would-be Goldstone boson, G Z ≡ G X , emerges when the gauged U (1) X symmetry is spontaneously broken. G Z can be properly identified by noticing that it has components along the scalars charged under U (1) X , i.e. Φ u,d and ϕ (see Table I), as well as it must be orthogonal to A and G Z , giving By comparing Eq. (16) and (14), we find the unambiguous U (1) X charge relations: Without loss of generality, we normalise the U (1) X charges by setting: X ϕ = 1.
We would like to emphasise that once the orthogonality among the Goldstone bosons is imposed, the U (1) X charges of the scalars in Table I cannot be chosen freely. This feature leads to a very distinctive implication to the extended gauge sector phenomenology: no tree-level mass mixing between the SM and U (1) X gauge bosons is generated, as discussed in the next section.
Finally, we can proceed to the last CP-odd state, the (pseudo-)Goldstone of the anomalous U (1) P Q symmetry, the axion a ≡ G P Q , which can be easily obtained by requiring it to be orthogonal to A, G Z and G Z : Once again, with the aid of Eq. (14), we identify the P Q charges of the scalars in terms of the vevs: and, as a normalisation condition, we can adopt P Q σ = 1.
Alternatively, the axion field above can be also identified by adding an explicit U (1) P Q -breaking term of the form κσ n + h.c. to the potential. Since σ is not charged under the other Abelian symmetries, see Table I, U (1) Y and U (1) X remain conserved. With the introduction of the new term, the CP-odd spectrum would contain two massive fields, one of which gets a mass proportional to κ. Then, in the limit κ → 0 -i.e. recovering the U (1) P Q -invariant potential in Eq. (2) -the κ-dependent mass goes to zero, so that the associated state can be identified with the Goldstone boson emerging from the spontaneous breaking of U (1) P Q . Such a field is exactly the axion given by Eq. (18). We will return to the axion field in Sec. VI to study its main properties and phenomenological features.
The relevant terms giving rise to gauge boson masses are where the covariant derivatives are given by Notice that the kinetic term for σ has being omitted since σ carries no local charge, i.e. D µ σ = ∂ µ σ.
When the scalar fields acquire a vev, the gauge bosons become massive via the Higgs mechanism.
The charged gauge boson, , whose associated would-be Goldstone boson, G ± , is defined in Eq. (7), gets the mass m W ± = gv/2. The neutral gauge bosons, in the basis , share the following squared mass matrix with the dimensionless parameters A and B given by where to obtain the right-hand sides of the equations above we have used the charges in Eq. (17) and the normalisation condition X ϕ = 1, which follow from the imposition of orthogonality among the Goldstone bosons. Due to the vanishing of A, the mass matrix in Eq. (22) becomes block diagonal. The upper 2 × 2 block mixes W µ 3 and B µ Y and is precisely what one gets in the SM, so that its diagonalisation generates the massless photon field, A µ , and the massive Z µ . On the other hand, the U (1) X field, Z µ ≡ B Xµ , remains unmixed 1 , and its mass is The interesting observation that Z remains unmixed is not exclusive of our construction. In fact, we expect this to happen in other U (1) X gauge extensions once one imposes that all Goldstone bosons must be orthogonal. This, however, shall be explored elsewhere.
From now on, we assume that Z is a heavy vector boson with m Z /g X v ϕ = 10 TeV.
With this benchmark choice, our model's predictions evade the current collider constraints on m Z , obtained from the analysis of dilepton final states at the LHC [47,48]. Cosmological constraints related to the effective number of extra relativistic species ∆N ef f ≤ 0.285 [10] are also relevant when selecting this benchmark. The cosmological constraint on ∆N ef f can be translated into a lower limit on m Z /g X (v ϕ ) since the light right-handed neutrinos, ν R , may thermalise with the SM fields in the early universe, via Z -mediated interactions, and then contribute to ∆N ef f . Whilst a detailed calculation of ∆N ef f is beyond the scope of the present work, we do not expect that our model's contribution to ∆N ef f will vary greatly with respect to that in Ref. [49] whose analysed model, a gauged U (1) B−L construction, shares important features with ours. Thus, the results in Ref. [49] have also been taken into account when selecting the benchmark v ϕ = 10 TeV.

V. FERMION SECTOR
We now turn our attention to the fermion sector, starting with the Yukawa interactions. The field content and its symmetry transformations, as shown in Table I, allow us to write the following renormalisable Yukawa Lagrangian: + y e ab L aL Φ u e bR + y n ab L aL Φ d n bR + y α ab ϕ n aL ν bR + y β ab σ n aL n bR + h.c. .
The structure of the Yukawa Lagrangian is similar to that of the so-called flipped or type-Y 2HDM in that the e R and u R couple to the same Higgs doublet: Φ u . Once the scalar fields acquire vevs, according to Eqs. (3) and (4), all fermions become massive, as detailed in the next subsections.

A. Charged fermion spectrum
We start by noticing that, similar to the flipped 2HDM, both the charged leptons and the uptype quarks get masses proportional to Φ u = v u / √ 2. The charged lepton masses can be obtained a 3×3 matrix, which can be diagonalised by performing the bi-unitary transformation: (U e L ) † M e U e R = diag(m e , m µ , m τ ). Likewise, the 3 × 3 up-type quark mass matrix is given by and its diagonalisation follows from the bi-unitary transformation The remaining quarks, d aL,R and k aL,R , when put together in the basis D L,R = (d, k) L,R , share the following 6 × 6 mass matrix Here, we adopt the ansatz in Refs. [51,52], which allows us to approximate the unitary matrices U D L,R as In this section, we show how neutrinos get naturally small masses via a type-I Dirac seesaw mechanism, which receives contributions from all the three energy scales present in the model.
The neutral lepton masses come from the last three terms in Eq. (25). When writing N L,R = (ν, n) L,R as the basis, we have the following 6 × 6 mass matrix with each element representing a 3 × 3 block.
The mass terms above are strictly of a Dirac type since U (1) B−L is conserved. This can be easily understood by noticing that, as discussed in Sec. II, the field charges under U (1) B−L are obtaining by taking l Q L = 1/3 and l Φu = l ϕ = 0 in the last column of Table I. Therefore, because no scalar field is charged under U (1) B−L , this symmetry is not broken spontaneously, implying that neutrinos, as well as all the other fermions, are necessarily Dirac fermions.
The texture of the Dirac mass matrix M N and the assumed vev hierarchy indicate that a type-I seesaw mechanism is in place. The diagonalisation of M N is achieved by a bi-unitary transformation (U N L ) † M N (U N R ) = diag(m ν 1 , m ν 2 , m ν 3 , m n 1 , m n 2 , m n 3 ), which can be divided into two steps by writing U N L,R = R N L,R V N L,R [51,52], given in Appendix B, similar to what has been done with the U D L,R matrices in Eq. (29). When the first transformation (with R N L,R ) is performed, we obtain the seesaw-suppressed mass matrix for the active neutrinos whose diagram is shown in Fig. 1 For large v σ , as in our benchmark above, the associated invisible axion can also account for the observed dark matter relic density. Therefore, the origin of small neutrino masses, a solution for the strong CP problem and the nature of dark matter go hand in hand in our construction. To illustrate the viability of our model, in Sec. VI (see Fig. 2), we identify a parameter space region within which the above-mentioned issues can be solved and that can be probed by forthcoming experimental searches. Furthermore, small neutrino masses, being proportional to v ϕ /v σ , also rely on the existence of a moderate U (1) X scale (v ϕ = 10 TeV). Therefore, TeV-scale U (1) X signatures, mediated by e.g. the extra gauge boson, Z , may also be within the reach of current or near-future experiments, such as the high-luminosity LHC.

C. Fermion couplings to vector bosons
In the previous sections, we have shown that our construction extends the SM field content not only by adding an extra gauge boson but also extra neutral leptons (n L,R , ν R ) and down-type quarks (k L,R ) which mix with their SM siblings. In what follows, we show that, due to these extra ingredients, fermions and vector bosons couple in a non-standard way.
The kinetic Lagrangian for fermions, giving rise to the fermion-vector boson interactions, can be, as usual, represented by where F spans through all the fermion fields in Table I   Since the extra neutral leptons and down-type quarks mix with the SM fields, flavour-changing neutral currents (FCNCs), mediated by both Z and Z , appear and are governed by the factors where U F L,R are the 6 × 6 unitary matrices that diagonalise the generalised down-type quark and neutral lepton mass matrices in Eq. (28) and (31), respectively, and I 3 (0 3 ) is the 3 × 3 identity (zero) matrix. The full vector and axial couplings are presented in Table II, in which we used the U (1) X charge definitions in Eq. (17), with the normalisation X ϕ = 1, as well as X Q L = 1/3. Table II and Eqs. (29) and (34), we find that the most relevant source of FCNC is given by the following Lagrangian involving the (known) down-type quarks and the Z boson Table II GeV. Therefore, the flavour-violating contributions come with an estimated suppression factor of at least 10 −6 , which is well below the current limit (< 10 −3 ) obtained from processes such as [54,55]. For the neutral leptons, the flavour-violating contributions, suppressed by v σ = 10 12 GeV, are much smaller and can be also safely neglected.

For instance, from the coefficients in
When it comes to the charged vector bosons, W ± , we also expect corrections to the SM contri- where where F varies through all the fermion mass bases and eQ F represents the electric charge of F .

VI. AXION PHYSICS
In this section, we bring our attention back to the axion field and discuss some of its relevant properties. To do so, it is convenient to rewrite the axion field, according to Eq. (14), as [56] a = 1 where the P Q charges, defined in Eq. (19) with P Q σ = 1, as well as f P Q , reparametrised in terms of two angles and where β and θ are defined as Notice that the first angle follows the conventional definition of β in 2HDM scenarios [35], whilst θ is expected to be small due to the assumed vev hierarchy v σ v ϕ v. Thus, the axion field is predominantly projected along the CP-odd component of the singlet σ: a a σ , and f P Q v σ . It is also worth pointing out that the charges in Eq. (40) satisfy the relation which arises from the only non-hermitian term in the scalar potential, i.e. λ 4 (Φ † d Φ u )(σϕ). In order to solve the strong CP problem through the Peccei-Quinn mechanism, the U (1) P Q symmetry must yield a nonzero [SU (3) C ] 2 × U (1) P Q anomaly coefficient C ag . This leads to the effective interaction for the axion field with the gluon field strength G b µν : in which α s = g s /(4π), where g s is the strong interaction coupling constant, andG b,µν ≡ µνσρ G b σρ /2 is the dual field strength. As we shall see in Eq. (49), the model has indeed a nonvanishing anomaly coefficient C ag = 3P Q σ = 3. From this, we define the axion decay constant f a = f P Q /N DW as well as the domain wall number N DW = C ag = 3 for the model. The same domain wall number occurs in the DFSZ-type model with the non-Hermitian term Φ † u Φ d σ in the scalar potential, but contrasts with the DFSZ-type model having instead the non-Hermitian term Φ † u Φ d σ 2 , which leads to N DW = 6 (see [57], for example).
The axion mass arises from nonperturbative QCD effects [14]. The leading order axion mass is given by in which m u (m d ) is the up-quark (down-quark) mass, m π the pion mass, f π 93 MeV and f a the pion and the axion decay constants [14]. Corrections to this formula, including higher orders in chiral perturbation theory and through lattice simulations, were obtained in [58][59][60]. Taking into account the benchmark v σ = 10 12 GeV for the U (1) P Q symmetry scale, defined in Eq. (41), we have that f a v σ /N DW , which leads to the axion mass m a 17 µeV.
Despite being very light, axions can account for a part of or even the total cold dark matter in the Universe, with their production realised through the vacuum re-alignment mechanism [18][19][20].
If the U (1) P Q symmetry breaking happened before or during the inflationary period of the Universe and was not restored afterwards, it is estimated that the axion field gives rise to a contribution to the dark matter relic density given by [57,61,62] Ω a h 2 ≈ 0.12 F θ 2 i f a 9 × 10 11 GeV where θ i is the initial misalignment angle which assumes values in the interval [−π, π], and the factor F accounts for anharmonicities that may occur in the axion potential. For example, with the benchmark f a 3.3 × 10 11 GeV, the axion could comprise the totality of the observed cold dark matter, i.e. Ω CDM h 2 = 0.12 [10], if F θ 2 i ≈ 3.2. This is consistent with cosmological observations if the Hubble expansion rate during inflation satisfies the constraint H inf 10 7 GeV [63], which follows from the non-observation of isocurvature fluctuations in the cosmic microwave background arising from quantum fluctuations of the axion field.

A. Axion coupling to photons
The interaction between an axion and two photons can be expressed as where the coupling g aγ is with α being the fine-structure constant. The factor −1.95 is model independent and comes from the ratio between the up-and down-quark masses, showing up in the calculation due to the mixing between axions and pions [23]. Moreover, we have the model dependent contribution C aγ /C ag -oftentimes defined as E/N in the axion literature -where C aγ and C ag are, respectively, the This result, C aγ /C ag = 2/3, is precisely what one gets in the type-II DFSZ model. Therefore, although our model has extra quarks which contribute to both C aγ and C ag , the ratio C aγ /C ag remains the same as in the DFSZ case. ADMX [67][68][69][70][71], HB bound [72,73], HDM [74][75][76]), whereas the dashed lines indicate projected experimental sensitivities (ABRACADABRA [77,78], MADMAX [79], IAXO [80], ALPS II [81][82][83]). The shaded orange and orange+yellow bands identify the preferred regions for neutrino masses once the limits from the Planck Collaboration [10,84] and KATRIN experiment [85] are considered, respectively, with effective Yukawas in the range 2.9 × 10 −6 Y ν ef f 1.2 × 10 −3 -see text for details.
The shaded region spans to the right of the plot up until the effective Yukawa reaches its minimum value, (Y ν ef f ) min , and the corresponding neutrino mass reproduces the largest possible value for m ν 3 according to constraints coming from the Planck Collaboration (orange)j m ν j < 0.12 eV -and KATRIN (orange+yellow) -m ee < 1.1 eV.
Finally, we would like to point out that, within these shaded regions, the model's prediction, i.e. the solid green line, also falls within the projected sensitivities of axion dark matter searches (MADMAX and ADMX) -dashed lines in brown.

B. Axion couplings to fermions
The charged leptons and the up-type quarks are charged universally under U (1) P Q , and, as a result, their couplings to axions are necessarily flavour-conserving and axial and can be written as The axion couplings to the neutral leptons and the down-type quarks bear similarities. In addition to the SM fermions, both sectors contain new fields which transform differently under the U (1) P Q symmetry, inducing axion-mediated FCNCs [86,87]. Due to flavour-violating interactions, scalar and axial couplings appear, which can be generically written as with F = N , D , and the mass basis, F = F L + F R , is related to the flavour basis via F = U F L,R F . The scalar and axial coefficients for the neutral leptons are given by with X N L,R defined according to Eq. (34), and the m N n are the eigenvalues of the neutral lepton mass matrix in Eq. (30). As for the down-type quarks, we have where m D m are the masses of the down-type quarks, i.e. the eigenvalues of the mass matrix in Eq. (28); X D L,R leads to FCNCs and follows from the definition in Eq. (34).
C. Constraining f a with a flavour-violating process It is possible to set constraints on the range of the axion decay constant by confronting our model's predictions with experimental bounds on flavour-violating processes involving an axion, such as the decays of heavy mesons. The most stringent constraint comes from the process K + → π + a and is given in terms of its branching ratio: Br(K + → π + a) 7.3 × 10 −11 [88].
At tree level, this branching ratio is evaluated as [86,87] Br where Γ tot 5.3×10 −17 GeV is the total decay width of K + , m K = 493.677 MeV and m π = 139.57 MeV are the kaon and pion masses [61], respectively, and the model-dependent contribution can be written as The relevant flavour-violating terms can be obtained from Eq. (34) and the Appendix B, leading to 12 . As derived in Appendix B, the matrices B D represent the mixing between standard and new quarks and, as such, are proportional to suppression factors given in Table III. By singling out the suppression factors, we can rewrite the X D 12 terms as represent all the matrix element products. Finally, we can substitute these expressions back into Eq. (55) and when comparing it with the experimental limit, Br(K + → π + a) 7.3 × 10 −11 , we find that f a must satisfy GeV .
Considering the benchmark assumed in previous sections for the scales in the model, i.e. (v d , µ, v ϕ ) = (10 2 , 10 3 , 10 4 ) GeV, we find that (µ 2 /v 2 ϕ ) sin 2 θ 8.3 × 10 −7 , which leads to the weak constraint: GeV. Thus, we observe that, since the FCNC contributions arise from the mixing with non-SM heavy quarks, the constraints on f a coming from flavour-violating processes are weakened. As a result, in this construction, supernovae limits on the axion decay constant [57] turn out to be more stringent.

VII. CONCLUSIONS
We have proposed a gauged two-Higgs-doublet model featuring an axion. Dangerous tree-level FCNCs, common in 2HDMs, are forbidden by the extra U (1) X gauge symmetry. Our construction suggests that solutions for the important issues of the nature of dark matter, the origin of neutrino masses and the strong CP problem may arise from the axion-neutrino interplay.
The extended scalar sector counts with two singlets, σ and ϕ, in addition to the Higgs doublets, Φ d and Φ u . As presented in Sec. III, when all scalars acquire vevs, satisfying the hierarchy v σ v ϕ v = (v 2 d + v 2 u ) 1/2 , spontaneous symmetry breaking takes place giving rise to five Goldstone bosons -four of which are, in fact, would-be Goldstone bosons absorbed by the gauge sector, while the last one is identified with a pseudo-Goldstone boson, the axion. At low-energies, the physical spectrum contains, besides the usual 2HDM degrees of freedom, the ultralight axion field, a, and a TeV-scale CP-even field. Moreover, we have shown that the imposition of orthogonality amongst the Goldstone bosons fixes the physical values for the U (1) P Q and U (1) X charges. As a direct consequence of this procedure, the tree-level mass mixing between the SM and the extra gauge boson, Z , vanishes identically.
The Yukawa sector bears similarities with the flipped (or type-Y) 2HDM in that the right-handed charged leptons and up-type quarks couple to the same Higgs doublet Φ u , while the down-type quarks and the neutral leptons couple to Φ d . To ensure gauge anomaly cancellation and generate small neutrino masses, we have introduced extra quarks, k aL,R , chirally charged under U (1) X , righthanded neutrinos, ν aR , and extra neutral leptons, n aL,R . The extra quarks get masses proportional to v ϕ and mix with the SM down-type quarks, leading to FCNCs mediated not only by Z but also a, the axion field. For the neutral leptons, a Dirac seesaw mechanism generates small masses for the active neutrinos, see Fig. 1. Neutrino mass suppression is controlled by the ratio v ϕ /v σ , where we have taken v σ = 10 12 GeV as the large PQ scale, suggesting that the U (1) X symmetry is broken at a much lower scale, v ϕ . On the other hand, since Z and the quarks k a , get masses around v ϕ , this scale is constrained from below by experiments and cosmology. These features suggest that v ϕ lies within a phenomenologically rich region, and we have assumed v ϕ = 10 4 GeV as a benchmark.
Finally, we have studied the properties of the axion in Sec. VI. The axion has components along all the scalars but is mostly projected along the CP-odd component of the singlet σ with a decay constant: f a v σ /N DW = v σ /3. For the chosen benchmark of v σ = 10 12 GeV, which leads to naturally small neutrino masses, the associated axion may account for the totality of the observed cold dark matter, considering the pre-inflationary scenario for axion production. Despite the presence of extra fermions, the axion coupling to photons depends on the same anomaly coefficient ratio as in the type-II DFSZ model: C aγ /C ag = E/N = 2/3. Fig. 2 shows that the preferred region for neutrino masses and axion dark matter may be tested by forthcoming axion experiments looking for axion-photon interactions. Furthermore, we have derived the axion couplings to fermions and noticed that a rather weak lower bound for f a is obtained from flavour-violating processes, such as charges: (q n L , q Q L , q Φu , q ϕ , q σ ), given in Table I. Among the global symmetries, only those which lead to vanishing anomaly coefficients in Eq. (1) can be safely promoted to local. To find the anomaly-free solutions, we calculate the coefficients using the charges in the U (1) global column and set them to zero, i.e.
whereq = 3q Q L + q n L − q Φu − q ϕ ; the sum in I takes only quarks into account, II considers only fermion doublets, whereas for the remaining coefficients all fermions contribute. It is easy to see that the equations above are simultaneously satisfied when q σ =q = 0, reducing the number of free parameters from five to three. Finally, by renaming the three independent charges of this subset as (l Q L , l Φu , l ϕ ), we obtain the U (1) af ree column of Table I.

Appendix B: Diagonalisation
Our task is to diagonalise the fermion mass matrices present in this work. To this end, we will first see how to block diagonalise a Hermitian matrix. Then, we will block diagonalise a non-Hermitian matrix as a generalisation of the previous procedure. Lastly, we will show how all the fermion mass matrices in this work can be fully diaogonalised.

Block diagonalisation of Hermitian matrices
Consider a 6×6 Hermitian matrix M . If we want to block diagonalise this matrix, it is sufficient to find a unitary matrix R, such that where M i , i = 1, 2, are 3 × 3 matrices. The procedure outlined here will work whenever the matrix M can be written as where, A 1 , A 3 are Hermitian matrices, the entries of the matrices A i , i = 1, 2, 3 are dimensionless and of order unity, whereas m i represent mass scales that satisfy the hierarchy m 3 m 1 , m 2 . To achieve the block diagonalisation, we parameterise the unitary matrix R as where B is a general complex 3 × 3 matrix to be determined. The square roots should be seen as series expansions in B [51, 52] In turn, B can also be expanded as B = B 1 + B 2 + B 3 + ..., where B n is of order n in the expansion parameter = m 2 /m 3 1. An alternative parameterisation for the matrix R can be found in Ref. [53].

Block diagonalisation of non-Hermitian matrices
Now, consider a 6 × 6 non-Hermitian matrix M F . Our next task is to block diagonalise the down-type quark and the neutral lepton mass matrices. As neither of them is Hermitian, instead of the unitary transformation in Eq. (B1), a bi-unitary transformation is required to block diagonalise each of them, where, R F L,R are unitary matrices and F 1,2 are the components of the basis F = D, N . The problem is solved if we find the matrices R F L,R . One way of doing it is to break the bi-unitary transformation into two unitary transformations. Multiplying Eq. (B6) by its Hermitian conjugate gives Similarly, inverting the product order, we multiply the Hermitian conjugate by Eq. (B6) and find The matrices D F L,R are block diagonal since M F block also is. The matrices M F M F † and M F † M F are Hermitian, as a result, the Eqs. (B7) and (B8) represent unitary transformations, and we can make use of Eq. (B3) to find the unitary matrices R F L,R that block diagonalise M F , up to the desired order. Thus, we can write the unitary matrices that diagonalise M D and M N , in Eqs. (28) and (30), respectively, in terms of the contributions in Table III.
ϕ y d y µ † y k y k † −1 v d vσ y n y β † y β y β † −1 B F R µ vϕ y µ † y k y k † y k −1 vϕ vσ y α † y β y β † y β −1 Table III. Approximation for the 3 × 3 matrices B F L,R , with F = D, N , up to the first order term in the expansion. The orders of magnitude of B F L,R are given by the vev relations for Yukawa matrices of order unity.

Diagonalisation of the mass matrices
We are ready to diagonalise the mass matrices of all fermions present in this work. Consider a 6 × 6 block diagonal matrix, M F block , which is composed of two non-Hermitian 3 × 3 blocks, M F 1,2 . This is the structure for the down-type quarks and neutral leptons of the model after the block diagonalisation. To completely diagonalise M F block , we need another bi-unitary transformation, where M F is the diagonal matrix for the mass eigenstates in the basis F = D , N , V F L,R = diag V F 1 L,R , V F 2 L,R are 6×6 unitary matrices, and V that is, the combined effect of two diagonalisation steps is equivalent to a complete diagonalisation performed by the unitary matrices Therefore, we can summarise all the individual diagonalisation procedures as where F = e, N, u, D, and U F R,L has the same dimension as the corresponding mass matrix M F .