A class of invisible axion models with FCNCs at tree level

We build a class of invisible axion models with tree-level Flavor Changing Neutral Currents completely controlled by the fermion mixing matrices. The scalar sector of these models contains three-Higgs doublets and a complex scalar gauge singlet, with the same fermionic content than the Standard Model. A horizontal Peccei-Quinn symmetry provides a solution to the strong CP problem and predicts the existence of a very light and weakly coupled pseudo-Goldstone boson, the invisible axion or familon. A phenomenological analysis is performed taking into account familon searches in rare kaon and muon decays, astrophysical considerations and axion searches via axion-photon conversion. Drastic differences are found in the axion properties of different models due to the strong hierarchy of the CKM matrix, making some of the models considered much more constrained than others. We also obtain that a rich variety of these models avoid the domain wall problem. A possible mechanism to protect the solution to the strong CP problem against gravitational effects is also discussed.


Introduction
The recent discovery by the ATLAS [1] and CMS [2] collaborations of a Higgs-like particle with a mass around 125 GeV represents one of the greatest achievements of physics in the last decades and constitutes an indisputable success of the Standard Model (SM) of particle physics. This discovery reinforces the SM as the theory of electroweak and strong interactions, no significant deviations from this framework have been observed so far at LHC, in precision experiments at flavor factories nor in electroweak precision test at LEP. However, the SM presents several unanswered questions which might be a hint of physics beyond the SM. One of such open questions is the so-called strong CP problem [3].
The strong CP problem is tightly related to the U(1) A problem of Quantum Chromodynamics (QCD). In the limit of massless quarks the QCD Lagrangian shows a chiral U(1) A symmetry.
The fact that after chiral symmetry breaking its associated pseudo-Goldstone boson was not found experimentally proved that this symmetry should be broken or not realized in nature [4].
This led to an apparent contradiction between theory and experiment which was termed as the U(1) A problem. The solution to this issue came from the realization by t'Hooft that nonperturbative QCD effects explicitly break this symmetry [5]. However, with the resolution of where g s is the strong coupling constant and G µν a andG µν a are the QCD field-strength tensor and its dual tensor, respectively. This way the QCD vacuum angle, θ QCD , together with g s remain as the only free parameters of massless QCD. If along with QCD the electroweak (EW) sector is introduced, one should take into account that the quark masses are complex in general.
To get the Lagrangian in the physical basis a chiral U(1) A transformation should be performed.
As a result, the QCD vacuum angle in Eq. (1) is substituted in the full theory byθ defined as being M the quark mass matrix. Forθ = 0, Eq. (1) introduces a violation of P and T but not C and consequently a violation of CP. However, the present bound on neutron dipole moment, |d n | < 2.9 × 10 −26 e cm [6], set a stringent bound on this angle θ < ∼ 10 −11 [7]. The reason why this parameter, coming from the strong and the electroweak sectors, is so small is unknown and gives rise to the Strong CP problem.
An elegant solution to the Strong CP problem was given by Peccei and Quinn [8]. This solution, commonly referred as the Peccei-Quinn (PQ) mechanism, consists on the introduction of a global chiral U(1) PQ symmetry with mixed anomalies with QCD. This symmetry effectively replaces the CP-violating phase by a CP-odd field, the so-called axion, which correspond to the pseudo-Goldstone boson resulting from the spontaneous breaking of the PQ symmetry [9]. The implementation of the PQ mechanism requires the extension of the matter content of the SM.
In its original formulation, the scalar sector is enlarged to a two-Higgs-doublet model (2HDM) with the PQ charges implementation enforcing Natural Flavor Conservation (NFC) [10]. This way the severe experimental constraints from Flavor Changing Neutral Currents (FCNCs) [11] are avoided. In this model the axion has a mass of few hundreds keV and presents large couplings to matter [9]. This formulation was soon ruled out by experimental data.
In order to satisfy the experimental constraints one needs to decouple the PQ symmetry breaking and the EW scales. This is achieved by the introduction of a gauge singlet field that acquire a vacuum expectation value (vev) that breaks the PQ symmetry at a scale much higher than the EW scale. This results in invisible axion models where the mass and the couplings of the axion are suppressed by the vev of the scalar singlet and, therefore, are naturally small. In this class of models the axion possesses several interesting features. For instance, the invisible axion is a promising candidate for cold dark matter [12]. Additionally, the type I seesaw mechanism [13] can be naturally implemented in these models allowing for the possibility to explain the smallness of the active neutrino masses and providing a dynamical origin to the heavy seesaw scale [14].
Two models stand as benchmark invisible axion models: the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) [15] and the Kim-Shifman-Vainshtein-Zakharov (KSVZ) [16] models. In the KSVZ model one adds to the SM particle content a heavy color triplet and SU(2) L singlet vector-like quark and a complex scalar gauge singlet. The SM fields carry no PQ charge in the KSVZ model. On the other hand, in the DFSZ model one introduces an additional Higgs doublet and a complex scalar gauge singlet while the PQ symmetry enforces NFC just like in the original PQ model. In this article we consider models of the DFSZ type where one only enlarges the scalar sector (possibly adding also right-handed neutrinos). In a recent paper [17], the authors presented an invisible axion model where the PQ symmetry is not family universal but rather a horizontal symmetry. As the PQ symmetry cannot be used now to implement NFC, potentially dangerous tree-level FCNCs might be present. The approach followed in this case to avoid large flavor violating scalar couplings was to implement the flavored PQ symmetry in the same fashion as in the Branco-Grimus-Lavoura (BGL) model [18]. This way FCNCs appear at tree-level but they are controlled by the Cabibbo-Kobayashi-Maskawa (CKM) [19] and the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) [20] matrices. This axion model is characterized by several interesting features. Among these, we stress the possibility of avoiding the domain wall problem [21][22][23], also the presence of Flavor Changing Axion Interactions (FCAI) can introduce experimental constraints stronger than the astrophysical ones in some cases [24]. Invisible axion models with a horizontal PQ symmetry have been built previously in Refs. [25][26][27][28] and in the context of horizontal gauge symmetries in Ref. [29].
The present work is devoted to the extension and detailed analysis of the model presented in Ref. [17]. In Sec. 2 we introduce the notation and briefly review the main aspects of the BGL model in the two Higgs scenario. Sec. 3 is dedicated to the determination of the required conditions for this symmetry to be a chiral PQ symmetry, therefore demanding it to be QCD anomalous. We show that this condition cannot be fulfilled in the two-Higgs-doublet BGL model and that an extension of the scalar sector is required. In Sec. 4 we present a three-Higgsdoublet implementation with a horizontal PQ symmetry enforcing a BGL-like suppression in the FCNCs, we refer to this as the three-Higgs flavored Peccei physics signatures related to the Higgs sector can also be found in this section. We summarize our results and conclude in Sec. 8.

Notation and the Branco-Grimus-Lavoura model
In this section we introduce the notation used throughout the article and present the so-called Branco-Grimus-Lavoura (BGL) model [18]. We consider a 2HDM with the Higgs doublets parametrized as Here v 1 > 0 and v 2 > 0 generate the quark masses. We also set α 1 = 0 and α 2 ≡ α without loss of generality. Due to the presence of an additional Higgs doublet the Yukawa Lagrangian takes the general form whereΦ j = iσ 2 Φ * j , σ 2 being the Pauli matrix. In order to study some of the new phenomena present in this framework it is convenient to work in the Higgs basis, where the Goldstone bosons G + and G 0 are singled out and only one Higgs doublet acquires a non-vanishing vev.
For that, we perform the following transformations with Expanding the Yukawa Lagrangian in the Higgs basis one obtains The matrices M q and N 0 q (q = u, d) encode the flavor structure in the 2HDM, these are given by and The quark mass matrices M u,d determine the Yukawa couplings of the scalar field H 0 while the matrices N 0 u,d determine the Yukawa couplings of the scalar R and the pseudoscalar I. Note that the fields {H 0 , R, I} are not mass eigenstates in general, the physical neutral scalar bosons will correspond to a linear combination of these fields.
The quark mass matrices can be diagonalized through the bi-unitary transformations: chosen appropriately so that These transformations guarantee diagonal quark couplings for H 0 but, in general, so that R and I have flavor violating couplings at tree-level. These are the sources of dangerous FCNCs at tree level in the 2HDM. The most common solution to this problem is the NFC condition. This scenario is nothing more than the requirement of simultaneous diagonalization of M u,d and N 0 u,d , or equivalently, the simultaneous diagonalization of the Yukawa textures in each sector. In the two Higgs doublet models the NFC condition can be implemented in two ways: • Through a discrete or continuous symmetry which restricts the number of Yukawas in each sector to one [10].
• Using the alignment condition, where the Yukawa matrices in the same sector have the same flavor structure up to an overall factor [30]. This can be seen as an effective theory of a larger model with the first condition imposed at the UV level [31]. It can also be seen as a first order expansion in a minimal flavor violating scenario [32][33][34].
There exists, however, a different scenario where NFC is only imposed in one sector and FCNCs present in the other sector are under control, this is known as the BGL model [18]. The model uses an abelian symmetry to impose the Yukawa textures This implementation is also known as the top-BGL, since it singles out the top quark. In this case the flavor matrices responsible for the FCNCs take the form with V = U † uL U dL the CKM quark mixing matrix. This simple implementation introduces no FCNC effects in the up-quark sector. In the down-quark sector, one has tree-level FCNCs, however, those are highly suppressed. We can see from Eq. (14) that the second term of N d will introduce FCNCs which are suppressed by: • The down-type quark masses.
• The off-diagonal CKM matrix elements.
This way, the model implements controllable FCNCs at tree level within the 2HDM. As shown in Refs. [35,36] this implementation is unique, up to trivial permutations, in models with abelian symmetries. A detailed phenomenological study of the experimental constraints on this model was presented in Refs. [37][38][39].
Although the BGL model presents several unique features, it still suffers from a few problems. The first problem is present in the scalar potential of the model. While the abelian flavor symmetry used to get the desired textures can be implemented through a discrete group, the scalar sector will exhibit an accidental global U(1) symmetry leading to a Goldstone boson after spontaneous symmetry breaking [18]. Alternatives to eliminate the accidental global symmetry have been discussed in Ref. [18], adding soft breaking terms to the scalar potential or extending the scalar sector with gauge singlet fields could protect the model against the phenomenologically dangerous Goldstone modes. On the other hand, the strong CP problem is not addressed in this scenario. While there are no large contributions to electric dipole moments in the BGL model [40], this is based on the assumption of a vanishing or very small θ term [3].
In this work we suggest that these apparent problems of the BGL model could be solved in an unified way if the required Yukawa textures are imposed by a global chiral U(1) PQ symmetry, bringing also other advantages we will discuss in the following sections. The model then provides a dynamical solution to the strong CP problem via the PQ mechanism while an axion appears in the spectrum, which could in principle account for the dark matter of the Universe.

The anomalous condition for a BGL-like model
In this section we shall find the anomalous condition for the abelian continuous symmetry that imposes the BGL Yukawa textures, extending the analyses done in Ref. [17]. We are then interested in finding the abelian generators under which the fields must transform, i.e.
with S L =diag(e iX uL θ , e iX cL θ , e iX tL θ ) , S d R = diag(e iX dR θ , e iX sR θ , e iX bR θ ) , and with These field transformations will induce the following constraints with k = 1, 2. The Yukawa texture patterns are dictated by the way the fermion fields transform, the Higgs field transformation will only select one of the allowed textures [33,35,36,40].
The best way to find these fermion transformations is to study the Hermitian combinations Γ k Γ † k and Γ † k Γ k (and similarly for ∆ k ). The symmetry constraints on these combinations give The above equations are nothing more than the commutation of a diagonal matrix S L,R with a Hermitian matrix. In order for these matrices to commute S L,R must share the same eigenvectors as the Hermitian combination, or have degenerate eigenvalues for the non-shared eigenvectors. We then get three scenarios: i) The matrix S L,R has only one phase. The Hermitian combination has no restriction; ii) The matrix S L,R has two different phases. The Hermitian combination must be block diagonal, with the 2 × 2 block in the same sector as the degeneracy in S L,R ; iii) The matrix S L,R has three different phases. The Hermitian combination must be diagonal.
From Eq. (13) we see that the ∆ j Yukawa textures are block diagonal in the up-charm sector. The hermitian combinations ∆ k ∆ † k and ∆ † k ∆ k will also share the same form. This in turn implies that the symmetry generators for the left-and right-handed fields must belong to case ii), i.e. the abelian generators have only two different phases where we set one of the charges to zero using a global phase transformation. Notice that the charges should satisfy the conditions Condition A: X tL = 0 and in order to stay in the scenario ii). When the left-handed quark doublet and the right-handed up-type quark transform under this symmetry, the phases appearing in the Yukawa term are with the additional condition To the matrix Θ u we call the up-quark phase transformation matrix. The generators in Eq. (21) This choice makesΦ j associated with the ∆ j of Eq. (13). We can now build the phase transformation matrix for the down-quark sector. The left-handed transformation is the same, since it is shared by the two sectors. Concerning the right-handed generator of Γ j (see Eq. (13)), it belongs to the case i) and, therefore, has the form The down-quark phase transformation matrix is then given by Eqs. (21) and (26) together with the first part of condition A are the minimal set of required conditions necessary to obtain the BGL textures in the down sector. In order to pick the desired textures we would need the scalar transformation To make the BGL textures in the up and down sectors compatible without introducing additional textures that spoil the nice features of the BGL-type models we need to impose some extra charge conditions, we call them texture matching conditions. They guarantee that the only non-BGL textures present in the Yukawa sector are the null textures, Since we introduced a chiral symmetry we get the anomaly free condition for the PQ symmetry with the QCD currents This anomaly free condition makes both S up Φ and S down Φ equal for X dR = −X uR , making the BGL implementation consistent and anomaly free with two Higgs doublets. If we want this symmetry to be anomalous then X tL = −1/2(X uR − X tR ) and the model must be extended.
We shall pursue a possible anomalous implementation in the multi-Higgs framework, making the three-Higgs doublet model the minimal extension. In the three Higgs implementation we can just join the scalar generators S up Φ and S down Φ into a single one In this three-Higgs doublet model implementation we get the following Yukawa textures with the charge constraints Texture Matching Conditions: Anomaly condition: Since we extended the Higgs sector, in principle it is no longer necessary to have Γ BGL 1 and ∆ BGL 1 coupling to the same Higgs doublet. We can have another three different implementations: . This implies −X dR = X tR − X tL ; However, the first two implementations violate the charge restrictions in Eq. (29) whereas the third one is a safe implementation and gives so we get the following Yukawa textures implementation For this symmetry to be anomalous and in order not to introduce additional textures that spoil the desired behavior of the model we need to guarantee, in analogy to the previous case, the following charge restrictions Texture Matching Conditions: Anomaly Condition: The BGL 2HDM model needs condition B in its anomaly free implementation. However, when extending it to a three Higgs scenario, with the possibility of null couplings, this condition no longer needs to be satisfied. Relaxing this condition by setting X tR = X uR + X tL , we get a new type of texture in the up sector (the combination of ∆ BGL 1 and ∆ BGL 2 ). Three new possible implementations become available: • New texture coupling to Γ BGL 1 . This implies X dR = −X uR ; • New texture coupling to Γ BGL 2 . This implies X dR − X tL = −X uR ; • New texture coupling to a null texture.
The first two cases violate the charge conditions in Eq. (29), this is the reason why there is no BGL 2HDM with this texture. However, the third possibility gives a safe implementation in a three Higgs scenario with the same scalar charge assignments as in the previous case (i.e. Eq. (34)). The Yukawa textures implementation is then given by For the symmetry to be anomalous and to guarantee that we introduce no additional Yukawa textures, the following charge conditions apply Texture Matching Conditions: Anomaly Condition: In conclusion, in this section we have shown that it is not possible to build an anomalous two-Higgs-doublet modelà la BGL and we have found three different implementations of the PQ symmetry for the three-Higgs-doublet model, up to permutations in the family or in the up-down sectors. These three cases are built from the generators in Eqs. (21) and (26). They read as follows: • Case I: where the Yukawa textures are given by Eq. (32), satisfies conditions A and B, and also the texture matching and anomaly conditions in Eq. (33).
The charges associated with the Higgs fields are • Case II: with the Yukawa textures shown in Eq. (35), satisfies conditions A and B, and also the texture matching and anomaly conditions in Eq. (36).
The charges associated with the Higgs fields are • Case III: with the Yukawa textures shown in Eq. (37), satisfies the constraint X tR = X uR +X tL , condition A, and also the texture matching and anomaly conditions in Eq. (38).
The charges associated with the Higgs fields are the same as in case II.

The three-Higgs-doublet class of anomalous models
In the previous section we have shown that the Yukawa textures in the BGL 2HDM cannot be imposed by a chiral PQ symmetry. We also derived the necessary conditions to build three-Higgs doublet models with FCNC at tree-level completely determined by the fermion mixing matrices. In the latter scenario, we obtained all the possible Yukawa texture implementations imposed by a PQ symmetry and determined the restrictions to the PQ charges in each case.
We provide details about the quark Yukawa sector of these type of models in Sec. 4.1. The scalar potential of this class of models is discussed in detail in Sec. 4.2. The extension of the models considered to the leptonic sector is discussed in Sec. 4.3.

The Yukawa quark sector
In similar a fashion to what was done in Sec. 2, we shall build the relevant flavor matrix combinations that mediate the FCNCs. The Yukawa Lagrangian in the three Higgs scenario is now written as where we just keep the same notation as in Eq. (3), but for j = 1, 2, 3 in this case. We go once more to the Higgs basis, by preforming the following transformations with In the Higgs basis the mass and Yukawa interactions are given by with the flavor matrices given by and These last flavor matrix combinations are the ones mediating the FCNCs in our framework.
We can now evaluate them for each of the three cases. In the basis where the quarks are mass eigenstates we get: • Case I: • Case II: • Case III: As expected, in all cases the FCNCs will be mediated by quark masses and off-diagonal elements of the CKM quark mixing matrix. This is virtually the same type of suppression as the one obtained in the BGL 2HDM implementation. The difference lies in the vevs ratios that we get in front of each term. This actually contrasts with the anomaly free three Higgs BGL implementation [33]. In that scenario the Yukawa textures, which differ from the 2HDM implementation, cannot give such a strong suppression to |∆S| = 2 processes as compared to the original BGL implementation. One generally gets suppressions of the order of (V * cd V cs ) 2 ∼ λ 2 (λ 0.225), requiring heavy neutral scalar fields. However, the fact that we kept the same Yukawa textures in passing from the two to the three Higgs implementation allows us to have suppressions of the type (V * td V ts ) 2 ∼ λ 10 for |∆S| = 2 processes, just like the original BGL scenario.

The scalar potential
Current experimental limits exclude axions coming from a PQ symmetry broken at the EW scale [27,[41][42][43]. To obtain a viable axion model the PQ symmetry must be broken at a scale much higher than the EW scale. The axion is then called invisible since its mass and couplings are suppressed by the large PQ symmetry breaking scale. We can achieve this in a similar way as in the DFSZ and KSVZ invisible axion models, that is, by introducing a complex scalar singlet which acquires a very large vev 0|S|0 = The new complex field S will have the following symmetry transformation The introduction of the complex scalar singlet increases the number of independent charges in one unity. From the Yukawa sector alone, with fermion charges chosen in order for the symmetry to be anomalous, we are able to reduce the number of independent PQ charges to just three. In this way the number of independent charges increases to four.
The scalar doublets transform as in Eq. (17) with the charges X Φi expressed in terms of the three quark charges, their explicit form will depend on whether we are working in case I or II/III (as detailed in the previous section). We shall split the potential in two parts: the phase The phase blind terms do not constrain the charge assignments, they are given by The parameters λ ΦS i and λ ii,jj run for all i, j = 1, 2, 3, while the parameter λ ij,ji run for i = j.
This part of the potential possesses a U(1) 4 global symmetry. The role of the phase sensitive part is to introduce terms which break (explicitly) this symmetry down to With this symmetry we will have two complex phases to which the scalar potential will not be sensitive, one will be the neutral Goldstone boson and the other the axion. This will introduce two new additional constraints, reducing the number of independent charges down to two.
We shall now present the possible phase sensitive terms that we may built and their constraints in terms of the PQ charges. We note that any term of the form Φ † i Φ j (or any combination where this is the only phase sensitive part) implies the charge relation X Φi = X Φj , which is automatically excluded by the charge conditions, see Eqs. (39) and (40). Also, terms that are only sensitive to phases of one single field such as S k , Φ † i Φ i S k , etc. would imply a discrete phase, which is not allowed in our framework.
In Table 1 we present all the possible, renormalizable and gauge invariant, phase sensitive terms (up to hermitic conjugation). We now have to check all the possible combinations of two terms from (1) to (6). Combining just the first three cases will lead to a constraint of the type X Φi = X Φj , which is excluded. When combining cases (1) to (3) with cases (4) to (6) all of these last three cases will be allowed simultaneously. After finding all the possible combinations and using the information about the explicit forms of X Φi in terms of the quark charges we get the following charge constraints: • Case I: • Case II/III: We must also have C I = 0, −1, −1/2, C II = 0, 1, 1/2 and C III = −1, 0, 1 2 , 1, 3 (see Eqs. (33), (36) and (38), respectively) in order to preserve the Yukawa textures and the symmetry to be anomalous. In Table 2 we present all possible values for C I,II,III and C I,II,III S in each possible phase sensitive potential implementation.

Term
Combination  At this point we have two free charges which we choose to be X uR and X S , for all cases. We can normalize all charges to the scalar singlet charge, without loss of generality, just by setting the condition X S = 1. This allows the PQ quark charges to be written in terms of the values C I,II,III S , C I,II,III and one free charge, X uR . They will now take the form: • Case I: • Case II: • Case III: In this section we have found up to 12 distinct phase sensitive potential implementations, see Table 2. For case I, T 7 and T 10 implementations are not compatible with the flavor PQ symmetry in the fermionic sector. In case II, the incompatible implementations are T 9 and T 12 . Finally, case III has the same incompatible implementations as case II plus T 1 , T 8 and T 11 implementations. As an illustrative example, let us choose the implementation T 2 . The scalar potential would take the form with λ dimensionless and µ with mass dimension. Under this particular potential implementation, and with our normalization, the PQ quark charges read While the scalar charges are: The fact that the scalar charges are the same for all the three cases should not be surprising. The scalar potential itself knows nothing about the distinct Yukawa implementations, that information enters only when we use the explicit expression of the scalar charges in terms of the quark ones.
Therefore, the scalars charges will only depend on the the distinct potential implementations.

The Yukawa leptonic sector
In this section we shall only be interested in the Yukawa couplings of the charged leptons and therefore we will say nothing on the Dirac or Majorana nature of the neutrinos. We will assume that the final neutrino mass matrix texture contains enough freedom, such that, in combination with the lepton mass matrix accommodates the full low-energy neutrino data. However, note that the neutrino Yukawa textures should satisfy some conditions such that the BGL quark and lepton textures are not spoiled through radiative corrections [44]. In Ref. [17] a particular model implementation has been presented were the neutrino sector has been worked out. However, since we are mostly interested in the axion properties in this class of models, we can just focus our attention to the charged lepton implementation.
The Yukawa leptonic Lagrangian will be of the form In a similar way as it happens in the quark Yukawa sector, it is convenient to rewrite the Yukawa lepton Lagrangian by rotating the Higgs doublets to the Higgs basis (see Eq. (42)) and by diagonalizing the lepton mass matrices through the bi-unitary transformations The Yukawa Lagrangian now reads as where, as it happened with the quarks, N e and N e will mediate the FCNCs. These flavor combinations will have the same expression, in the flavor basis, as N d and N d present in Eq. (46) with the replacement Γ i → Π i .
Regarding the PQ symmetry transformations, the scalar field transformations are given in Eqs. (39) and (40) for cases I and II/III respectively. We now need to determine the PQ charges of the leptonic fields. In general these will transform under the continuous symmetry as with A global phase transformation allows us to set X eL = 0 without loss of generality, just as we did in the quark sector.
We could proceed with the symmetry implementation just like in the quark sector, however, we can also combine the BGL-like textures in the quark sector with NFC for the charged lepton such that we have several phenomenological models available. We shall then split these implementations into two classes: (1) With FCNCs in the charged lepton sector. This is the extension to three Higgs doublets of the symmetry implementation in Ref. [44].
In this case, in order to have the FCNCs under control we choose the implementationà la BGL, i.e.
Just like in the quark sector, we need the other sector mass matrix (i.e. neutrino mass matrix) to be block diagonal, in order to have the PMNS mediating the FCNCs. The way to achieve this will depend on the Dirac or Majorana nature of neutrino and is out of the scope of this paper (see Ref. [17] for more details). The symmetry implementation is just like the one in the quark sector, i.e. X eL = X µL ≡ X l L and X eR = X µR = X τ R ≡ X lR .
The constraints are The equivalent to conditions A and B in the quark sector also apply to the lepton charges.
Since we have set X eL = 0, the charged lepton charges become completely defined by the known scalar charges, i.e.
In this case there are six implementations possible, as it was shown in Ref. [36]. Using the information that all the charges of the scalar fields are different we get (a) In this scenario both left and right generators must be fully degenerate, i.e. X eL = X µL = X τ L ≡ X lL and X eR = X µR = X τ R ≡ X lR . This implies the following constraint In this scenario both left and right generators must be two-fold degenerate, i.e.
X eL = X µL ≡ X l L and X eR = X µR ≡ X l R . This implies the following constraints (c) In this scenario the left and right generators have the same form as in the previous one. However, the constraints are (d) In this scenario the left and right generators must have no degeneracy. The constraint is given by (e) In this scenario the left and right generators are the same as before. The constraints are given by In this scenario the left and right generators are the same as before. The constraints are given by In general, we have only information on the difference between left-and right-handed charged lepton charges. The condition X eL = 0 allows us to have the charged lepton charges fully determined by the known scalar charges only in cases (1) and (2a). For the other cases we would need to know the neutrino sector implementation. Nevertheless, as we shall see in the next section, the knowledge of the difference is enough to get most of the axion properties.

Axion properties
The anomalous U(1) PQ symmetry of the class of models built in the previous sections is spontaneously broken by the vev of the singlet field S at a very high scale, just like in the standard DFSZ and KSVZ models. Non-perturbative QCD effects induce a potential for the axion field, allowing us to shift away the strong CP phase and also give a small mass to the axion [9], the physical one (denoted by a). In the following we derive the most relevant axion properties for our model. We start by writing the relevant Lagrangian for the physical axion where L aγγ is the axion interaction to photons, we will shown in Sec. 5.1, and L aψψ the axion interaction to matter, we will present in Sec. 5.2. The axion mass is given by [9]: with m π 135 MeV and f π 92 MeV the pion mass and decay constant, respectively. The parameters z and w denote the quark mass ratios z = m u /m d 0.56 and w = m u /m s 0.029.
The quantity C ag is determined by the chiral color anomaly of the current associated with the U(1) PQ transformation [45], in our model it is given by This quantity turns out to be independent of the free charge and can be expressed as C M ag (with M = I, II, III) and is given by The quantity C M ag is therefore only dependent of the scalar implementation once the Yukawa textures are specified.

Axion-photon coupling
The axion two-photon interaction is described by the effective Lagrangian with α = e 2 /4π 1/137, F µν is the electromagnetic field strength tensor and F µν its dual tensor. The effective factor C eff aγ takes the form [46]: where the second term is a model independent quantity which comes from the mixing of the axion with the π 0 and the η while C aγ and C ag are model dependent quantities associated to the axial anomaly. These are determined in terms of the fermion charges by while C ag was already introduced in Eq. (81). The quantity C aγ can be expressed as with M = I, II, III. Here we have introduced the parameters C I,II,III and C I,II,III S specified in Table 2 and a new combination of fermionic charges A I,II,III defined in Table 3. The charged lepton combinations are denoted by a vector (i, j, k), which represents the Higgs doublet that is coupled to the left-handed leptons (e, µ, τ ). For example, the case (1, 1, 3) tell us that Φ 1 is coupled to e L and µ L while Φ 3 couples to τ L . This can correspond to the charged lepton Yukawa implementations (1), (2c) or (2e). Note also, that the case (3, 3, 1) is not a relabeling of the scalar fields since we keep the quark sector unchanged and, therefore, we will get a distinct result. Once the choice on the Yukawa textures is made, the parameter C M aγ will depend on the potential implementation and the way the charged leptons transform under the PQ symmetry.

Cases
A

Axion couplings to matter
The interactions of the axion with fermions are described by with η = (1 + z + w) −1 . When calculating the axion couplings to matter one should redefine the axion current in such a way that it does not mix with the neutral Goldstone boson associated with the spontaneous symmetry breaking of the electroweak gauge symmetry. This redefinition results in a shift of the original scalar charges [46], which in terms of the fermion charges reads with = {e, µ, τ } and The explicit expression for Z will take the same form in cases II and III (they share the same Higgs charge assignments) but a different form in case I, i.e.
for case II and III .
We now define the shifted charge matrix as which take the explicit form X uR =diag(X uR , X uR , X tR ) , X dR = X dR I , X eR = diag(X eR , X µR , X τ R ) .
These charge matrices determine the couplings of the axion to fermions in the flavor basis. By going to the mass basis the fermion fields are rotated through the unitary transformations in Eq. (10) and in Eq. (60). These transformations will change the charge matrix to In this new basis the quark charge matrices take the explicit form where we have used X tL = X tL . For the charged leptons we have in scenario (1) where we have used X τ L = X τ L , in scenario (2) we get The axion vector and axial couplings to matter are then given by From the above equations we can see that the flavor changing axion interactions will be mediated by the off-diagonal elements of X dL (and X eL in case (1)). This is a common property of Goldstone bosons in flavor models, however it is an additional feature for the axion compared to the standard DFSZ and KSVZ scenarios.
Regarding the vectorial couplings, it is interesting to remark that the PQ symmetry implementation that we have used until now is defined up to a global vectorial phase transformation which allows us to remove some of the vectorial couplings. As a result, the Lagrangian will remain invariant if we redefine the PQ charges by performing the following transformation with X X defined in Eq. (93) and α an arbitrary constant. In the DFSZ model this transformation is enough to remove all the vectorial couplings. However, this is not the case in the models we are presenting. For example, by setting α = −X uR /2 the transformed quark charges read as such that now C V au 11 = C V au 22 = 0. Similarly, we could have set α = − (X tL + X tR ) /2 to make C V au 33 = 0 but there is no value of α that makes C V au = 0 for the three families simultaneously. A similar procedure can be applied to the lepton sector.
Additionally, we can fix the value of the free PQ charge, X uR , in order to remove extra vectorial couplings. For instance, setting X uR = X tR + X tL in Eq. (100) would also make C V au 33 = 0 for α = −X uR /2. However, the condition X uR = X tR + X tL can only be satisfied in cases I and II while in case III it would violate the charge constraints in Eq. (22). The same happens for other values of α and thus one cannot use the freedom in X uR to remove vectorial couplings in case III.
Finally, note that the vectorial transformation and the freedom to fix the value of X uR only affect the diagonal vector couplings while the off-diagonal ones remain unchanged. In any case, it is simple to see that in the case of on-shell fermions the axion-fermion interaction is purely pseudoscalar for fermions of the same flavor with ψ representing a fermionic specie (up quarks, down quarks and charged leptons). Therefore, the nature of the axion interaction in the up quark sector is purely pseudoscalar for on-shell quarks. However, due to the presence of FCNCs in the quark sector, the axion interaction will no longer conserve flavor, reflecting a scalar (beside the pseudoscalar) nature of the axion field in models with FCNCs.
The axion axial couplings to light quarks are explicitly given by for case II and III , for case II and III , Below the chiral symmetry breaking scale the axion nucleon interactions can be parametrized by with σ 3 the Pauli matrix in the isospin space and N = (p, n) T the nucleon doublet. The isoscalar and isovector couplings are given in Refs. [26,46]. The couplings to protons and neutrons are given by the combinations g p ≡ g 0 + g 3 =(g u − 2ηC ag )∆u + (g d − 2ηC ag z)∆d + (g s − 2ηC ag w)∆s , g n ≡ g 0 − g 3 =(g u − 2ηC ag )∆d + (g d − 2ηC ag z)∆u + (g s − 2ηC ag w)∆s , with ∆u = 0.841 ± 0.020, ∆d = −0.426 ± 0.020 and ∆s = −0.085 ± 0.015 [47].
The coupling to electrons in scenario (1) is given by The scenario (2) can be obtained in the limit |V τ 1 | 2 → 0. The explicit form of g e will depend on the scalar potential implementation and the vevs of the doublet fields.

The domain wall problem
During the evolution of the Universe the PQ symmetry gets broken in different ways. In the early Universe the PQ symmetry is spontaneously broken by the expectation value of the S field.
At this stage the potential has the mexican-hat shape, the angular part of the field becomes a Goldstone boson and the Lagrangian remains global phase invariant. As the Universe cools down non-perturbative instantonic effects at the QCD scale take place and the PQ symmetry gets explicitly broken by the QCD gluon anomaly [SU(3) C ] 2 × U(1) PQ [9]. This is the PQ mechanism for the resolution of the strong CP problem.
However QCD instantons only break the symmetry down to a discrete Z N subgroup. This can be easily seen from Eq. (1) and the fact that the θ term is invariant under the shift θ → θ + 2πk. While before the QCD scale the shift a/v PQ → a/v PQ + α was allowed for any α, the presence of the axion coupling to gluons restricts the phase to the values α k = 2πk/|C ag | (with k = 0, 1, . . . , |C ag | − 1), which just reflects the original θ periodicity. Therefore, the order N of the discrete group is given by the color instantons effects to be N = |C ag |. As pointed out by Sikivie [22], these models will have N DW degenerate disconnected vacua. This in turn leads to an unwanted domain wall structure in the early universe [21][22][23].
In general, the domain wall number, N DW , coincides with the order of the discrete group, i.e. N DW = N = |C ag |. Nonetheless, in some cases only a subgroup of Z N acts non-trivially on the vacuum. To examine the vacuum structure one should analyze the gauge invariant order parameters of the theory. In this way the domain wall number will coincide with the dimension of the higher order subgroup of Z N which acts non-trivially on at least one of the order parameters. For the models we are discussing, it suffices to notice that for the singlet as we set X S = 1 the vacuum periodicity is N , i.e. all elements of the residual Z N act nontrivially on the vacuum. As a result, we have N DW = |C ag | for the class of models studied in this article. 1 Many axion models suffer from the domain wall problem. In particular, the well known DFSZ invisible axion model has a domain wall number N DW = 2N g or N DW = N g depending on the scalar potential implementation, with N g the number of quark generations. In Table 4 we present the values of the domain wall number for the different implementations of the models we are presenting. As we can see, while some of the implementations also suffer from the domain wall problem, others have N DW = 1, for which the resulting domain wall structure is harmless [48].  Even for N DW = 1, some solutions to the domain wall problem can be found in the literature.
It is possible to avoid the domain wall problem by assuming that inflation has occurred after the PQ symmetry breaking. In this case, one can derive limits on the inflationary scale based on the observation of isocurvature fluctuations in the cosmic microwave background [49]. Also, there have been several attempts to introduce an explicit breaking of the PQ symmetry that also breaks the Z N discrete group in such a way that the PQ solution to the Strong CP problem is protected [22,50]. This explicit breaking could come from gravity, giving rise to Planck scale suppressed operators [51,52]. However, it was argued that this solution would give rise to long lived domain walls which introduce cosmological problems [53]. Additionally, gravity violations of the PQ symmetry should be controlled in order not to spoil the PQ solution [54][55][56][57][58][59][60]. This issue will be considered in the next section.
In the absence of gravity, the axion potential coming from the instantonic contributions can be written as [8] V axion −Λ 4 QCD cos which has a minimum at a phys =θ = 0 and where the estimated axion mass is m a Λ 2 QCD /v PQ . On the other hand, when including gravitational effects, the axion potential will change. For example in our invisible axion model, we should expect higher dimensional PQ violating terms of the type 1 M n−2 with the Planck scale denoted by M Pl . Let us consider, for simplicity, the second term in the above equation. By introducing this term in the Lagrangian the axion potential gets modified and takes the form [58] V axion −Λ QCD cos The parameter c is just a coupling constant and δ a CP violating phase. The problem with this new axion potential is that the minimum is no longer at a phys = 0, but rather at which in general will be far from zero. The axion mass will also be affected by these gravitational effects, taking the form Therefore, in this simple scenario we see that gravitational effects will in general spoil the strong CP solution coming from the PQ symmetry.
Fortunately, over the years several solutions to this problem have emerged, which allows us to preserve the PQ solution of the strong CP problem. GUT motivated models [57,58], extra dimensional [61] ones and even models having neutrinos playing a big role in gravity protection [62] can be found in the literature. However, many of these solutions need a significant extension of the original PQ model.
In this section we will focus on the use of gauge discrete symmetries to protect the PQ solution against gravity [63]. This solution has the interesting feature that the low energy spectrum of the theory does not need to be extended. Gauge discrete symmetries, which arise through the spontaneous symmetry breaking of a gauge symmetry, are not broken by gravity and can provide natural suppression to the harmful gravitational effects. The idea proposed is to have a large discrete abelian symmetry Z P forbidding, up to a given order, these unwanted terms [64,65]. For example, if the symmetry only allows terms of the form S m /M m−4 Pl for m ≥ 13, we will just get irrelevant contributions to the axion mass and its vev [66].
In what follows we will identify the PQ symmetry as an accidental global symmetry at low energies, associated with the spontaneous breaking of a gauge symmetry, U(1) A , at high energies down to a discrete subgroup. We shall follow Ref. [64], using a discrete gauge symmetry to stabilize the axion without enlarging the low energy particle content. To this end, we shall use the discrete version of the Green-Schwarz anomaly cancellation mechanism [67].
From the effective theory point of view, since we have at low energies the SU(3) C ×SU(2) L × U(1) Y gauge group, there are several possible anomalies we must consider: The Green-Schwarz anomaly cancellation conditions are then given by with δ GS a constant that cannot be specified by the low energy theory and k 1,2,3,A the levels of the Kac-Moody algebra [68] which are integers for non-abelian groups. The equality involving the hypercharge currents, A 1 , give no useful constraints since the associated level k 1 is not an integer in general [63]. Similarly, the anomaly A A can be canceled by the Green-Schwarz mechanism but with no useful constraints due to the arbitrariness in the normalization of U(1) A .
Finally, the anomaly A G give no useful constraint either.
When the U(1) A is broken down to a Z P , the effective low energy theory will satisfy the discrete version of the Green-Schwarz cancellation condition [63,69] with m and m integers. The model under discussion is non-supersymmetric, nevertheless, the Green-Schwarz mechanism should still be available since the breaking of supersymmetry can happen at the scale much higher that the weak scale.
Our goal is to build a U(1) A symmetry that contains a discrete subgroup capable of solving the strong CP problem. The U(1) PQ group is anomalous and, therefore, capable of giving such a solution (as it was seen in the previous sections). However, U(1) PQ cannot be identified with U(1) A as the PQ symmetry alone is, in general, not enough to satisfy the Green-Schwarz anomaly conditions. Fortunately, the model also presents Baryon number conservation (+1 charge for quarks, −1 for anti-quarks), which is QCD anomaly free but it is anomalous under SU(2) L . We shall then try to see if the combination U(1) PQ + γU(1) B is suitable to be our axial symmetry. As the lepton charges depend on the specific representation in the neutrino sector, we will focus on the quark sector. The generalization to the lepton sector will be discussed at the end. From the U(1) PQ charge assignments in Eqs. (55), (56) and (57) we can find the anomaly coefficients with M = I, II, III. The factor γ is then found to be Using the simplest realization of the Kac-Moody algebra, i.e. k 2 = k 3 = 1, we get for each case Normalizing the combination to have all charges integer number we then define the axial abelian symmetry as The charges under this new symmetry are given in Table 5. Finally, note that the inclusion of PQ lepton charges would modify A 2 in the following way while A 3 would remain unaltered. This accounts to a correction of γ of the form which transforms the quark charges in Table 5 to C III S C III S 9 9 9 and leaves the lepton and scalar charges unchanged.
In Table 6 we present the axial symmetry and its discrete Z 13 version in the T 6 scenario, for each case I, II and III. The discrete anomaly coefficients are A 3 = −2 and A 2 = 24 for case I, A 3 = −15/2 and A 2 = 12 for case II, and A 3 = 3 and A 2 = 45/2 for case III. In each case, by construction, the anomaly coefficients satisfy the discrete Green-Schwarz cancellation condition Eq. (114).
In this example the phase sensitive terms are explicitly given by Due to gravity effects we expect the most relevant U(1) PQ breaking contributions to be of the with k integer. The largest contribution will be the one coming from the O 0 operator. Due to the Z 13 symmetry this contribution will only take place for k = 13, i.e. S 13 /M 9 P l . This operator will give a contribution to the axion mass squared of the order v 11 PQ /M 9 P l ∼ [10 −72 , 10 −39 ] GeV 2 for Case II: Case III: Z 13 2 11 x + 2 x + 11 −x + 3 x x + 12 x + 8 9 Table 6: Particular example with the phase sensitive scalar potential T 6 .
These are extremely small contributions, making the model safe against large gravitational corrections.
In this section we have shown how we could avoid large contributions to the axion mass, as well as to theθ parameter, just by invoking a discrete gauge symmetry. However, there are many more effective operators that will be induced by gravity than those presented above. this additional contribution. This is not so strange taking into account that in our framework no information on the Yukawa hierarchy is given. We know that O(y u , y d , · · · ) 1 and in our framework this is imposed by hand. In a more complete model, these hierarchies could be made dynamical and there we should also take attention to these additional gravity induced terms.

Model variations
The models presented in the main sections of the paper had FCNCs in the down-quark sector and the top quark was singled out. However, there are many other possible implementations that will still give the same minimal flavor violating scenario. These model variations can be found by performing any of the two operations: (i) Symmetric permutations in the flavor space; (ii) Changing up and down right-handed generators.
We can apply these two operations to the models previously studied in order to get all possible model variants.

Type (i) operation
The permutations in flavor space will change the textures in the sector with no FCNCs, i.e the up sector if we apply this operation in the original formulation. The symmetry generators take now the form with P and P 3 × 3 permutation matrices. The 2 by 2 block in the NFC sector will change structure, we get Where P ij permutes the lines i and j (when applied on the left) and columns i and j (when applied on the right). Besides the textures the only changes due to (i) are in the axion-matter couplings. The permutation matrices single out other flavors. Therefore, the action of the permutation matrices will change the CKM and PMNS elements entering in Eq. (95) and Eq. (96), respectively. We get the following redefinition Consequently, the couplings u, d, s and e are appropriately changed.

Type (ii) operation
We change the sector where the FCNCs are present, that can be accounted with the following symmetry generators S L = diag(1, 1, e iX tL θ ) , S u R = e iX dR θ I , S d R = diag(e iX uR θ , e iX uR θ , e iX tR θ ) .
We have switched the S u R and S d R generators, keeping the same labels for the charges. Thus, in this scenario the S u R is completely degenerate, but instead of labeling the charge X uR we keep it labeled as X dR , just as in the previous case. By keeping the same label we can easily compare this new case with the previous scenario where the FCNCs where in the down sector. The Higgs charges, for the three cases, take the same form as in the original scenario (see Eqs. (39) and (40)) but with an overall minus sign.
The quark charges take the same form as in Eq. (55), (56) and (57) as long as in the scalar sector the role of the S field is substituted by the S * , keeping the X S = 1 normalization. This will account for the overall minus sign coming from the Higgs charges. While the lefthanded quark charges have the same expression as in the original scenario, the right-handed ones switched sectors. The coupling to gluons will not change, however the coupling to photons changes since the up and down electric charges are different. It will be given by This will have the same form as Eq. (86), but now the coefficients take the form given in Table 7 Cases Table 7: Charge combinations A I , A II and A III entering in the description of the axion coupling to photons. The numbers in the first column label the Higgs doublet giving mass to the charged leptons (e, µ, τ ); for example (1,1,2) stands for the case where Φ 1 gives mass to e and µ while Φ 2 to τ .
The axion coupling to matter will also change, since we have changed the FCNC sector. The shift we need to perform in order to account for the orthogonality between the axion and the Goldstone boson will get an overall minus sign with respect to Eq. (91). Because we changed sectors without changing the labels of the charges, the charge shifts will coincide with the ones in Eq. (89). The shifted charge matrix takes the form X uL =X dL = diag(0, 0, X tL ) , X eL = diag(X eL , X µL , X τ L ) , X uR =X dR I , X dR = diag(X uR , X uR , X tR ) , X eR = diag(X eR , X µR , X τ R ) .
In the mass basis they take the explicit form The axial couplings to the light quarks are now given by for case II and III , for case II and III . (132)

Model variations dictionary
In this section we present a compilation of the most significant changes resulting from applying the operations (i) or (ii) to the original formulation. These can be found in Table (8).

Discussion
In the previous sections we have characterized a class of invisible axion models with tree-level FCNCs. We have also detailed the most relevant properties of the axion in these models. In this section we will analyze the different constraints on these models due to familon searches in kaon and muon decays, astrophysical considerations, as well as axion searches that rely on the axion-photon conversion mechanism. We will separate the discussion as follows: In Sec. 7.1 we discuss the constraints that can be extracted from the axion-photon coupling. In Sec. 7.2 we consider constraints on the axion from its flavor diagonal couplings to matter (nucleons and electrons). In Sec. 7.3 we discuss constraints derived from familon searches in rare kaon and muon decays. Our main results regarding the constraints on the axion are summarized in    frameworks considered, for this we feel it is important to discuss also the possible decoupling limits of the scalar sector.

Constraining the axion coupling to photons
The axion couples to photons through the dimension 5 operator in Eq. (83). The axion-photon coupling constant is defined by In various extensions of the SM, weakly coupled light pseudoscalar particles emerge naturally.
However, the axion possesses (due to QCD effects) an inherent correlation between the photon coupling and its mass (m a /1 eV) 0.5 ξ g 10 , where g 10 = |g aγ |/(10 −10 GeV −1 ) and ξ = 1/|C eff aγ |. The dimensionless coefficient ξ is in many axion models of order 1. In the well-known DFSZ (type II and flipped) and KSVZ models ξ takes the approximate values 1.4 (0.8) and 0.5, respectively. 2 In our 3HFPQ scenario we can get, besides the standard values, an additional set of discrete values allowing us to cover a large range of axion-photon couplings. We present in Fig. 1 the values of the model dependent quantity C aγ /C ag for those 3HFPQ models with N DW = 1, where the different values for each potential implementation represent distinct models for the charged leptons. It is well known that the evolution of stars place strong constraints on the axion coupling 2 For the KSVZ model we have taken the benchmark of X em Q = 0, with X em Q denoting the electric charge of the exotic color triplet Q. For the DFSZ model there are two possible implementations of NFC, the Higgs doublet coupling to l R can couple either to down-type quarks (type II) or to up-type quarks (flipped).
to photons. A strong bound can be derived from globular-cluster stars [70]. These are homogeneous gravitationally bound systems of stars formed around the same time, allowing for detailed tests of stellar-evolution theory. The actual experimental bound gives g 10 < ∼ 1 for axion masses up to 30 keV. Recently, the analysis of the evolution of massive stars lead to the bound g 10 < ∼ 0.8, based on the fact that Cepheid variable stars exist [71]. An even more stringent bound from an updated analysis of 39 Galactic Globular Clusters has been reported [72], setting the limit g 10 < ∼ 0.66 at 95% CL.
Several helioscope and haloscopes experiments are currently involved in probing the g aγ coupling. The most powerful axion helioscope experiment is the CERN Axion Solar Telescope (CAST), which searches for solar axions via axion-photon conversion using a dipole magnet directed towards the sun. The CAST experiment achieved the limit g 10 < ∼ 0.88 for m a < ∼ 0.02 eV, while slightly weaker bounds were obtained for heavier axions [73]. Still the astrophysical bounds represent a slight improvement over the CAST results. It is expected that the next generation of axion helioscope experiments, such as the International Axion Observatory (IAXO) [74], will provide better bounds on the axion-photon couplings in the future.
Microwave cavity haloscopes, including the Axion Dark Matter experiment (ADMX), exclude a window for the axion around a few µeV [75][76][77]. These experiments search for cold dark matter axions in the local galactic dark matter halo.

Constraining the axion couplings to matter
We define the axion-electron coupling constant as h aee = |g e |m e /v PQ , with g e given by Eq. (105).
The axion-electron coupling is bounded from astrophysical sources. In globular clusters, energy losses in red-giant stars due to axion emission would delay helium-ignition and make the redgiant branch extend to brighter stars. Helium ignition arguments in red-giant branch stars place the following upper bound on the axion-electron coupling, h aee < ∼ 3 × 10 −13 [78]. A more restrictive bound comes from white-dwarf (WD) cooling due to axion losses [70,79], The bounds that can be extracted on the PQ symmetry breaking scale, or alternatively, on the axion mass, are very model dependent for this observable. The value of g e given by Eq. (105) not only depends on the particular charge assignments of the model considered but also on the vevs of the Higgs doublets. In some regions of the parameter space it is even possible to obtain g e 0 so that WD cooling arguments would not place a strong bound on the axion mass. 3 Taking the benchmark point |g e |/N DW = 10 −1 for example, implies the upper bound m a < ∼ 15 meV.
Axion-nucleon interactions are constrained by the requirement that the neutrino signal of the supernova SN 1987A is not excessively shortened by axion losses [70,80]. We find these constraints to be similar than those coming from the WD cooling arguments in general.
However, the SN 1987A limit involves many uncertainties which are not easy to quantify [70].
Once more, the bound extracted on the PQ scale from the SN 1987A will depend on the vevs of the Higgs doublets.
The axion couplings to matter can also be tested in terrestrial laboratories, with promising prospects of probing unexplored regions of the axion parameter space. Dark matter axions can cause transitions between atomic states that differ in energy by an amount equal to the axion mass. By tuning the atomic states energy using the Zeeman effect it is possible in principle to detect axion dark matter candidates in the 10 −4 eV mass range [81]. The axion can also be tested in dedicated laboratory experiments looking for oscillating nucleon electric dipole moments (EDMs) [82][83][84], and, oscillating parity-and time reversal-violating effects in atoms and molecules [83,84]. The proposed Cosmic Axion Spin Precession Experiment (CASPEr) for example, could cover the entire range of axion dark matter masses m a µeV by looking for oscillating EDMs in a nuclear magnetic resonance solid-state experiment [85].

Constraining flavor changing axion interactions
In our framework, and contrarily to what happens in the DFSZ and KSVZ models, the axion couples differently to different flavors and has flavor changing interactions at tree level. Pseudo-Goldstone bosons arising from the spontaneous breaking of a horizontal symmetry are known as familons [24,86]. In our class of models this familon is the axion, and it will have non- The leptonic decay µ + → e + a can in principle be used to constrain charged lepton flavor violating interactions of the axion. In Ref. [87] the authors reported the experimental bound Br(µ + → e + a) < 2.6 × 10 −6 at 90% CL. This result however relies on the assumption that the positron is emitted isotropically to avoid large backgrounds from the ordinary muon decay µ + → e + ν eνµ . This assumption would be valid if we only had vectorial couplings but not axial ones. In our scenarios, the lepton flavor violating axial and vectorial couplings are equal, so the assumptions behind the µ + → e + a bound do not apply. In this case the best process to bound the charged lepton flavor violating axion couplings is the radiative decay µ + → e + a γ. With this process it is possible to extract limits which are independent of the chirality properties of the axion couplings [88]. The most stringent experimental bound at the moment is Br(µ + → e + a γ) < 1.1 × 10 −9 at 90% CL [89], obtained at the Los Alamos Meson Physics Facility (LAMPF) using the Crystal Box detector. From this process we can extract with the axion-lepton couplings g V µe = g A µe = C A ae 21 in the scenario with FCNCs in the charged lepton sector. We obtain a robust bound from this process since the flavor changing couplings g V, A µe are completely determined by elements of the PMNS lepton mixing matrix due to the underlying PQ symmetry. The bound extracted on the PQ scale, or equivalently the axion mass, from µ + → e + aγ does not vary much between all the models with FCNCs in the charged lepton sector. The reason being that that the PMNS matrix is very anarchical, that is, Obviously, models without tree-level FCNCs in the charged lepton sector avoid the constraints coming from µ + → e + aγ.
Models with FCNCs in the up-quark sector do not receive strong constraints from flavor observables, all the relevant observables involve a double insertion of the axion couplings in this case. On the other hand, models with FCNCs in the down-quark sector are strongly constrained by limits on K + → π + a. To the best of our knowledge, the strongest bound on K + → π + a decays has been set by the E787 Collaboration at Brookhaven National Laboratory, achieving Br(K + → π + a) < 4.5 × 10 −11 at 90% CL [90]. The partial decay width for this process is given by Γ(K + → π + a) = 1 64π with β = 1 − m 2 π /m 2 K and g V sd = C V ad 21 . The relevant hadronic matrix element can be extracted in the limit of exact SU(3) flavor symmetry. At the zero momentum transfer the form factor has the fixed normalization F 1 (0) = 1 [91]. We have π + (p )| sγ µ γ 5 d |K + (p) = 0 because K + and π + are pseudoscalar mesons. From this result we can extract a lower bound Just like in the charged lepton sector, the coupling g V sd is fixed in terms of elements of the CKM quark mixing matrix due to the underlying PQ symmetry. A robust bound can then be extracted on the axion mass which is independent of the many free parameters of the model. Future improvements on the µ + → e + aγ bounds are difficult to achieve with present facilities, see discussion in Ref. [92]. On the other hand, improvements on the K + → π + a limits can be expected from the NA62 experiment at CERN [93].
In Figs. 2 and 3 we summarize all the constraints discussed so far on the axion properties.
We do not show explicitly the limits from massive stars on g aγ though these are similar to that from CAST. In Fig. 2 we show constraints on models with FCNCs in the charged lepton sector and in the down-quark sector which select the top-quark. The strongest bound from flavor observables arises in this case from µ + → e + aγ because of the strong suppression factor |V * ts V td | entering in K + → π + a decays. The wide yellow band represents the prediction scanning over all the 3HFPQ models of this type. For this type of models astrophysical bounds from WD cooling put in general a stronger limit on the axion mass than flavor processes, the WD bound however depends strongly on the vevs of the Higgs doublets while the flavor limits do not. This is precisely what occurs for the model analyzed in Ref. [17], which corresponds to a case I model with scalar potential implementation T 11 and leptonic implementation (3,3,2 In Fig. 3 we show the constraints on those models with FCNCs in the down-quark sector which select the up or charm quark. The most relevant limit on the axion mass comes now from K + → π + a due to the value of the product of CKM matrix elements |V * ud V us | 2 ∼ |V * cd V cs | 2 |V * td V ts | 2 lifting the decay rate. For some models of this type the bound from kaon decays can be as strong as m a < ∼ 2 × 10 −5 eV. This is one of the main results of our work. Among all the models with FCNCs in the down-quark sector, those which select the top quark are much less constrained because of the very effective CKM suppression entering in K + → π + a decays.

Higgs physics
The Another degree of freedom corresponds to the axion which, up to corrections of order O (v/v PQ ), is given by the phase of the scalar gauge singlet. The other 10 degrees of freedom become physical scalar fields, leaving 2 electrically charged and 6 neutral physical scalars. It is not our intent to present a detailed analysis of the Higgs phenomenology in this class of models. We will, nevertheless, say a few words on some of these aspects. However, before discussing the Higgs phenomenology it is necessary to have some basic grasp of the decoupling structure of the kind of models considered.

Decoupling in the scalar sector
After the scalar fields acquire a vev, mixing among the scalars with the same charge is induced.
Due to the large hierarchy between the vevs, i.e. v PQ v, the radial part of the gauge singlet acquires a large mass and we can treat the mixing as SU(2) L -conserving. The scalar potential of the SU(2) L doublets is then given by with the square mass matrix taking the form The couplings λ i are associated with a phase sensitive terms (i) of Table 1. In the case where a phase sensitive term has mass dimension, such as in cases (4), (5) and (6) with k i = 1, we parametrize it as µ i = v PQ λ i . Additionally, note that for models T 1 to T 9 the PQ symmetry forces one of the couplings to be zero, that is λ k = 0 with k = 4, 5 or 6, while for models T 10 to T 12 all the couplings are expected to be non-zero.
Because of the large value of the PQ symmetry breaking scale, the scalar sector of invisible axion models usually presents a decoupling scenario. This way only the axion and a SM-like Higgs remains at the EW scale while the others acquire a mass of order v PQ . However, it is possible to achieve specific values for the parameters in order to avoid the decoupling limit in such a way that two (or three) Higgs doublets get masses around the EW scale.
In what follows, we analyze the different decoupling limits and give a possible texture reproducing each scenario. In the textures we use the parameters b, c ∼ O(1) and ∼ O v 2 /v 2 PQ . The last parameter, , has been introduced in order to show how EW corrections coming from the terms we have neglected in Eq. (140) can lift the zero eigenvalues to the EW scale. We also distinguish between the case where one of the λ-couplings is zero as in models T 1 to T 9 and the case where all the couplings are non-zero, corresponding to models T 10 to T 12 .
• One doublet at the electroweak scale.
This scenario is characterized by the presence of a SM-like Higgs at the EW scale, with the other two scalar doublets having masses at the PQ symmetry breaking scale. As a result, the infrared theory will correspond to the SM plus the axion (whose properties and couplings were discussed in Sec. 5) supplemented by higher dimension operators suppressed by the PQ breaking scale which can be neglected. 4 We list two textures generating this scenario This texture can be implemented in models T 1 to T 9 , even though the zero has been located in the position of λ 5 . The same mass spectrum is generated by permuting the value of the parameters appropriately. On the other hand, for the models T 10 to T 12 one possible texture is given by where the constraint b = c needs to be satisfied to have just one doublet at the EW scale.
• Two doublets at the electroweak scale.
In this decoupling scenario we obtain a 2HDM with tree-level FCNCs controlled by the CKM and PMNS matrices, the quark masses and the ratio of the vevs of the two Higgs doublets whose masses are at the EW scale (in a similar fashion as in the BGL 2HDM).
However, the values of the flavor changing scalar couplings cannot be determined in general as it will depend on the specific implementation of the scalar parameters in M 2 .
In any case, as these couplings present the same structure as in the BGL models, similar constraints in the parameter space are expected.
For models T 1 to T 9 it is not possible to reproduce this scenario unless some of the parameters are ultraweak, i.e. of order O v 2 /v 2 PQ . In this case one possible texture is given by which is only valid for models T 1 , T 2 and T 7 . The equivalent texture for models T 3 to T 6 , T 8 and T 9 can be directly obtained from the previous texture by permuting the entries in the matrix. Finally, one texture reproducing this scenario in models T 10 to T 12 is • Three doublets at the electroweak scale.
Having the three doublets at the EW scale is only possible if we force all the parameters in Eq. (141) to be ultraweak, that is if every term in M 2 take values around the EW scale.
As we have discussed in Sec. 4, this scenario gives rise to FCNCs which are suppressed by the CKM and the PMNS matrices with the explicit scalar flavor violating couplings depending on the model implementation.
This simple analysis of the possible decoupling scenarios is by no means a full and detailed study of the scalar spectrum. The textures above are just illustrative and many other textures with different degrees of tuning might be present for any of the three relevant decoupling scenarios. Finally, it should be noted that the scalar sector in this class of models suffers from a fine tuning problem (commonly known as the hierarchy problem), just like most models where more than one scale is present in the theory. A solution for this problem is out of the scope of the present work. However, some promising directions have been pursued in the literature within the framework of invisible axion models [94].

Higgs phenomenology
If Neglecting mixing effects among the scalar fields, the phenomenology of these scalars would be basically the same than in the BGL 2HDMs analyzed in Refs. [37][38][39]. For example, dangerous |∆S| = 2 contributions to K 0 −K 0 mixing due to neutral scalars would be very suppressed in the top BGL models because the flavor changing couplings are proportional to |V * ts V td |, allowing the mass of these scalars to be at the weak scale [37].
A classification of flavor observables which receive important contributions in the BGL 2HDMs and a comprehensive phenomenological analysis of this models was presented in Ref. [37]. , τ → ππ and µ − e conversion in nuclei. The previous processes arise in the SM at the loop level and receive strong suppressions due to the GIM mechanism or the smallness of neutrino masses, this makes these processes very sensitive to small new physics contributions. Charged scalars will also contribute at tree-level to semi-leptonic pseudoscalar meson decays (M → ν, B → D ( * ) ν) and leptonic τ decays (τ → ν ν τ ), possibly causing observable violations of lepton universality. Neutral and charged scalars will contribute at the loop level in processes likeB → X s γ and 1 → 2 γ and will in general dominate over the SM contribution which appears at the same level. The discovery of additional scalars at the LHC and characteristic decay signatures of the non-standard scalars in the BGL 2HDMs have been analyzed in Ref. [38,39]. The main results of these analyses is that within BGL 2HDMs additional charged and neutral scalars can be as light as 150 GeV while being compatible with present 125 GeV Higgs, flavor, electroweak precision and collider data [37][38][39].

Conclusions
In this work we have built a class of invisible axion models with FCNCs at tree level which are controlled by the fermion mixing matrices, therefore extending the work done in Ref. [17]. The scalar sector contains three-Higgs doublets and a complex scalar gauge singlet field. A flavored Peccei-Quinn symmetry provides a solution to the strong CP problem via the Peccei-Quinn mechanism, giving rise to an invisible axion which could account for the cold dark matter in the Universe. The main features of such three-Higgs flavored Peccei-Quinn (3HFPQ) class of models are summarized in Table 9, making the relevant comparisons with the KSVZ and DFSZ axion models. Experimental limits on the axion have been analyzed taking into account familon searches in rare kaon and muon decays, astrophysical considerations and axion-photon conversion experiments.  The most important findings of our analysis are:

Models
• Models with tree-level FCNCs in the down-quark or charged lepton sectors receive important constraints on the PQ scale from familon searches in kaon and muon decays. These bounds are very robust for the class of models considered in this work since the flavor changing axion couplings are completely controlled by elements of the fermion mixing matrices due to the underlying PQ symmetry.
• Models with tree-level FCNCs in the down-quark sector for which the top quark is singled out receive the strongest upper bound on the axion mass from white-dwarf cooling arguments in general, though this bounds depend strongly on the vevs of the Higgs doublets. Bounds from K + → π + a are very weak due to the strong CKM suppression. Fig. 2 summarizes all the constraints on this scenario.
• Models with tree-level FCNCs in the down-quark sector for which the up (or charm) quark is singled out receive the strongest upper bound on the axion mass from K + → π + a decays since in this case the flavor changing couplings are not as suppressed |V * us V ud | ∼ |V * cs V cd | |V * ts V td |. Fig. 3 summarizes all the constraints on this scenario.
• Constraints from µ + → e + aγ are very similar in all the models with FCNCs in the charged lepton sector due to the anarchical structure of the PMNS matrix. The bounds derived from µ + → e + aγ are stronger than those obtained from K + → π + a in models with tree-level FCNCs in the down-quark sector for which the top quark is singled out.
• The axion of models without FCNCs in the down-quark and charged lepton sectors do not receive important constraints from flavor observables. In this case the strongest bounds on the axion can be derived from the axion-photon coupling and white-dwarf cooling arguments.
• A large variety of the models considered have N DW = 1, avoiding the domain wall problem.
Allowed values for the model dependent quantity C aγ /C ag (see Eq. (85)) in these models was presented in Fig. 1, large deviations on the axion-photon coupling compared with the DFSZ model are obtained in some cases. One interesting aspect is the fact that we are able to mimic the DFSZ axion coupling to photons and have at the same time N DW = 1.
A zero C aγ can be achieved but only in models with N DW > 1.