Selective enhancement of the QCD axion couplings

We present a mechanism wherein the QCD axion coupling to nucleons, photons, or electrons, can be enhanced selectively without increasing the axion mass. We focus in particular on the axion-nucleon couplings, that are generally considered to be largely model-independent, and we show how nucleophilic axion models can be constructed. We discuss the implications of a nucleophilic axion for astrophysics, cosmology and laboratory searches. We present a model with enhanced axion couplings to nucleons and photons that can provide an excellent fit to the anomalous emission of hard X-rays recently observed from a group of nearby neutron stars, and we argue that such a scenario can be thoroughly tested in forthcoming axion-search experiments.

suddenly became one of the most serious puzzles of the standard model (SM). An elegant solution, known as the Peccei-Quinn (PQ) mechanism [4,5], was quickly put forth, and it is intriguing that, after more than four decades, it is still widely considered as the most likely explanation of why CP is a good symmetry of QCD. A striking consequence of the PQ mechanism is that an ultralight and very feebly coupled pseudo-scalar field, the axion, must exist [6,7].
In the first and simplest realisation of the PQ mechanism, the so-called Weinberg-Wilczek (WW) model [6,7], the axion couplings to SM fields were not sufficiently suppressed and the model was soon ruled out by laboratory experiments (for a historical account of early searches for the WW axion see e.g. Section 3 in Ref. [8]). Two types of models ensuring a sufficient suppression of all axion couplings were then put forth, the Kim-Shifman-Vainshtein-Zakharov (KSVZ) [9,10] and the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) model [11,12]. Although these types of axions were initially dubbed 'invisible' because of the feebleness of their couplings, Sikivie showed that search strategies exploiting the axion coupling to photons (g aγ ) could still allow to reveal these elusive particles [13]. However, three decades of experimental efforts have kept probing the axion-photon coupling without yielding to its discovery. Interestingly, recent new developments led to the possibility of searching for the axion by exploiting other couplings besides g aγ . In particular, CASPEr-Wind [14] exploits the axion-nucleon coupling to search for an axion Dark Matter (DM) wind, originating from the relative motion of the Earth with respect to the Galactic DM halo [15]. Another detection strategy is implemented by the ARIADNE collaboration [16,17], which use nuclear magnetic resonance techniques to probe an axion-mediated monopole-dipole force, sourced by a macroscopic unpolarised material and detected via a polarised sample of nucleon spins. Similar approaches involving electron spins are pursued by QUAX-g e [18] and QUAX-g p g s [19].
Presently, the sensitivity of these experiments is still far from the parameter space region of canonical QCD arXiv:2010.15846v1 [hep-ph] 29 Oct 2020 axion models. 1 It is then natural to ask whether these experiments could be already probing other types of noncanonical QCD axion models which, although they lie in parameter space regions away from the canonical benchmarks, can still provide a solution to the strong CP problem.
In this work we discuss a construction wherein the QCD axion coupling to nucleons, photons, or electrons, can be enhanced selectively without increasing the axion mass m a . 2 QCD axion models based on this construction can then populate regions in the mass-couplings parameter space which are generally believed to be accessible only to (light) axion-like particles (ALPs), although a peculiar characteristic of such QCD axions is that in most cases they are endowed with flavour-violating interactions. Our construction takes inspiration from the clockwork mechanism [20][21][22][23][24], however, it differs from the clockwork axion model discussed in Ref. [25] in that it introduces n+1 Higgs doublets and a single SM singlet complex scalar Φ, and also because, similarly to DFSZ types of scenarios, the QCD anomaly is due to the SM quarks rather than to new heavy coloured states. The construction is in fact more similar to the types of models presented in Refs. [26,27] in which Higgs doublet clockwork gears were used to obtain a 2 n enhancement of the axion-photon or axion-electron coupling.
To illustrate the main features of our mechanism we first focus on the axion coupling to nucleons (g aN , with N = p, n) which are generally considered to be the most model independent of all couplings, and we show that various modifications, and in particular large enhancements, are instead possible. Note that a first step in the direction of constructing axion models with modified axion-nucleon couplings was made in [28], where it was shown that variant axion models characterised by generation-dependent PQ charges can lead to a strong suppression of g aN . The possibility of enhancing g aN was instead considered in Ref. [29], where the value of the axion-nucleon coupling was decoupled from that of the axion mass by assigning U(1) PQ charges to SM quarks such that the latter do not contribute to the QCD anomaly. An exponential enhancement of the axionnucleon coupling is then achieved via the introduction of several complex scalars Φ k hosting in their orbital modes the axion, coupled via a clockwork-like potential. This construction, however, requires effective dimension five operators in order to eventually couple the axion to the light quarks, while the QCD anomaly of the PQ current is instead due to new KSVZ-like coloured 1 In the case of ARIADNE, some extra assumptions about the structure of CP violation are also required to yield a detectable signal. 2 By this we mean that the QCD relation between ma and the axion decay constant fa is not modified. This also implies that the axion coupling to the neutron electric dipole moment (nEDM), which depends on this relation, is also unaffected.
fermions. Here we show that a similar result can be obtained with just one SM singlet scalar Φ, and without the need of non-renormalizable interactions, by introducing additional Higgs doublets. This has also the advantage of allowing to enhance different axion couplings rather than g aN .
As regards the nucleophilic axion, the possibility of having very light axions with large couplings to the nucleons implies a rich phenomenology, and opens up a parameter space region that can be largely probed by the next generation of axion experiments. Indeed, for axion masses below the µeV, searches at ABRA-CADABRA phase 1 with resonant signal readout [30], and at CASPEr-Wind [14] will cover all scenarios with more than 15 additional doublets. In particular, ABRA-CADABRA can have sufficient sensitivity to probe models with just 5 extra doublets for a neV axion mass. For larger masses, projected sensitivities of KLASH [31], CAPP [32], and MADMAX [33] will completely cover the mass range between ∼ 1 µeV and ∼ 500 µeV. The rest of the parameter space at larger masses could be finally tested at ARIADNE [16] under the assumption of maximal CP violation. The next generation axion helioscope IAXO, will also be able to probe a wide mass range [34] significantly beyond the limit from SN1987 cooling, while some regions in parameter might be accessible already by BabyIAXO [35].
Finally, to highlight the flexibility of axion models based on our construction, we address a specific issue which is related to the observation of an excess of hard Xrays emitted from a group of nearby neutron stars (NS) referred to as the "Magnificent Seven" (M7) [36]. As was argued in Ref. [37] interpreting this excess as due to axions produced in the neutron star (NS) core and converted into photons in the NS magnetic field requires a sufficiently light axion mass (below ∼ 10 µeV) and at the same time couplings to both nucleon and photons considerably stronger than the ones predicted by canonical QCD-axion models. This is precisely the type of axion that our construction can accommodate, and we show that a nucleophilic axion can in fact provide a very good fit to the observed anomaly.
The paper is organised as follows: in Section II we recall the form of the axion interactions with the SM states by writing down the relevant effective Lagrangian, and we introduce the notations. In Section III we describe the details of the construction, focusing on the necessary ingredients to obtain a nucleophilic axion, and illustrate the reasons why one can generally expect flavour violating axion interactions. In Section IV we explore the phenomenological consequences of our scenario and the prospects for experimental probes in the next-future. In Section V we draw our conclusions. Two Appendixes complement this paper. In Appendix A we present two alternative constructions yielding respectively a photophilic and electrophilic axion. In Appendix B we discuss some issues related to the structure of the quark Yukawa matrices that can be enforced by the PQ symmetry of our clockwork-inspired multi-Higgs model.

II. AXION EFFECTIVE LAGRANGIAN
In order to set notations, let us recall the expression of the effective Lagrangian describing the axion interaction with photons and with matter fields f = p, n, e: where α is the electromagnetic (EM) coupling constant and with C a = 0.038(5) c 0 s + 0.012(5) c 0 c + 0.009(2) c 0 b + 0.0035(4) c 0 t the sea quarks contribution. In Eq. (2) E and N are respectively the EM and QCD anomaly coefficients that are defined in terms of the anomalous PQ current with α s the strong interaction coupling constant. While the axion couplings to the quarks c 0 q with q = u, d, s, c, b, t appearing in Eq. (3) and (4) are defined by the Lagrangian term Taking a Yukawa term q L q R H q , a simple expression for c 0 q in terms of PQ charges is [27] where in the last step we have replaced the fermion PQ charges X q L,R with the charge X Hq of the corresponding Higgs doublet. Finally, a common way to rewrite Eq. (1) which will be used in Section IV is where the second term is obtained by integration by parts and we have defined

III. ENHANCEMENT MECHANISM AND NUCLEOPHILIC AXIONS
Enhancing selectively a given coupling of the axion typically requires a mechanism to sizeably increase the axion coupling to the nucleons, to the electrons or to the photons, without increasing at the same time the coefficient of the QCD anomaly. In this section, we will illustrate how this can be realised in order to generate a nucleophilic axion.

A. Generation-dependent PQ charges
The key ingredient of our construction is the existence of a large hierarchy among the PQ charge differences X q L −X q R for the quarks of different generations. We will start by assuming that the charge differences for the first generation have hierarchically large values. The overall contribution of the first generation to the coefficient of the QCD anomaly, however, vanishes if the value for the up quark is equal in size but opposite in sign to that of the down quark. The anomaly coefficient then remains determined by the charges of the quarks of the other two generations, which we assume to have O(1) charge differences. Given that for each quark flavour the L-R charge differences must match the charge of the Higgs doublet responsible for the mass of that specific quark, see Eq. (8), it is clear that to realise this scenario the Higgs sector must be extended to include additional Higgs multiplets. Hence we start by assuming that the fermion Yukawa couplings involve three Higgs doublets H 0 , H 1 and H n with hypercharge Y = − 1 2 and PQ charges X 0 , X 1 and X n . For the third Higgs doublet we assume a hierarchically large charge value: where, without loss of generality, we have taken the difference between the first two PQ charges to be positive, X 1 − X 0 > 0. A detailed mechanism that can produce the charge hierarchy in Eq. (11) will be discussed in Section III C. Next, we assume that the three Higgs doublets couple to the quarks via the following generationdependent Yukawa operators: whereH = iσ 2 H * . This structure realises the conditions described above: the axion couplings to the u and d quarks are proportional to the large charge X n and they do not get particularly suppressed by the coefficient of the QCD anomaly that it is fixed in terms of the small charges of the quarks of the third generation: Consequently, also the axion couplings to the nucleons in Eqs. (3)-(4) get significantly enhanced by the large couplings of the light quarks |c 0 The Cabibbo-Kobayashi-Maskawa (CKM) mixing angles can be generated by adding to Eq. (12) intergenerational operators as for example (c L u R H 1 ) + (t L u R H n ), together with additional terms consistent with the charge assignments implied by the latter two operators plus those in Eq. (12). We anticipate that the axion will be eventually contained in the neutral 'orbital modes' of the doublets H 0,1,n and thus, since the charge assignments in Eq. (12) are generation-dependent, once the quarks are rotated to their mass basis the axion couplings to the quarks will in general be flavour-violating, see Section III B.
As regards the leptons, they couple to the complex conjugate Higgs doubletsH 0,1,n . The EM anomaly coefficient E is then readily obtained as: It is clear that depending on which specific doublet the leptons are coupled to, the EM anomaly could be also enhanced by the large charge X n , or could remain of the order of the QCD anomaly. For the nucleophilic axion we are interested in the latter possibility, so we will assume that the leptons couple universally toH 0 : The E/N factor is then given by: It is important to remark at this point that the pattern of Yukawa couplings in Eq. (12) can be straightforwardly re-arranged in a different way to yield: • An axion dominantly nucleophilic: g aN g ae , g aγ . This can be obtained for instance by coupling the up and down quarks and the τ lepton to H n andH n respectively. We obtain E, N ∼ O(1) but enhanced coupling to nucleons and to the τ lepton.
• An axion dominantly photophilic: g aγ g aN , g ae . This is easily obtained by coupling for instance only the third generation leptons toH n . We obtain E ∝ X n , N ∼ O(1) and enhanced coupling to the τ .
• An axion dominantly electrophilic: g ae g aN , g aγ . This can be obtained by coupling among the leptons only the electron toH n , and for example b and c respectively toH n and H n . In this way the 'large' contributions to E and N cancel out, and g an is not enhanced. We obtain E, N ∼ O(1) but enhanced coupling to e (and to b and c).
Two examples of nucleophilic axions will be discussed in Section III D, while two model realisations for a photophilic and an electrophilic axion will be presented in Appendix A.
The large hierarchies in the couplings described above imply that large radiative contributions to suppressed couplings are possible. For example, a large g aγ would generate, via a triangular loop, a large radiative contribution to the lepton couplings that could be relevant when the tree-level value c 0 is not particularly large [38,39] where Q 2 = 1, µ IR is the IR scale at which the coupling is evaluated and the second equality holds for E/N 1. Radiative corrections to g ae from large axion couplings to the quarks are instead a two-loop effect and are more suppressed. Similarly to the leptons, the axion coupling to the quarks would also receive a large radiative contribution from an enhanced g aγ analogous to the one in Eq. (18). The converse is instead not true: large axion-couplings to the fermions do not yield large radiative contributions to g aγ . This is because in the effective field theory limit m f → ∞, the axion-photon coupling is solely fixed in terms of the ratio of anomaly coefficients E/N and does not renormalize. The contribution of fermion loops in fact requires a helicity-flip, and for finite m f is suppressed as m 2 a /m 2 f and thus completely negligible. Finally, when the couplings of the first generation quarks c 0 u,d are not particularly large, while at least one quark of the heavier generations has a much larger coupling, the axion-nucleon coupling g aN might become dominated by the sea quark contribution C a . From the expression of C a given below Eq. (5) it can be seen that if c 0 t /c 0 u,d > ∼ 250 even the contribution of the top quark would exceed that of the valence quarks.

B. Flavour violating axion couplings
The set of Yukawa operators in Eq. (12) implies that the corresponding model belongs to the class of multi-Higgs doublet models with no natural flavour conservation [40], in which the exchange of Higgs scalars can represent a dangerous source of flavour changing neutral currents (FCNC). These effects, however, can be safely suppressed by assuming the so-called decoupling limit [41] which ensures that a single neutral Higgs scalar with properties indistinguishable from that of the SM-Higgs boson survives in the low-energy spectrum, so that the model can be rendered consistent with limits on FCNC processes and with the LHC measurements of Higgs properties.
Another source of flavour violation which does not decouple in the same limit is represented by flavour violating axion interactions, which appear when the quark fields in Eq. (12) are rotated to the mass basis. Besides axial-vector couplings analogous to the ones given in Eq. (7), off-diagonal axion-quark interactions are characterised also by vector interactions. In matrix notation, the axion-quark couplings can be written as where q are vectors containing the three types of up-or down-type quarks, and c 0V,A q (with elements c 0V,A qiqj ) are 3 × 3 matrices of couplings defined as: where X q L,R are diagonal matrices of the PQ charges of the quarks, e.g. X u L = (X u L , X c L , X t L ) T , and U q L,R are the quark unitary rotation matrices. Note that flavour violating effects will depend on the differences between the charges of the same-type (up-or down-, L or R) quarks of different generations, and get enhanced when these differences are large. The most relevant bounds on axionquark flavour-violating effects arise from mesons decays, which yield the following limits [42,43]: where for simplicity we have used for the quark mass eigenstates the same labels q = d, s, b, . . . . With the charge assignment implied by Eq. (12), we can use Eq. (20), obtaining for example the down-type quarks of the first and second generation: where we have defined x d L = (X d L − X 0 )/(X n − X 0 ) and similarly for x d R . A suppression of c V sd can be straightforwardly obtained for example by fixing U d L 12 ∼ U d R 12 ∼ 0 and generating the Cabibbo angle solely from the upquark sector, that is by assuming that the large PQ charges are associated with particularly small mixing angles. In the following, we work for simplicity in the approximation in which inter-generational mixing effects can be neglected, so that the current basis coincides, to a good approximation, with the mass basis. 3

C. Clockwork enhancement of the PQ charges
In order to generate a hierarchy in the PQ charges as given in Eq. (11) we add to the three doublets H 0,1,n an additional set of n − 2 scalars H 2 , H 3 , . . . , H n−1 , for a total of n + 1 Higgs doublets all with hypercharge Y = − 1 2 . Besides the extra doublets, we also introduce an electroweak singlet scalar field Φ with PQ charge X Φ and vacuum expectation value (VEV) v Φ v. We assume that Φ couples to H 0 and H 1 via one of the following two renormalizable terms (14)), the QCD potential has the same periodicity than the axion field, hence there is a single potential minimum and the number of domain walls (DW) [44] is N DW = 1. With the second choice 2N = 2X Φ , there are two physically distinct but degenerate minima, and N DW = 2. It should be noted that this result crucially depends on the fact that the quarks that determine the anomaly and the field Φ couple to the same pair of Higgs doublets H 0,1 , and it would not hold if, for example, Φ is coupled to a different pair of doublets or if, maintaining the scalar couplings in Eq. (23), X n contributes to the QCD anomaly.
The n+ 2 scalar fields {H i , Φ} carry a U (1) n+2 rephasing symmetry H k → e iα k H k , Φ → e iαΦ Φ. We will assume that, in addition to one of the operators in Eq. (23), the scalar potential also contains the following set of quadrilinear terms so that U (1) n+2 is broken explicitly to U (1) PQ × U (1) Y .
Since H 0,1,n need to pick-up a VEV to generate the quark masses, even in the case when all the additional doublets H k (k = 2, . . . , n − 1) have positive mass square terms, they will still acquire induced VEVs because they appear linearly in the terms in Eq. (24). This feature can be used to generate a significant hierarchy between the VEVs. Indeed, assuming that H 0,1 acquire respectively the VEVs v 0 and v 1 , and that all the other doublets have positive mass square terms µ 2 k > 0 the induced VEVs will read so that small VEVs can be typically expected if the masses of the k > 2 doublets are large µ k v 0,1 or if the couplings of the operators in Eq. (24) are small.
The arrangement of the quadrilinear Higgs couplings in Eq. (24) implies that the PQ charges X (H k ) = X k satisfy: that is X n is exponentially enhanced with respect to ∆ X . 4 In the presence of many Higgs doublets carrying PQ charges, identifying the physical axion and deriving its couplings to the fermions involves some subtleties. To ensure that the axion has no component in the longitudinal mode of the Z boson, one has to impose an orthogonality condition between the PQ and hypercharge currents where v 2 = n i=0 v 2 i ≈ 246 GeV is the electroweak breaking VEV. From Eq. (27) we see that the values of the PQ charges are determined in terms of the charge difference ∆ X and of the structure of the Higgs doublets VEVs. Let us define X 0 and X n can then be written as where the second equation makes use of Eq. (26). From these two equations we obtain which makes clear that the hierarchy between the PQ charges is determined by the parameter κ. In principle κ could range in the interval [0, 1] with the smallest values requiring a strong suppression of the VEVs of the doublets with the largest PQ charges. However, as we will see in the next Section, the phenomenologically allowed range is in fact slightly narrower. The value of κ will be crucial to determine the values of the axion couplings. Note that the anomaly coefficients E, N do not depend on κ, as can be explicitly verified by inserting in Eq. (17) X n − X 0 = 2 n−1 ∆ X . Namely, E, N are insensitive to the particular VEV structure, as could have been expected since anomalies do not depend on IR physics.
The physical axion field is defined in term of the VEVs of the Higgs doublets and of the electroweak singlet and of their neutral pseudo-scalar components a i and a Φ as Note that due to the exponential enhancement of the Higgs charges X i for large values of i, it might not be always accurate to approximate a ≈ a Φ . What remains true, is that the scale that suppresses all axion couplings is bounded from below by v a > v Φ , and hence for sufficiently large values of v Φ all current limits on the axion couplings can be easily evaded.

D. Axion couplings to matter
The axion coupling to the SM quarks and leptons depend on the particular Higgs doublet to which the fermion couples, and on the value of κ, that is on the vacuum structure. From Eq. (12) and Eq. (16) we see that the second generation quarks, the b-quark and the leptons interact with H 0 . Their couplings to the axion are then given by (see Eq. (29)): The quarks of the first generation couple to H n , so that from Eq. (30) we have We see from these equations that the parameter κ is crucial to determine a hierarchy in the couplings strength, especially if its value can approach the boundaries of the interval [0, 1]. Let us denote by v f the VEV of the Higgs doublet coupled to the fermion f . Perturbativity of the Yukawa couplings y Therefore, while κ could in principle vanish for v 0 → v, such configuration is clearly forbidden by non-zero quark masses. In particular, from Eq. (12) and (36) the perturbativity of the top Yukawa in the n + 1 Higgs theory translates into (v 1 /v) 2 (m t /(v √ 2π)) 2 ≈ 0.08. Then a lower bound on κ can be readily obtained by retaining only the first term of the sum Eq. (28) Due to the suppressing exponential factor 2 n−1 the lower value for κ does not depart significantly from zero. Note that the only contribution to κ that has no exponential suppression comes from v 2 n /v 2 , but for this contribution the perturbative limit comes from the down quark mass, and it remains below 10 −10 . As was remarked after Eq. (31), the maximum value of κ is obtained when v n is maximum, that is when v 2 n ≈ v 2 − v 2 1 . From this we obtain the upper limit where the last relation holds in the limit of large n. All in all the phenomenologically allowed range for κ is In terms of κ the axion couplings to the SM fermions and to the photon read (model A): g ae = 2 n−1 κ + δc 0 (m e /f a ) .
In order to highlight the exponentially enhanced contributions, in writing these equations we have made some approximations: in the first three relations we have neglected the model independent contributions to the couplings (the pure numbers in the left-hand side of Eqs. (2)-(4)) which are clearly subdominant, and for the axionphoton coupling we have also omitted the factor of 8/3 appearing in Eq. (15). However, in Eqs. (41) and (42) we have included the s and c sea quark contributions which can also get exponentially enhanced. Finally, let us recall that the couplings to the fermions also receive a radiative contribution from triangle loops involving g aγ . In view of the upper limit on κ in Eq. (39) this correction is irrelevant for the axion-nucleon couplings so it has been neglected in the expressions for g ap and g an , but it can become important for g ae in the limit κ → 0, so in Eq. (43) we have included the corresponding correction δc 0 which is given in Eq. (18).
For generic values of κ within the range given in Eq. (39) the structure of the couplings in Eqs. (40)- (43) naturally favours an enhancement of the axion interaction with the nucleons and with the photon. However, as was anticipated in Section III A, we can easily arrange a pattern of Higgs-fermion couplings different from the ones given in Eq. (12) and Eq. (16), and produce other types of unconventional axions, for instance dominantly photophilic or dominantly electrophilic. Since the corresponding models can also be of phenomenological interest, we discuss some examples in Appendix A. Here we discuss the possibility of generating a certain suppression of the axion couplings to the nucleons with respect to the coupling to the photon, which in Eqs. (40)- (42) have similar enhancements. This might in fact be desirable in view of the strong limit on g aN from the duration of the SN1987A neutrino burst [45] and, in particular, it is required if one attempts to fit the observed excess of hard X-ray emission from a group of nearby NS [37] in terms of axion emission from the NS core, a possibility that will be analyzed quantitatively in Section IV. A simple way to suppress to a certain extent g aN is to couple the 'large charge' Higgs H n to second generation quarks. The axion-nucleon interaction then receives the dominant contribution from the s and c sea quarks, which have additional suppression factors, and we have (model B): with in this case 10 −6 κ 0.92, where the lower limit corresponds to require that the charm Yukawa coupling remains perturbative. 5

IV. PHENOMENOLOGICAL IMPLICATIONS
The models discussed above have a rich phenomenology. The enhanced coupling to nucleons allows experimental searches through novel methods that have been put forth in recent years, and it also implies an efficient productions in stellar environments at low values of the axion mass. Both models A and B in fact predict enhanced couplings to the photons and, depending on the choice of the parameter κ, the couplings to electrons might get enhanced as well, so that these models can be called astrophilic. To satisfy the good agreement between stellar evolutionary models and observations, the axion decay constant f a should then be increased to very large values to counterbalance the effects of the large couplings induced by the exponentially enhanced PQ charges. Thus, stellar evolution forces astrophilic axions to be unusually light. For example, astrophysics bounds allow the well-studied DFSZ axion to have a mass up to about 10 meV, whereas for our nucleophilic axions, for n 15 the mass is constrained to lie below 1 µeV (see Figs. 1 and 2). The large couplings/small mass feature of these models is quite interesting from the experimental prospective since, as will be discussed in Section IV B, several proposed experiments can have enough sensitivity to probe large portions of the low mass region not yet excluded by astrophysics observations.

A. Astrophysics
Stellar evolutionary theoretical studies, combined with accurate observations of stellar populations, lead to strong constraints on the axion couplings to matter and radiation (see, e.g., [27,[46][47][48][49]). The bounds from stellar evolution for model A, corresponding to the couplings in Eqs. (40)- (43), are shown in Fig. 1 and 2. In these plots the x-axis corresponds to the axion mass, and the y-axis to the number of Higgs doublets (i.e. the number of clockwork gears) that in our construction are responsible for the enhancement of the couplings. In Fig. 1 the parameter κ has been fixed to its minimum value κ = 0 which implies that g ae is determined solely by the radiative contributions and hence strongly suppressed. Fig. 2 corresponds to to the maximum value allowed by perturbative unitarity on the charm Yukawa coupling κ = 0.92, in which case g ae is exponentially enhanced.
The hatched grey area in the figures represents the region of parameters excluded by the duration of the SN1987A neutrino signal. Historically, observation of the SN1987A neutrinos has provided the strongest bounds on the axion-nucleon couplings [45,[50][51][52]. In the plots we use a state-of-the-art determination of the SN1987A limit from Ref. [45]. The lower edge of the grey re- gion corresponds to axions so weakly coupled that, while they can freely stream out of the supernova, the amount of energy they can carry away does not shorten sufficiently the neutrino burst duration. In this regime, the limit applies to the following combination of couplings g 2 an + 0.6g 2 ap + 0.5g an g ap . The upper edge is determined by the onset of the trapping regime, in which axion interactions are sufficiently strong that their mean free path is smaller than the proto-NS radius so that, like neutrinos, axions remain trapped inside the star for a sufficiently long time not to affect significantly the neutrino burst duration [45,51]. Both bounds are afflicted by several uncertainties and should be taken with a grain of salt. However, the upper edge of the region is already covered by the CAST results (see Section IV B) so that determining the precise values of the couplings for the onset of the trapping regime is not crucial. In contrast, a reliable assessment of the validity of the limit for the free streaming regime is an important issue. In the pictures we use the bound derived in [45]. However, it was recently claimed that axion-pion interactions may contribute more than previously thought to the axion emissivity, in which case the limit could be sizeably stronger [53].
At low values of the axion mass, the region in which axions can account for the M7 anomaly is in strong tension with several astrophysical observations. These include the recent NuSTAR search for hard X-rays emission from Betelgeuse [54] and from the Quintuplet and West-erlund 1 super star clusters [55]. Both analyses exclude the region m a 10 −10 eV for g aγ > ∼ few×10 −12 GeV −1 . A similar region is also excluded by the Fermi LAT bound on diffuse gamma rays that would result from conversion of SN axions into photons in the Galactic magnetic field [56]. At even smaller masses, m a 10 −11 eV, CHANDRA observations of NGC1275 set a slightly stronger bound on the axion-photon coupling [57]. The strongest constraint is, however, implied by the nonobservation of gamma rays by the Solar Maximum Mission (SMM) at the time of SN1987A explosion [58][59][60], and corresponds to the purple hatched region in the figures. We see from the pictures that this bound is especially relevant for n > ∼ 20, and this is because the enhancement of the axion-nucleon couplings in this regime strongly enhances the SN axion emissivity.
The orange line labelled HB in the figures represents the horizontal branch (HB) star bound [61,62], which constraints the axion-photon coupling from the ratio of the observed number of stars in the HB and RGB in globular clusters (the bound shown in the figures is from the latest analysis in Ref. [62]).
The strongest astrophysical bound on the axionelectron coupling g ae is derived from observations of red giant branch (RGB) stars [63][64][65]. The RGB bound in Fig. 1 is very weak since in model A, κ ≈ 0 implies g ae ≈ 0 at tree level. However, as can be seen from Fig. 2, g ae gets exponentially enhanced for κ = 0.92, and then the RGB bound dominates over all other limits, including that from SN1987A. (The same can occur also in other constructions, as for example in model B at large κ. In this specific case, however, g aN also gets enhanced, and the RGB bound is only slightly more constraining than the SN1987A limit.) Finally, the yellow area in the figures corresponds to the hint from a 2 σ fit to the observed excess of hard X-ray events from the Magnificent Seven (M7) group of NS [37]. The origin of this anomaly is not understood, and it was speculated in Ref. [37] that the excess might be attributed to axion-like particles (ALPs) with couplings to both photons and neutrons. The axion-neutron coupling would be responsible for the production of ALPs in the hot NS core, while the coupling to photons would allow the ALP to be converted into photons in the strong magnetic fields surrounding the NS. The resulting photon flux would then be detected by X-ray detectors, such as XMM-Newton and Chandra, and would correspond to the excesses observed by these instruments. Notice that the observations demand quite large couplings and yet a small mass, in order for the axions to be efficiently converted into photons in the magnetic field. More specifically, Ref. [37] found that the mass should not exceed ∼ 20 µeV, while the couplings should satisfy the relation g aγ g an ∼ a few 10 −21 GeV −1 . Given the upper limit on the mass, such couplings are prohibitively large for canonical QCD axion models, such as DFSZ (see the red vertical segment in Fig. 3), and that is why one would The gan-gaγ parameter space for ma = 10 µeV. The vertical red segment corresponds to the DFSZ axion with gan within the phenomenologically allowed range. The blue strip corresponds to model B for 0 κ 0.92. The grey hatched areas are excluded respectively by the SN1987A and HB bounds. The (2σ) region for which the M7 anomaly can be explained in terms of axion/ALP emission/conversion [37] corresponds to the yellow band. more generically invoke an ALP. However, our construction is versatile enough to provide a QCD axion that fits well the data even in this extreme case. In fact, thanks to the exponential enhancement of the axion couplings to nucleons and photons, axion production and their conversion into photons can proceed with sufficient rates even in the small mass window. This is shown in Fig. 3, where we present the axion parameter space for m a = 10 µeV. The red line refers to the DFSZ axion, constrained by perturbative unitarity. The blue band represents the parameter space for model B with κ varying in the allowed range 0 κ 0.92 and 2 ≤ n ≤ 40. We see that the photophilic and nucleophilic axion of model B can have sufficiently large couplings to account for the M7 anomaly even at this low mass value. As regards model A, the combination of the SN1987A and RGB bounds always excludes the region of couplings hinted by the M7 anomaly. In fact, in order to evade the SN1987A bound, model A requires a large κ which, in turn, is strongly constrained by the RGB bound.

B. Experimental bounds and perspectives
In this section we discuss the potential of current and next generation axion experiments to probe the nucleophilic models. For definiteness, we show the results for model B, though many of our conclusions apply to In green are the cavity experiments, with ADMX in darker green. The other green areas comprise the projected sensitivities of KLASH [31], CAPP [32], and MADMAX [33]. The region in light purple indicates the reach of ABRACADABRA phase 1 with resonant signal readout [30]. The sensitivity of CASPEr-Wind phase 1 and 2 to gaN [14] corresponds to the two red lines. The reach of ARIADNE is that enclosed in the brown line, at relatively large masses [16], and the projected sensitivities of IAXO [34] and BabyIAXO [35] are given by the two blue lines.
model A as well.
Current axion searches are probing (mostly) the axion coupling to photons. However, some proposals for future search strategies suggest exploiting also the couplings to nucleons, and in particular to neutrons (see [13,66] for comprehensive reviews). Presently, one of the tightest experimental bounds come from the CERN Axion Solar Telescope (CAST), which has probed the axion-photon coupling down to g aγ = 0.66 × 10 −10 GeV −1 for masses up to 0.02 eV [67].
In Fig. 4 we show the viable and the excluded regions for model B for κ ≈ 0. The region excluded by CAST is shown in blue. The region probed by the cavity haloscope ADMX [68,69], which covers so far only a narrow mass window, is shown in dark green. Next generation cavity experiments are expected to probe a considerable wider mass range from ∼ 300 µeV to ∼ 0.5 meV, as shown by the light green region. The region includes the KLASH [31], CAPP [32], and MADMAX [33] propos-als. The region at low masses can be probed by ABRA-CADABRA [30] and CASPEr-Wind [14]. As shown in the figure, these experiments are expected to cover a large portion of the sub µeV region of the nucleophilic axion model. Finally, the ARIADNE proposal [16] would allow to probe the large mass region through the axion coupling to neutrons. However, the corresponding limits will hold only under the assumption that the amount of CP violation is maximal and saturates the neutron electric dipole moment bound.
The next generation of solar axion searches will have a considerable higher sensitivity than CAST, and will be able to set limits for a wide range of axion masses. The International Axion Observatory (IAXO) [34], will have enough sensitivity to test the axion explanation of the M7 anomaly anomaly. BabyIAXO [35], an intermediate experimental stage of IAXO which is expected to become operative by the mid of the current decade, will already expand considerably the region probed by CAST, and can already probe a large part of the M7 region. Finally, the forthcoming light-shining-through-walls experiment ALPS II [70] (not shown in Fig. 4), which is expected to take data starting from 2021, will also surpass CAST, probing the parameter space almost to the level of sensitivity of BabyIAXO, although in a smaller mass window m a 0.1 meV.
Let us note that, with the exception of CASPEr-Wind and ARIADNE, all the axion experiments included in Fig. 4 probe the axion-photon coupling, which has the same form in model A than in model-B, and in particular it does not depend on κ. Hence, similar results can be expected for model A and for other choices of κ.

C. Cosmology
Exponentially enhanced axion couplings may lead in the primordial Universe to a thermal population of hot axions, and this would modify the effective number of neutrinos with respect to that inferred by the Planck collaboration from CMB measurements [71]. Given that the axion couples strongly to first generation quarks, we will first consider the regime in which axions decouple from the thermal bath at T d 1 GeV. In this regime quarks are bounded into nucleons whose number density is Boltzmann suppressed, and into pions that (for T dec > ∼ m π ) are as abundant as photons, so that axion coupling to pions is the relevant quantity. The axion thermal production rate Γ aπ has been estimated in [72]: where f π = 93 MeV, h is an exponentially decreasing function satisfying with h(0) = 1 and is tabulated in [72], and the axion-pion coupling is: Taking as an example model A, and neglecting the modelindependent contribution (the third term in parenthesis) we have: C aπ ≈ −2c 0 u /3. Using this coupling together with Eq. (48), and recalling that the decoupling temperature is defined by the condition Γ aπ (T d ) The first numerical bound corresponds to the 2σ measurement of N eff from the Planck collaboration [71], while the lower limit on the axion decay constant for model A (f A a ) has been computed in the following way: the Planck upper bound on ∆N eff implies g eff (T d ) > ∼ 15.3 which corresponds to a limit on the decoupling temperature T d > ∼ 66 MeV (see Fig.1 in Ref. [73]). For this temperature the Hubble parameter is H(T d ) 2.4 · 10 −18 MeV, and the limit is then obtained from the out-of-equilibrium condition Eq. (48) using C A aπ = −2c 0 u /3. Note that this limit is sub-dominant compared to the SN1987A bound on the nucleon couplings.
In the case where couplings to second and third generation are enhanced, thermalisation can occur via the qg → qa and qq → ga processes [74]. However, the predicted value of ∆N eff will only be in reach of future CMB-S4 experiments, leaving to a distant future the possibility of deriving constraints on the axion decay constant by using these processes.

V. CONCLUSIONS
We have presented in this work a simple structure to obtain a nucleophilic QCD axion. The key ingredient of this construction is the presence of flavourdependent Yukawa interactions of Higgs doublets with the SM quarks. In particular, one Higgs doublet with a very large PQ charge must interact with both the up and down quarks, leading to a cancellation of its contribution to the QCD anomaly and, at the same time, a large axion-nucleon coupling. We have further constructed an explicit realization of this scenario via a clockwork mechanism directly at the level of the Higgs doublets.
Interestingly, we have shown that this construction can in fact generate an exponential enhancements of almost any axion coupling (barring the nEDM one), thus realising a scenario in which the parameter region of QCD axion models can be extended to overlap with mass/couplings regions that are generally considered viable only in specific ALP models. The two main restrictions of this setup is that one can only enhance the axion couplings with respect to the standard QCD axion cases, and that the model typically predicts flavour-dependent (as well as flavour violating) interactions. Additionally, due to the gradation of the doublets PQ charges in steps of 2 i the model provides strong constraints on the allowed Yukawa couplings, similar to those found in simpler two Higgs doublet (2HD) models (see Appendix B). It would be interesting to examine in more details to which point our clockwork inspired multi-Higgs doublet model departs from the results found in these 2HD setups.
We have then analysed the phenomenology of two such nucleophilic models. First we considered a simple setup (model A), where the axion couplings to first generation quarks are strongly enhanced, along with the axionphoton interaction. We then studied a variation (model B) where the couplings to second generation quarks are instead boosted. The nucleon-axion couplings are then mostly generated by the axion interaction with the sea quarks in the nucleon. We emphasise that both scenarios can be easily tested in the near future by a large number of experiments, such as CASPEr-Wind, ABRA-CADABRA, or ARIADNE. In particular, the proposed model B provides an elegant solution for the excess of X-ray events originating from the "Magnificent Seven", a group of isolated NS, thanks to the light mass of the axion and its enhanced couplings to both photons and nucleons. We have further shown that such explanation will be probed by various upcoming axion experiments in the near future. X 22 = −X 0 and X 33 = −X 1 , and choosing Y 23 to be non-vanishing (e.g. to generate a particular CKM mixing) means that X 23 must match the charge of one Higgs doublet, e.g. X 23 = −X j . Y 32 would then be a texture zero unless X 32 also matches the charge of some Higgs H k , that is: Thus (X j , X k ) = (X 0 , X 1 ) or (X 1 , X 0 ) allows for both Y 23 , Y 32 = 0. Only two Higgs doublets are involved, and in fact this is the same solution one has in the 2HD case. For j, k ≥ 1, we see from the structure of the charges given in Eq. (26) that the coefficient of ∆ X can never be matched since 2 j−1 + 2 k−1 = 1 has no solution. Hence in the clockwork inspired multi-Higgs doublet scenario there are no additional possibilities with respect to the 2HD case to avoid one zero texture in the (23) Yukawa submatrix. Repeating this argument for the (13) subma-trix one would obtain 2 j−1 + 2 k−1 = 2 n−1 + 1 which has again the same solutions (j, k) = (1, n) or (n, 1) in terms of just two Higgs doublets than the 2HD case. Only for the (12) submatrix, for which the j, k ≥ 1 condition reads 2 j−1 + 2 k−1 = 2 n−1 , the solution j = k = n − 1 opens up a new possibility for avoiding zero textures in terms of three Higgs doublets {H 0 , H n−1 , H n }. Clearly, the same considerations will also hold for the other possible PQ charge assignments leading to model B, C and D.