Axion couplings to electroweak gauge bosons

We determine the model-independent component of the couplings of axions to electroweak gauge bosons, induced by the minimal coupling to QCD inherent to solving the strong CP problem. The case of the invisible QCD axion is developed first, and the impact on $W$ and $Z$ axion couplings is discussed. The analysis is extended next to the generic framework of heavy true axions and low axion scales, corresponding to scenarios with enlarged confining sector. The mass dependence of the coupling of heavy axions to photons, $W$ and $Z$ bosons is determined. Furthermore, we perform a two-coupling-at-a-time phenomenological study where the gluonic coupling together with individual gauge boson couplings are considered. In this way, the regions excluded by experimental data for the axion-$WW$, axion-$ZZ$ and axion-$Z\gamma$ couplings are determined and analyzed together with the usual photonic ones. The phenomenological results apply as well to ALPs which have anomalous couplings to both QCD and the electroweak bosons.


Introduction
The term "axion" denotes any (pseudo)Nambu-Goldstone boson (pNGB) of a global chiral U (1) symmetry which is exact at the classical level but is explicitly broken only by anomalous couplings to the field strength of a confining group. 1 Axions are the characteristic byproduct of solutions to the strong CP problem based on an anomalous U (1) axial symmetry, usually called Peccei-Quinn (PQ) symmetry [1]. Whatever its mass, an axion necessarily has anomalous gluonic couplings and it only deserves that name if its presence solves the strong CP problem. When the number of axions in a given theory outnumbers the total number of distinct instanton-induced scales other than QCD, one (or more) light axions remain.
In scenarios with an axion a coupled only to QCD, two states of the spectrum have anomalous couplings, the axion and the η . Only one eigenstate can then acquire a mass of the order Λ QCD : this is identified as the physical η , while the axion would remain massless in the limit of vanishing quark masses. The latter are relevant, though, and their impact on the mixing of pseudoscalars results in an axion mass given by [2,3] where m π , f π , m u , m d denote the pion mass and decay constant, and the up and down quark masses, respectively. The axion scale to which all axion couplings are inversely proportional is denoted by f a . The expression in Eq. (1.1) is characteristic of all usual "invisible" axion constructions, characterised by QCD being the only confining gauge group. Paradigmatic examples are the DFSZ and the KSVZ models [4][5][6][7] and their variations. These invisible axions are also called QCD axions. The main phenomenological constraints for QCD axions are obtained from their couplings to photons, gluons and fermions. For instance, the effective coupling to photons g aγγ is defined as δL ⊃ 1 4 g aγγ a F µνF µν , (1.2) where F µν is he electromagnetic field strength, α em the corresponding fine structure constant, and E and N denote respectively the electromagnetic and color anomaly coefficients. The ratio E/N is model-dependent and has been recently computed in Refs. [8,9] for different representations of exotic fermions in standard "invisible axion" models. This ratio contains all direct loop contributions of PQ and electrically charged fermions, including the up and down and other SM quarks if PQ charged. The second term in the parenthesis in Eq. (1. 3) is instead model-independent. [3,[10][11][12][13] It stems from the mixing of the axion with the η pseudoscalar -as both couple to G µνG µν with G µν being the QCD field strength-and also with the pion through the η -pion mixing due to the up-down quarks mass difference. Obviously, such model-independent mixings are relevant for axions lighter than the QCD confining scale, while if the axion is much heavier this effect should become negligible.
Light enough axions (that is, below O(100 MeV)) can participate in astrophysical phenomena. [14][15][16] The constraints that follow from their non-observation in photonic processes lead to very high values for the decay constant, f a > 10 8 GeV. It follows then from Eq. (1.1) that m a ≤ 10 −5 eV for QCD axions.
Nevertheless, in specific QCD axion models the coupling to photons may be suppressed [9,17] and moreover large uncertainties may hover over the purely hadronic constraints. [18] It is thus important to analyze the axion couplings to other electroweak gauge bosons, as they may become the phenomenologically dominant couplings in certain regions of the parameter space for those models. Couplings of axions and also of axion-like particles (ALPs) to heavy gauge bosons are increasingly explored [19][20][21][22][23][24] in view of present and future collider data, and also in view of rare meson decay data. For instance, in addition to LHC-related signals, recent works [25,26] consider the one-loop impact of aW W couplings on rare meson decays.
Here we first determine the model-independent components of the mixing of QCD axions with electroweak gauge bosons, which result from the mixing of the axion with the pseudoscalar mesons of the SM. In other words, we determine the equivalent of the 1.92 factor in the photonic coupling in Eq. (1.3), for the couplings of the axion to W and Z gauge bosons. A chiral Lagrangian formulation will be used for this purpose, determining the leading-order effects. The heavy electroweak gauge bosons will be introduced in that Lagrangian as external -classical-sources. Our results should impact the analyses for light axions of theories which solve the strong CP problem. They are novel and relevant in particular whenever the axion is lighter than the QCD confining scale and is onshell in either low-energy or high-energy experiments. They also impact the comparison between the data taken at experiments at low and high-energy. For instance, a null result in NA62 data for K → πa does not imply the absence of a signal at high energy in an accelerator such as that from off-shell axions at LEP or at a collider. This is because model-independent contributions are present at the low momenta dominant in rare decays (and cancellations may then take place), while at high energies they are absent.
As a second step, we will extend the analysis to heavy axions which solve the strong CP problem. Axions much heavier than Λ QCD and with low axion scales are possible within dynamical solutions to the strong CP problem, at the expense of enlarging the confining sector of the Standard Model (SM) beyond QCD. [27][28][29][30][31][32][33][34][35][36][37][38] These theories introduce a second and large instanton-induced scale Λ Λ QCD to which the axion also exhibits anomalous couplings, resulting in the bulk of its large mass. Precisely because the axion mass typically lies well above the MeV regime, these heavy axions avoid the stringent astrophysical and laboratory constraints and present and future colliders may discover them. The transition between the light and heavy axion regime will be explored for the coupling of the axion to the photon and to the electroweak gauge bosons.
Finally, the phenomenological part of the analysis will be carried out on a "two-coupling-ata-time" basis: it will take into account the simultaneous presence of a given electroweak coupling and the axion-gluon-gluon anomalous coupling (essential to solve the strong CP problem). For the analysis of data, for the first time the experimentally excluded areas for the EW couplings g aW W , g aZZ and g aγZ will be identified and depicted separately, besides the customary ones for the g aγγ coupling. Furthermore, the relations among the exclusion regions stemming from electroweak gauge invariance will be determined and exploited. Model predictions will be illustrated over the experimental parameter space.
Aside from the main focus of the paper on true axions, our analysis applies to and calls for a timely extension of the ALP parameter space. Very interesting bounds on ALPs from LEP and LHC [19-24, 39, 40] assume often just one electroweak coupling for the axion, and no gluonic coupling. The path to consider any two (or more) couplings at a time will change the experimental perspective on ALPs.
What is the difference between a heavy axion and an ALP with both anomalous electroweak and gluonic couplings? The key distinction is that the former stems from a solution to the strong CP problem while a "gluonic ALP" may not. Both exhibit anomalous couplings to QCD and in both cases there is an external source of mass besides that induced by QCD instantons and mixing. However, for a true heavy axion that extra source of mass does not induce a shift of the θ parameter outside the CP conserving minimum (and thus the solution to the strong CP problem is preserved), while for a generic gluonic ALP such a shift may be induced. Nevertheless, this important distinction is not directly relevant for this work, as the novel aspects and phenomenological analysis developed below are valid for both true heavy axions which solve the strong CP problem and for gluonic ALPs. To sum up, all results below for heavy axions apply directly to gluonic ALPs as well. In addition, the conclusions based purely on EW gauge invariance have an even larger reach: they hold for all type of axions and for generic ALPs (that is, ALPs with or without gluonic couplings).
The structure of the paper can be easily inferred from the Table of Contents.

The Lagrangian for the QCD axion
Without loss of generality, the axion couplings can be encoded in a model-independent way in an effective Lagrangian. At leading order in inverse powers of the scale f PQ at which the global PQ symmetry is broken, and at energies above the electroweak (EW) scale, it reads whereâ denotes the axion eigenstate at energies well above the confinement scale and L SM is the SM Lagrangian, In this expression, W µν and B µν are respectively the SU (2) L and U (1) Y gauge field strengths, while Φ, Q L , L L , U R , D R and E R respectively encode the Higgs doublet, the fermion and lepton doublets, the vectors of right-handed up-type, down-type quarks and right-handed charged leptons (neutrino masses will be disregarded), with ψ running over all those type of fermions andΦ ≡ iσ 2 Φ * . The three 3 × 3 matrices in flavour space Y D , Y U and Y E encode the Yukawa couplings for down quarks, up quarks and charged leptons, respectively. For the effective axion couplings, we first focus on the anomalous couplings of the axion and their impact on the pseudoscalar mass matrix. We will work in the basis in which the only PQ-breaking operators in the Lagrangian, δLP Q a , are the anomalous couplings of axions to gauge bosons. This choice is allowed by the reparametrization invariance of the effective axion Lagrangian (see Appendix C). The CP-conserving and PQviolating next-to-leading order (NLO) corrections due to axion physics then read 2 with OG, OW and OB denoting the anomalous axion couplings to gluons, SU (2) L and U (1) Y gauge bosons, respectively, where α s , α W and α B denote respectively the fine structure constants for the QCD, SU (2) L and U (1) gauge interactions, and N 0 , P 0 and L 0 are dimensionless operator coefficients. Customarily, the Lagrangian in Eq. (2.3) is rewritten as where the Lorentz indices of the field strengths are implicit from now on and The model-dependent group theoretical factors can be generically written in terms of the fermionic PQ charges X i as where X i is the difference between the right-handed and left-handed PQ charges: 3 ) are respectively the Dynkin indices of the fermionic representation R i under QCD and SU (2) L , and Y i denotes the hypercharge. The sum over "heavy" fermions and the subscript 0 indicate that the contribution to the anomalous couplings from all exotic heavy quarks, and/or heavy SM quarks (s, c, b, and t quarks) if PQ charged, is encoded in the N 0 , L 0 and P 0 operator coefficients. That is, the possible contribution of the SM first generation up (u) and down (d) quarks is not included in those coefficients. Indeed, for models in which they have PQ charges an additional PQ-invariant term must be considered, replacing the u and d Yukawa couplings in Eq. (2.2) by: which assumes as convention that the axion transforms under the PQ symmetry as a −→ a + f PQ . The X u,d dependence in Eq. (2.11) will be shown below to generate extra contributions to the physical anomalous couplings. In this equation Q 1 denotes the first family doublet, and flavourmixing effects as well as leptonic couplings are omitted from now on for simplicity. In all equations above, color and SU (2) L indices are implicit and the QCD θ angle has been removed from the Lagrangian via the PQ symmetry.
Among the most general set of purely derivative operators, additional couplings could also be considered, e.g.
where c Q 1 , c U 1 and c D 1 are matrices of arbitrary coefficients in flavour space. Nevertheless, the reparametrization invariance of the Lagrangian [43] allows to work in a basis in which these terms (which would seed pseudoscalar kinetic mixing) are absent and their impact is transferred to other axionic couplings. 4 From now on they will be disregarded all through the analysis on pseudoscalar mixing. In summary, the Lagrangian to be analyzed below when considering mixing effects reads (2.13) Below EW symmetry breaking and above confinement After electroweak symmetry breaking but before QCD confinement, the effective Lagrangian in Eq. (2.13) results in where In these equations s w and c w denote the sine and cosine of the Weinberg mixing angle and α em = α W c 2 w = α B s 2 w . For models in which the the first generation of SM quarks are not PQ charged, X u,d = 0. When those quarks are instead charged under PQ, their contribution to the anomalous gauge couplings has to be included in the group theory factors, which are replaced by

The Lagrangian below the QCD confinement scale
Three pseudoscalars mix once quarks are confined: the axion, the SM singlet η 0 and the neutral pion π 3 . The π 3 -η 0 mixing is due to the quark mass differences which break the global flavour symmetry. At leading order in the chiral expansion and in the two quark approximation, the mass Lagrangian for the pions and η 0 reads where B 0 can be expressed in terms of the QCD quark condensate qq as B 0 f 2 π = −2 qq , and M q denotes the quark mass matrix, (2.25) The matrix of pseudoscalar fields can be written as where the η 0 decay constant has been approximated by f π and In the presence of the axion, the anomalous QCD current GG couples to both the axion and the η 0 fields and mixes them. The two mixing sources combined result ultimately in an axion-pion mixing. For simplicity, we will first consider the case with X u,d = 0 in Eq. (2.11), and afterwards the case X u,d = 0.

SM light quarks not charged under PQ (X u,d = 0)
A popular example of this class of models are KSVZ ones, in which only heavy exotic quarks are charged under PQ. For any model in which the u and d SM quarks are singlets of the PQ symmetry, their quark mass matrix in the basis here considered is that in Eq. (2.25). In this case N = N 0 , as all contributions to the anomalous gluonic coupling are already included in N 0 . The potential for the three pseudoscalars is, at first order in the pseudoscalar masses, where the last term stems from the instanton potential with K ∼ Λ 4 QCD . [44][45][46] The resulting mass matrix for the three neutral pseudoscalars is given by The diagonalization leads to the well-known expressions for the pseudoscalar mass terms: (2.30) It follows from these expressions that K can be expressed in terms of the low-energy physical parameters as The corresponding mixing matrix is given by 32) or, equivalently, the mass eigenstates are given by a â + θ aπ π 3 + θ aη η 0 , (2.33) π 0 π 3 + θ π a a + θ πη η 0 , (2.34) Here, all the mixing angles are assumed small and Only the leading terms for each mixing entry have been kept in Eqs. (2.32)-(2.37). The complete matrix to first order in 1/f a and in quark masses -that is, O(m 2 π /m 2 η )-can be found in Appendix A: the impact of the extra terms may be comparable to that of next-to-leading operators in the chiral expansion and will thus not be retained here.
The results in Eqs. (2.33)-(2.37) illustrate that the physical low-energy axion eigenstate acquires π 3 and η 0 components and thus inherits their couplings to all gauge bosons, weighted down by their mixing with the axion. These results apply to any physical process in which the axion is on-shell and the axion mass is lighter than the confinement scale.
We are interested in identifying the model-independent contributions in the coupling to the electroweak gauge bosons for light axions and for the SM light pseudoscalars. We will first recover in our basis the customary axion-photon couplings, to set the framework.

Axion-photon coupling
The physical g aγγ is given by where the last two terms are the contributions induced by the model-independent axion-pion and axion-η QCD mixing. Denoting by q u and q d the electric charges of the up and down quarks, respectively, the photonic couplings of the SM light pseudoscalars are given by For the present case with X u,d = 0, that is E 0 = E and N 0 = N , it follows that resulting in the well-known expression which is valid to first order in chiral perturbation theory.

SM light quarks charged under PQ (X u,d = 0)
The quark mass matrix in Eq. (2.24) is to be replaced by The potential in Eq. (2.28) is now generalized to resulting in a pseudoscalar squared mass matrix M 2 {π 3 , η 0 ,a} which takes the form (2.44) The mixing angles in Eqs. (2.36)-(2.37) still hold but for the pion-axion mixing which is now given by 5 The coefficient in front of the mass-dependent term in this equation coincides with the strength of the physical gluonic couplings, 6 given by and in consequence The expressions for the mass of the physical pion, η and axion are the same than those in Eq. (2.30).

Axion-photon coupling
For the case in which the up and down quarks are charged under PQ, X u,d = 0, N is given by Eq. (2.46) resulting in In consequence In summary, the most general mass matrix leads to the same expression than for the case X u,d = 0 in Eq. (2.55) if taking into account in E also the contribution of the up and down quarks.

Axion couplings to EW gauge bosons
The description in terms of the effective chiral Lagrangian is only appropriate for energies/momenta not higher than the cutoff of the effective theory, 4πf π , that is, the QCD scale as set by the nucleon mass. In this context, the W and Z bosons can be considered as external currents that couple to a QCD axion whose energy/momentum is not higher than Λ QCD , for instance a light enough on-shell QCD axion. In other words, the W and Z bosons enter the effective chiral Lagrangian as classical sources, alike to the treatment of baryons in the effective chiral Lagrangian.
We extend here the results of the previous section to the interactions of axions with SM heavy gauge bosons. Eq. (2.55) is thus generalized for any pair of electroweak gauge bosons X, Y , Indeed, all axion couplings to EW bosons receive a model-independent component due to the pionη -axion mixings, in the regime in which the involved energy/momenta are smaller or comparable to the confinement scale. A relevant point when computing the couplings of the QCD axion to heavy EW bosons is the fact that, after confinement, a new type of SU (2) L -breaking effective interaction of the formâ is present in addition to those in Eqs. (2.4)-(2.6). It stems via axion-pion coupling from the π a W a µνB µν effective interaction. The details of the computation can be found in App. B. In terms of the physical photon, Z and W , the interaction Lagrangian for the QCD axion is then given by Eq. (2.55) is the known leading-order result [10,13] for the photonic couplings of the QCD axion, which holds in all generality for on-shell axions lighter than the QCD confinement scale. The contributions in Eqs. (2.56)-(2.58) are new and extend that result to the couplings of heavy gauge bosons in the appropriate energy range. Indeed, the last term in the parenthesis for each of these couplings encodes the impact of the mixing of the axion with the pion and η , see Eqs. (2.15)-(2.19) for comparison with the unmixed case. These are model-independent contributions, valid for any QCD axion, i.e., for any model in which the SM strong gauge group is the only confining force and thus the only source of an instanton potential for the axion. In other words, they are valid for any axion whose mass and scale are related by Eq. (2.30). Note that those corrections hold precisely because the axion mass is smaller than the confining scale of QCD, which is the regime in which the SM pseudoscalars lighter than the QCD scale are the physical eigenstates of the spectrum and the mixing effects are meaningful. Numerically, at leading order in the chiral expansion it follows that The numerical value of the model-independent term in Eq. (2.59) differs from the usual one [13] of 1.92 in Eq. (1.3), as the latter includes higher order chiral corrections, a refinement out of the scope of this paper and left for future work.
The model-independent results obtained here for the coupling of light QCD axions to the SM electroweak bosons may impact axion signals in rare decays in which they participate. For instance, in low-energy processes the axion could be be photophobic at low energies [17] (or more generally, EW-phobic), in models in which the terms in parenthesis cancel approximately, unlike at higher energies at which the model-independent component disappears and only the modeldependence (encoded in E/N , M/N , Z/N and R/N ) is at play.

Gauge invariance
As it was already enforced in Eq. (2.22), the couplings of the axion to the EW gauge bosons are not independent as a consequence of gauge invariance. Indeed, all four couplings stem from the two independent effective operators in Eqs. (2.5) and (2.6), plus that in Eq. (2.52) for a light QCD axion (m a Λ QCD ). In consequence, three physical couplings can be independent among the set {g aγγ , g aW W , g aZZ , g aγZ }, and the following relation must hold: Note that this result does not depend on the details of any particular axion model; it is independent of the presence or absence of gluonic couplings and it thus applies in all generality for a light pseudoscalar with only anomalous electroweak couplings. That is, it is also valid for generic ALPs which only have EW interactions. Furthermore, it suggests that it may be inconsistent to assume only one EW coupling: the minimum number of physical EW couplings for an axion or ALP is two. The relation in Eq. (2.63) sets an avenue to oversconstrain the parameter space which is promising. It allows to use the better constrained EW couplings to bound the fourth one.

Beyond the QCD axion
We discuss in this section the case of a "heavy axion": an axion whose mass is not given by the QCD axion expression in Eq. (1.1) but receives instead extra contributions. We have in mind a true axion which solves the strong CP problem, for which the source of this extra mass does not spoil the alignment of the CP conserving minimum. This is the case for instance of models in which the confining sector of the SM is enlarged involving a new force with a confining scale much larger than the QCD one [27][28][29][30][31][32][33][34][35][36][37][38]. This avenue is of particular interest as it allows to consider heavy axions and low axion scales (e.g. O(TeV)), and still solve the SM strong CP problem. The axion mass can then expand a very large range of values. It can become much larger than the EW scale or, conversely, be in the GeV range or lower. For the purpose of this work, the latter range is to be kept in mind as a general guideline, so as to remain in the range of validity of the effective Lagragian with confined hadrons. The procedure will also serve as a template to show how the mixing effects disappear from the axion-gauge couplings as the axion mass is raised.
In practice, the analysis below applies identically to a true heavy axion which solves the strong CP problem and to a gluonic ALP, that is, any ALP which has both electroweak and gluonic anomalous couplings, even if not related to a solution to the strong CP problem. For the sake of generality, consider the addition of an extra mass term to the effective Lagrangian obtained after EW symmetry breaking but above confinement in Eq. (2.14), no PQ charge, as it is straightforward to enlarge the analysis beyond this hypothesis, as shown in the previous section. After confinement, the pseudoscalar mass matrix in Eq. (2.29) is then replaced by In the limit f a K 1/4 (i.e. f a Λ QCD ) the expressions for the π, η and a mass eigenvalues become The corresponding axion-pion and axion-η mixing angles take a very simple form in the limit M m π , where again only the leading term on each entry of the mixing matrix has been retained. In fact, it can be checked that these equations hold even for small values of M , as long as M is non-degenerate with the pion or η mass. The comparison with Eq. (2.36) ilustrates that the M-dependent corrections may be important for axion masses near the pion mass or the η mass.
The divergences in Eq. (3.4) for an axion degenerate in mass with the pion or the η are an artifact of the approximations which in practice correspond to large mixing values, as expected in those ranges. The numerical results do not rely on that approximation and are illustrated in Fig. 1: the spikes correspond to an axion degenerate with either the pion or the η .

Heavy axion couplings to EW gauge bosons
The couplings of the heavy axion to the electroweak gauge bosons reflect the dependence of the mixing parameters on the extra source of mass M , as follows: The M -dependent corrections in these couplings can be relevant for heavy axions which solve the strong CP problem as well as for gluonic ALPs, as long as their mass parametrized by Eq. (3.3) is sensibly larger than that for the QCD (i.e. invisible) axion, M 2 > m 2 π f 2 π /f 2 a . These expressions hold as long as chiral perturbation theory is valid, that is M 1 GeV.
For values of M noticeably larger than the η mass the model prediction depicted is only indicative, as the effective Lagrangian in terms of pions and η is not really adequate and the description should the be done in terms of the couplings to quarks. At those energies (m a m η ) QCD is perturbative and it would be pertinent to compute the two-loop contribution of the gluonic coupling to the EW gauge boson couplings. For the case of photons, a qualitative estimation has been performed in Ref. [21] with the result: leading to a photonic axion coupling given by The derivation of the equivalent formula for the coupling of axions to heavy EW gauge boson couplings is left for future work. Nevertheless, the analysis presented here conveys the qualitative behaviour expected for the transition between the low and high axion mass regimes.

Gauge invariance
For high axion masses (i.e. M Λ QCD in Eq. (3.3)), the mixing of the axion with the SM pseudoscalars becomes negligible. For those scales, QCD enters the perturbative regime and Eq. (3.9) illustrates how the model-independent contributions to the EW couplings diminish. As the latter become negligible, the axion coupling to EW gauge bosons is parametrized by just the two effective interactions in Eqs. (2.5) and (2.6). In other words, the heavy axion EW couplings span a parameter space with two degrees of freedom (instead of three for light axions with m a Λ QCD , see Sect. 2.2). Two independent constraints follow for heavy axions: where we have chosen to express the couplings g aW W , g aZZ in terms of the overall better constrained ones g aγγ and g aγZ . 7 These powerful relations will be exploited in the next section to further constrain uncharted regions of the experimental parameter space.
Alike to the discussion after Eq. (2.63), the relations in Eqs. (3.11) and (3.12) apply not only to heavy axions and heavy gluonic ALPs, but also to generic ALPs which only exhibit EW interactions and are much heavier than nucleons. The corollary that at least two EW couplings -if any-must exist for any axion or ALP holds as well.

Phenomenological analysis
The impact of the results obtained above on present and future axion searches will be illustrated in this section. Both tree-level and loop-level effects will be taken into account. Indeed the latter are relevant when confronting data on photons, electrons and nucleons, as the experimental constraints on these channels are so strong that they often dominate the bounds on EW axion couplings.

Loop-induced couplings
The tree-level coupling of the axion to photons can be suppressed in some situations [9,17] (photophobic ALPs are also possible [24]). Additionally, many models have no tree-level couplings to leptons or suppressed couplings to nucleons. [8] However, all possible effective couplings will mix at the loop level. This affects the renormalization group (RG) evolution, via which all couplings allowed by symmetry will be generated even when assuming only a subset of couplings at some scale.
Before proceeding with the phenomenological analysis, we discuss in this subsection the loopinduced effective interactions arising from the direct coupling to electroweak gauge bosons. Because the experimental and observational limits are usually strongest for photons, electrons, and nucleons, the loop-induced contributions to these channels can give stronger constraints than those stemming from the tree-level impact on other channels.
Denoting the effective axion-fermion Lagrangian by it has been shown [21] that the coefficient c f 1 eff receives one loop-induced corrections from elec- 7 Obviously, the milder constrain in Eq. (2.63) also applies here.
troweak couplings, To obtain this result, the ultraviolet scale inside the loop diagrams has been identified with the axion scale f PQ , which is the cutoff of the effective theory. Note that the one-loop induced contributions to the fermion couplings are independent of the axion mass (other than a negligible dependence through the axion-gauge couplings such as g aγγ , see below). They are generically of order α/4π as expected, that is, over two orders of magnitude smaller than the original effective gauge coupling. Nevertheless, the experimental constraints are so strong that they will often provide the leading bounds on gauge-axion couplings. The most relevant fermionic limits are those on the coupling to electrons and light quarks. While Eq. (4.2) is directly applicable to leptons and heavy quarks, at low energies light quarks form hadrons: the loop-induced coupling to nucleons have highest impact. Following Refs. [8,13], the relation between nucleon and light quark couplings can be written as In Eq. The axion-photon coupling also receives one-loop corrections in the presence of couplings of the axion to fermions or gauge bosons. For energies or masses in the loops higher than Λ QCD , quarks are the appropiate propagating degrees of freedom, 8 and the one-loop contributions for an on-shell axion can be expressed as [21] g eff where the subscript f runs over leptons and heavy quarks. Here, τ i = 4m 2 i /m 2 a and Asymptotically, B 1 → 1 in the limit m a m f and B 1 → −m 2 a /(12m 2 f ) when m a m f . This means that the contribution of fermions heavier than the axion is strongly suppressed. Similarly, B 2 → 1 + π 2 /4 − log 2 m a /m W when m a m W , whereas B 2 → m 2 a /(6m 2 W ) in the m a m W limit.
It is worth noting that chaining the two previous one-loop contributions gives an approximate estimation of the two-loop contribution of a given heavy gauge boson coupling to the axion-photon coupling. As an example, consider g aW W in Eq. (4.2): it results in an effective fermion coupling c f 1 eff which, when subsequently inserted in the second term in Eq. (4.5), results in an effective axion-photon coupling. This can be compared with the third term which gives directly a one-loop contribution of g aW W to g aγγ : for m a m W , the second term in Eq. (4.5) may actually be numerically larger than the third term, that is, the two-loop contribution via fermionic couplings may dominate over the one-loop gauge one, as it was pointed out in Refs. [21,24]. Indeed, this two-loop contribution may be phenomenologically the most relevant one to constrain the axion couplings to heavy electroweak gauge bosons. A caveat is that only a true two-loop computation may settle the dominant pattern, but the analysis discussed is expected to provide an order of magnitude estimate.

Axion decay channels and lifetime
The plethora of couplings discussed, contributing either at tree or loop level, produces a rich variety of production and decay channels of the axion, depending on its mass and on the relative strength of the couplings. A quantitative evaluation of the lifetime and branching ratios is essential for assessing what experiments or searches are more adequate to test different regions of parameter space.
In order to determine the detection capabilities of a given final state channel, an important element is whether the axion can decay into it within the detector, or escapes and contributes to an "invisible" channel. For purely illustrative purposes, Fig. 2 compares the axion decay widths into different final states for a particular choice of the model-dependent parameters. 9 This figure serves to indicate the mass threshold for the different channels and is also a good indicator of the relative width of each decay channel. The determination of the areas experimentally excluded -to be developed in Sec. 4.3-will not depend on the value of the effective couplings assumed in this figure, though. 10 In the low mass region m a < 3m π , only decays to pairs of electrons, muons or photons are possible. The axion typically becomes long lived enough so as to be stable at collider and flavour experiments. Note that this region is particularly sensitive to a possible cancellation/suppression of the photonic coupling g aγγ (this happens for instance in models of axions in which the modeldependent parameter E/N partially cancels the model-independent contribution, see Eqs. (2.55) and (3.6)). This would suppress the decay width to photons and thus enhance the branching fraction to fermions, especially close to the respective mass thresholds.
The hadronic channel plays a central role as soon as it opens. It then dominates the decay of the axion due to the large gluonic coupling. The lightest possible hadronic final state is three neutral pions. At around the GeV scale many other final states become viable, but in this region chiral perturbation theory starts to break down and we refrain from making any quantitative predictions. 11 At high axion masses above 3 GeV the inclusive decay to hadrons can be safely estimated within perturbative QCD. 12 At much higher energies, tree-level decays to pairs of EW gauge bosons become possible and, though subdominant with respect to the hadronic one, will play a relevant role in collider searches.

Experimental constraints on the (heavy) axion parameter space
We have reinterpreted a number of axion searches into our framework. Far from an in-depth review, this study primarily intends to point out the relative strength of the different observables in probing different flavours and parameter regions of axion and ALP models. Interestingly and contrary to common lore, we find that some regions of parameter space can be better tested through the axion couplings to heavy gauge bosons rather than to photons.
The coloured areas in Fig. 3 show the regions experimentally excluded if taking into account exclusively the axion-gluon coupling g agg (which in axion models fixes the axion scale f a ). Although this work focuses on the case where also EW gauge boson couplings are present, this parameter space is also shown for reference.
The couplings of axions to EW gauge bosons (g aγγ g aW W , g aZZ and g aγZ ) will be instead explored within a two-operator framework: the axion-gluon coupling g agg and one electroweak gauge coupling are to be simultaneously considered. In other words, for each EW axion coupling 9 The value for gaγγ used in this figure corresponds approximately to fa = 10 TeV and E = 1 in axion models. Rescaling for other values of the couplings can be achieved by taking into account that Γi ∝ g 2 aXY . 10 The widths used in determining the coloured regions in Figs. 4 induced at tree-level derive directly from the effective gaXY value for each point of the parameter space. 11 After this paper was completed, Ref. [47] appeared which discusses the hadronic decays of axions in this region 1 GeV < ma < 3 GeV. 12 Note that heavy-flavour tagging can allow to distinguish final states involving heavy quarks, but this separation will not be taken into account here.  Figure 3: Coupling to gluons. The excluded areas rely on considering exclusively the coupling g agg by itself. The SN1987a limit was computed in Ref. [15], while the rare decays and LHC exclusions are obtained using the results of Ref. [48] and Refs. [20,49], respectively. Particular models are represented by the overlaid lines, the black one corresponding to the QCD axion and the blue one to heavy axion models (here for a benchmark choice of f a = 10 TeV). See Sec. 4.4 for more details. We provide exclusion plots without superimposed lines as auxiliary files.
g aXY , the regions experimentally excluded will be determined assuming the Lagrangian 13 This choice is mainly motivated by the focus on solving the strong CP problem, or alternatively as an ALP analysis which goes beyond the traditional consideration of only one operator at a time.
The regions experimentally excluded for axion-EW gauge boson couplings are then depicted in Fig. 4 as coloured areas. Note that we don't discuss any cosmological bounds. The reason for this is that models of heavy axions, typically containing an extended confining sector, are expected to significantly alter the standard cosmological picture that is usually assumed to obtain such exclusions. The study of heavy axion cosmology thus requires a full self-consistent study which is left for future work. The resulting greenish regions in Fig. 4a match well-known exclusion regions for g aγγ , although the overlap is not complete because the latter typically only take into account the effective axionphoton coupling; the additional presence in our analysis of the axion-gluon coupling g agg has a particularly relevant impact in the heavy axion region (see below).
Figs. 4b, 4c and 4d respectively for g W W , g ZZ and g γZ are novel. The possibility of measuring four distinct EW observables offers a multiple window approach and a superb cross-check if a signal is detected, given the fact that for axion masses much smaller (larger) than Λ QCD only three (two) couplings are independent among the set {g aγγ g aW W , g aZZ , g aγZ }, see Eq. (2.63) (Eqs. (3.11) and (3.12)).
For the majority of the regions excluded in Fig. 4, the experiment directly constraints certain regions of the parameter space {g aXY , m a } and no further assumptions are required; those constraints apply then also to ALPs which have no gluonic couplings. However, for some collider searches the interplay between the particular EW coupling g aXY under analysis and the coupling to gluons plays a relevant role, and it is necessary to assume their relative strength. This will be taken as given by This is well motivated by pseudo Nambu-Goldstone bosons with anomalous couplings generated by the triangle diagram with O(1) group theory factors. In any case, the results are largely insensitive to this assumption, since in the best limits stemming from LHC searches the production cross section times branching ratio scales as It is then enough to adopt the reasonable assumption that the coupling to gluons is larger than that to EW gauge bosons g agg g aXY . This assumption has been adopted exclusively for LHC searches in the axion mass region where the hadronic decay channels are open, i.e. m a > 3 m π .
In addition, the bounds obtained via loop contributions have a logarithmic dependence on the cut-off scale f a . A relation between f a and g aXY needs to be assumed in order to translate the bounds on fermionic or photonic couplings into bounds on EW gauge boson couplings. In these cases, Eq. (4.8) will be used again, which for axion models translates into where Eq. (2.54) has been used. In any case, as the cutoff dependence in loops is logarithmic this assumption has a minor impact on the exclusion plots.
To sum up, for each EW-axion coupling the experimentally excluded regions in Fig. 4 are depicted on the parameter space {g aXY , m a } without any assumption, except: A For LHC searches and m a > 3 m π , Eqs. (4.8) and (4.9) were used, which for most cases is equivalent to assume g agg g aXY .
B For the regions labelled as "photons", "electrons" and "nucleons" in Figs. 4b , 4c and 4d, the loop-induced bounds have a very mild dependence on the assumption in Eq. (4.10).
After having presented the general strategy that we use for the reinterpretation of constraints into our setup, we proceed to describe the origin of each of the exclusion regions coloured in Fig. 4.

Coupling to photons
The combination of astrophysical and terrestrial probes makes this search a particularly powerful tool to test the axion paradigm, especially for low mass axions. Even in the case of relatively large axion masses and/or situations where the coupling to photons can be suppressed, photons still place strong constraints both at tree-level and through loop-induced effects.
The experimental limits on g aγγ are compiled in Fig. 4a. At the lowest axion masses considered here, m a 10 meV, the most competitive bounds come from the CAST helioscope [50], and will improve in the future with the upgrade to the IAXO experiment. [51] At slightly larger masses up to m a ∼ 1 keV, g aγγ is constrained by an energy-loss argument applied to Horizontal Branch (HB) stars. [14] A similar argument applies to the supernova SN1987a and constrains larger masses up to the 100 MeV range, both using extra cooling arguments [52] and by the lack of observation of a photon burst coming from decaying emitted axions. [53] In the same mass range, larger couplings can be constrained using beam dump experiments, with these exclusions led at present by the     In panels b, c and d, the loop-induced constraints labelled as "photons", "electrons" and "nucleons" depend mildly on f a (see text). The label "g.i." illustrates the exclusion power of EW gauge invariance. The superimposed lines correspond to KSVZ-type QCD axion models (black line and grey band) and to one benchmark heavy axion model with f a = 10 TeV (blue lines). The parameter space for heavy axions corresponds to moving up and down that set of blue lines, see Sec. 4.4. The results apply as well to gluonic ALPs (we provide exclusion plots without superimposed lines as auxiliary files).
NuCal [54] experiment together with the 137 [55] and 141 [56] experiments at SLAC. We adapt here the constraints compiled in Ref. [57], noting that these bounds rely solely on the photon coupling.
For yet higher axion masses, where colliders provide the best limits, the gluon coupling plays a relevant role. As long as no hadronic decay channel is open, the LEP constraints based on Z → γγ and Z → γγγ searches obtained in [19,20] for ALPs without gluonic couplings are also applicable to heavy axions. However, for masses above 3m π , hadronic final states start to dominate and we refrain from claiming any exclusion, as a new dedicated analysis would be required to use these channels. This explains the white gap just left of the grey hatched "no human's land" region, which should be at least partially covered when the forementioned analysis is performed. It is nevertheless possible to exploit some collider searches, if a relation between the gluonic and the EW couplings is assumed. Assumption A is adopted here. Our reinterpretation of the analysis in Ref. [58], in which the L3 collaboration looked for hadronic final states accompanied by a hard photon, yields the limit labelled "L3" in Fig. 4a, though it is ultimately superseded by LHC exclusions. The region labelled "Flavour" is excluded by data from Babar [59] and LHCb [60], as computed in Ref. [40]. For high axion masses near the TeV scale, the limits from LHC are much stronger than those from LEP, because of the enhanced axion production via gluon-gluon fusion. We have included the limits obtained in this context in Refs. [61,62] using run 1 data. These limits are extremely strong and should improve with the addition of run 2 data, especially at higher energies.
Finally, the bounds on g aγγ described above have been translated -using assumption B-into competitive bounds for the other EW axion couplings, by means of their loop-mediated impact.

Coupling to fermions
Flavour blind observables involving fermions can be used to constrain gauge boson couplings via the impact of the latter at loop level, see Eq. (4.2). In order to fix the mild logarithmic dependence on the cutoff scale, assumption B will be adopted.
The most relevant constraints on flavour-blind axion-fermion interactions are of astrophysical origin and come from either electrons or nucleons. Firstly, a coupling of the axion to electrons allows for efficient extra cooling of some stars, which allows to place a bound on the axion-electron coupling g aee via the observation of Red Giants (RG) [14]. Secondly, and in a manner similar to the discussion above for photons, a too strong coupling of the axion to nucleons would have shortened the duration of the neutrino burst of the supernova SN1987a. We use the most recent evaluation of this bound calculated in Ref. [15]. These two observations (RG and SN1987a) give the strongest limits on the coupling of axions to gauge bosons for axion masses respectively below 10 keV and 10 MeV, as can be seen in Fig. 4.
In addition, the one-loop induced fermion couplings also play a role in many of the observables considered here. In particular, they open potential axion decay channels into pairs of fermions.

Rare decays
For axion masses in the MeV-GeV range, g aW W is best tested by its one-loop impact on rare meson decay experiments, where axions can be produced in flavour-changing neutral current (FCNC) processes. This search was first proposed in Ref. [25] (where ALPs either stable or decaying to photons were considered). Recently, these bounds have been recomputed in Ref. [24] in the context of photophobic ALPs, considering as well the potential decays of the axion to a pair of fermions due to the one-loop induced coupling in Eq. (4.2). We reinterpret these searches, taking into account in addition the effects of the gluonic axion coupling under the assumptions A and B. The main consequence is that, for axion masses m a > 3m π , the sensitivity is drastically reduced because of the opening of hadronic axion decay channels.
At low axion masses below 2m µ , the axion is long-lived enough so that it can be considered stable for experimental purposes. This means that in these regions the axion has to be looked for in invisible searches. The most stringent limits were placed by the E787 and E949 experiments testing the K + → π + X channel, with X invisible. [48] Following Ref. [25], we reinterpret this search in terms of axions coupled to W bosons, which yields the constraint shown in Fig. 4b. These bounds will be improved in the near future by the NA62 experiment.
Axions can also be produced from rare meson decays in proton beam dump experiments, where they can be looked for in searches for long-lived particles. The current best limits are set by the CHARM experiment [63], where the axion can be produced in Kaon and B meson decays and subsequently decays to a pair of electrons or muons. The framework developed in Ref. [64] has been recast to obtain the limit on g aW W shown in Fig. 4b. 14 Direct couplings to heavy gauge bosons LEP provides the best environment to directly test the g aγZ coupling for axion masses below m Z , as shown in Fig. 4d. The first constraint set assuming only the g aγZ coupling was placed in Ref. [23] exploiting the limit on the uncertainty of the total Z width, [66] Γ(Z → BSM) 2 MeV at 95% C.L., which allows to set a conservative bound on the process Z → aγ. Stronger limits can be placed by more specific searches, as studied in Ref. [24]. The best limit at axion masses low enough for the axion to be long-lived stem from the Z → γ + inv. search. For higher axion masses, the large hadronic branching fraction makes the Z → γ + had. search the more fruitful one to look for axions. Under assumption A for the relative strength of the gluonic and EW couplings, we exploit the results of the search performed by the L3 collaboration as presented in Ref. [58] to obtain strong limits for m a in the range from 10 GeV up to the Z mass. Note that, even if the search is the same than that used to constrain the photonic axion coupling, the exclusion for g aγZ has a larger reach due to the fact that the process is mediated by an on-shell Z boson, instead of a very virtual photon.
LHC allows to look for a plethora of processes sensitive to axions. In particular, for heavy axions it provides the best limits on the coupling to heavy EW gauge bosons. The drawback of restricting the analysis to processes that separately involve only one of the EW gauge boson couplings plus the gluon coupling is the reduced number of available searches. Nevertheless, the advantage is that it provides robust and model-independent constraints.
The authors of Ref. [23] studied the LHC phenomenology of axions that are stable on collider lengths and thus would manifest themselves as missing energy. In particular, mono-W and mono-Z final states where an axion is radiated as missing energy/momentum can set constraints on the three couplings g aW W , g aZZ and g aγZ , as shown in Fig. 4. For large axion masses m a > m Z , the authors of Ref. [24] suggested that triboson final states place the strongest bounds on ALPs coupling to massive gauge bosons, though the sensitivity of this search is hindered for axions because of the large hadronic branching ratio that we take into account. Adapting their constraints -with assumption A-leads to the exclusion of regions in parameter space near the TeV range, as shown in Figs. 4b and 4d for g aW W and g aγZ , respectively. Note that significant exclusions can also be placed through the loop-induced coupling to photons. However, the most promising LHC search is one that -to the best of our knowledge-has not been performed yet. We advocate [67] the use of pp → a → V V processes, which benefit from the large production cross section through the gluonic coupling together with the clean final states that the decay to EW gauge bosons produce. We foresee that this search will have a sensitivity to the couplings of axions to heavy EW gauge bosons similar to the photonic case presented in Fig. 4a. Though potentially very interesting, the detailed analysis that this study requires is beyond the scope of this work and is left for the future. [67]

Impact on (heavy) axion models and gluonic ALPs
The black oblique line in Fig. 3 corresponds to the linear relation between 1/f a and m a for the QCD axion, Eq. (1.1). The horizontal blue branch is one example of how that relation changes after Eq. (3.3) for an illustrative example of a true heavy axion.
The black, grey and blue lines in Fig. 4 illustrate possible {m a , f a } values for axions which have a gluonic coupling g agg (and thus may solve the strong CP problem) in addition to at least one coupling to heavy gauge bosons. Those model-dependent lines are superimposed 15 on the colored/white regions excluded/allowed by experiments for each one of the couplings in the set {g aγγ , g aγZ , g aZZ , g aW W }, as determined above. The examples chosen corresponds to KSVZ-type axions: either a standard QCD axion or a heavy axion as in theories with an enlarged confining gauge sector.   Fig. 4.
In each panel, for a given value of the model-dependent coupling: -The expectation for the pure QCD axion is depicted by grey and black lines. The bands in Fig. 4 delimited by grey lines correspond to just one exotic KSVZ fermion representation.
The values of the model-dependent parameters delimiting these benchmark bands [17] are summarized in Table 1. The black line is instead an illustrative case with two fermion representations such that the coupling to photons g aγγ cancels up to theoretical uncertainties. [17] The upward bending of the lines obeys the expected change of the prediction for axion masses larger than the η mass, a regime in which the last term in the parentheses in Eqs. , and it is often higher than for heavy axions lighter than the pion, for which the partial cancellation between the model-dependent and model-independent mixing effects may operate.
The figures illustrate that the M -dependent corrections may be relevant even for not very large M values. For instance, two close values of the model-dependent photon couplings E/N may give a very close g aγγ prediction for M values above the η mass, while that prediction can widely differ for M values smaller than the QCD scale. This is clearly reflected by the lines corresponding to the smaller values of E/N in Fig. 4a. The parameter space for heavy axion models spans in fact most of the region to the right of the oblique band for the QCD axion: parallel sets of horizontal lines above and below the blue ones depicted are possible and expected for other values of the heavy axion parameters. For a given f a , varying M (that is, varying m a ) is tantamount to move right or left on a horizontal blue line, while varying f a displaces up or down the set of horizontal blue lines. Finally, all these considerations for heavy axion models apply as well to gluonic ALPs, as argued earlier.

Gauge invariance
For heavy axions or any type of ALP with masses m a Λ QCD , couplings to EW gauge bosons are directly tested and gauge invariance imposes the two relations in Eqs. (3.11) and (3.12). Therefore, the combination of the experimental constrains on two of the operators in the set {g aγγ , g aW W , g aZZ , g aγZ } translates in model independent bounds on the other two couplings. For light masses m a ≤ Λ QCD , only Eq. (2.63) applies instead. These bounds based solely on EW gauge invariance have been depicted by black curves on the upper right corner of Figs. 4b, 4c and 4d. They susbtantially reduce the latter parameter space, especially in the cases of g aW W and g aZZ , whose current direct constraints are less powerful. They are to be taken with caution, though, since in each of the exclusion plots only one EW coupling was taken into account, while the relations deduced from gauge invariance involve several non-vanishing EW couplings. Furthermore, in a future multiparameter analysis where tree-level axion-fermion couplings are included, those relations could be corrected via one-loop effects.

Implications for heavy axion models
In existing models that solve the strong CP problem with heavy axions, often either where Λ Λ QCD is a new strong confining scale, and the primed fields denote exotic "pions" corresponding to the exotic fermions charged under the extra confining force. Let us assume the first option as an example. In this case, 16 Assume that an experiment detects an axion-photon signal and a certain value for the axion mass which correspond to a point in the white region of Fig. 4a and located to the right of the QCD 16 The mum d /(mu + m d ) 2 factor in the QCD contribution is not shown for simplicity.
band. For instance, let us consider a point located on one of the flat sections of the blue lines depicted. The interpretation in terms of a heavy axion depends on whether the axion is heavier or lighter than the pion and η , respectively: m a m η . The model-independent effects due to the mixings with SM pseudoscalars have become negligible, and g aγγ is a direct measure of the product (1/f a )E/N .
m a m π . In this case the measured g aγγ value is undistinguishable from that for the QCD axion, with the E/N and f a dependence given by Eq. (2.55). In other words, it would indicate either a heavy axion or a QCD axion with some degree of photophobia, as in that region they become undistinguishable. This is so at least for the lowest axion masses and/or without the help of other measurements involving heavy EW gauge bosons.
In both cases, in the framework of a KSVZ model, an additional measurement of the axion coupling to heavy gauge bosons would be enough to disentangle the values of f a and M , that is, to determine the high scale Λ . The model-independent corrections determined in Eqs. (3.5)-(3.8) can be essential when exploiting low-energy processes (e.g. rare decays), specially when they lead to the cancellation of a given coupling. Such cancellation in a channel in general will not apply to the couplings of the axion to other gauge bosons. Overall, the fact that only two or three axion-EW couplings are independent, while four channels can be explored, will allow to overconstrain the system.

Conclusions
Among the novel results of this work, we have first determined at leading order in the chiral expansion the model-independent components of the coupling of the QCD axion to heavy EW gauge bosons: g aγZ , g aZZ and g aW W . They stem from the axion-η -pion mixing induced by the anomalous QCD couplings of these three pseudoscalars. Our results extend to heavy EW gauge bosons the well known result for the photonic coupling of the axion g aγγ . They must be taken into account whenever an axion lighter than Λ QCD is on-shell and/or the energy and momenta involved in a physical process are of the order of the QCD confining scale or lower. As a previous step, we re-derived pedagogically the leading contributions to g aγγ for the case of the most general axion couplings (App. C), and then proceeded to the determination of the couplings to the SM heavy gauge bosons.
This analysis of the EW couplings of the QCD axion may have rich consequences when comparing the presence/absence of signals at two different energy regimes. For instance, the axion could be photophobic at low energies [17] or even EW-phobic (e.g. in rare meson decays) because of cancellations between the model-independent and model-dependent components, while an axion signal may appear at accelerators or other experiments at higher energies at which the model-independent component disappears.
We have next extended those results to the case of heavy axions which solve the strong CP problem. This has allowed to explore how the mixing of the axion with the pion and η evolves with rising axion mass, and in consequence how the model-independent contributions to all four EW axion couplings vanish as the axion mass increases above the QCD confinement scale. We have determined the modified expression for g aγγ relevant for heavy axions, which may have rich consequences: an hypothetical measurement of that coupling outside the QCD axion band could point to either a heavy axion or a photophobic QCD axion. The analogous expressions for g aγZ , g aZZ and g aW W have been also worked out.
On the purely phenomenological analysis, we developed a "two simultaneous coupling" approach in order to determine the regions experimentally excluded by present data for g aγγ , g aγZ , g aZZ and g aW W versus the axion mass. Each EW coupling has been considered simultaneously with the anomalous gluonic coupling essential to solve the strong CP problem. The allowed/excluded experimental areas have been depicted for each of those couplings as a function of the axion mass. This is the first such reinterpretation for g aγZ , g aZZ and g aW W . Even for g aγγ , the results of previous studies often did not apply and must be reanalysed: for instance the present bounds extracted from LEP and LHC data tend to focus on ALPs which would not have gluonic couplings, with very few exceptions [20,40,61,62]. Furthermore, we have included an estimation of the one-loop induced bounds for each EW axion coupling, which leads to supplementary constraints.
The expectations from KSVZ-type of theories have been then projected and illustrated over the obtained experimental regions, both for the QCD axion and for heavy axions. The compatibility of hypothetical a priori contradictory signals in high and low energy experiments in terms of a given axion has been pointed out. Furthermore, we discussed how to interpret an eventual signal (or null result) outside the QCD axion band in terms of the value of the new high confining scale generically present in heavy axion theories.
A simple point with far reaching consequences results from EW gauge invariance. In all generality, not all couplings of axions to EW gauge bosons are independent among the four physical ones in the set {g aγγ g aγZ , g aZZ and g aW W }. In particular, the relations obtained imply that at least two EW gauge couplings -if any-must be non-vanishing for any axion or ALP. These facts have been used to project the exclusion limits for the presently best constrained couplings onto the parameter space for the less constrained ones. In particular, for axions/ALPs much heavier than Λ QCD , those relations have been projected on the parameter space for g aW W , g aZZ and g aγZ , reinforcing their constraints. A future multiparameter analysis may correct them via loop corrections. Nevertheless, the results obtained here clear up the uncharted experimental regions and may be of use in setting a search strategy. More in general, the existence of four physical axion couplings to EW bosons at experimental reach constitutes a phenomenal tool to over-constrain the axion parameter space and to check the origin of an eventual axion signal.
Finally, all results obtained for heavy axions apply as well to ALPs which have both EW and gluonic anomalous couplings. The constraints stemming from EW gauge invariance extend even to generic ALPs which do not couple to gluons. In consequence, this work also extends automatically the usual parameter space for ALPs that do not intend to solve the strong CP problem, adding to the incipient efforts to go towards a multi-parameter strategy.
gional" (FEDER) through the project FPA2016-78645-P. The work of G. A. was supported through an ESR contract of the H2020 ITN Elusives (H2020-MSCA-ITN-2015//674896-ELUSIVES) and through a "La Caixa" predoctoral grant of Fundación La Caixa. The work of P.Q. was supported through a "La Caixa-Severo Ochoa" predoctoral grant of Fundación La Caixa.

A The neutral pseudoscalar mass matrix
In the expansion in which all terms appearing at first order in 1/f a and in quark masses are kept, Eq.

B Anomalous couplings of the pseudoscalar mesons to the EW gauge bosons
In addition to the axion couplings in Eqs. (2.4)-(2.6), new couplings involving gauge bosons are expected as soon as the extra pseudoscalar mesons appear in the spectrum below the confinement scale. In particular, for the η 0 and the pions, the following interactions with EW gauge bosons (considered as external currents) are possible below the QCD confinement scale: W a µνW aµν η 0 f π , B µνB µν η 0 f π , B µνB µν π 3 f π , W a µνB µν π a f π . (B.1) The lightness of pseudoscalar mesons in QCD appears as a natural consequence of the spontaneous breaking of the chiral symmetry U (2) L × U (2) R −→ U (2) V . The pions and eta mesons can be identified as the pseudo Nambu-Goldstone bosons of the broken symmetry. However, the chiral symmetry is only approximate. It is explicitly broken not only by the quark masses, but also by the electroweak gauge interactions. This fact famously explains the difference between charged and neutral pion masses but also allows the computation of the coupling of pseudoscalar mesons to the EW bosons through the anomaly, where q u and q d are the electromagnetic charges of the up and down quarks, we find, Recalling that the coefficients c u, d 2 correspond to the PQ charges, the combinations that appear in the above equation can be identified as the electromagnetic and QCD anomaly coefficients for the up and down quarks, Redefining the axion decay constant as usual f a = f PQ /N , we can express To sum up, from the most general mass matrix we have obtained the same result of Eq. (2.55) taking into account that in the E/N we have to sum over all fermions transforming under the PQ symmetry, that means including the up and down quarks E = E heavy + E u, d ,