Collider Probes of Axion-Like Particles

Axion-like particles (ALPs), which are gauge-singlets under the Standard Model (SM), appear in many well-motivated extensions of the SM. Describing the interactions of ALPs with SM fields by means of an effective Lagrangian, we discuss ALP decays into SM particles at one-loop order, including for the first time a calculation of the $a\to\pi\pi\pi$ decay rates for ALP masses below a few GeV. We argue that, if the ALP couples to at least some SM particles with couplings of order $(0.01-1) \mbox{TeV}^{-1}$, its mass must be above 1 MeV. Taking into account the possibility of a macroscopic ALP decay length, we show that large regions of so far unconstrained parameter space can be explored by searches for the exotic, on-shell Higgs and $Z$ decays $h\to Za$, $h\to aa$ and $Z\to\gamma a$ in Run-2 of the LHC with an integrated luminosity of 300 fb$^{-1}$. This includes the parameter space in which ALPs can explain the anomalous magnetic moment of the muon. Considering subsequent ALP decays into photons and charged leptons, we show that the LHC provides unprecedented sensitivity to the ALP-photon and ALP-lepton couplings in the mass region above a few MeV, even if the relevant ALP couplings are loop suppressed and the $a\to\gamma\gamma$ and $a\to\ell^+\ell^-$ branching ratios are significantly less than 1. We also discuss constraints on the ALP parameter space from electroweak precision tests.

New pseudoscalar particles with masses below the electroweak scale appear frequently in well-motivated extensions of the Standard Model (SM). Examples are axions [1][2][3][4][5][6][7][8] addressing the strong CP problem or pseudoscalar mediators of a new interaction between dark or hidden sectors and the SM [9]. Further, various anomalies can be explained by the presence of new spin-zero states with pseudoscalar couplings. Examples are the longstanding deviation of the anomalous magnetic moment of the muon from its SM value [10,11], or the excess in excited Beryllium decays 8 Be * → 8 Be+e + e − recently observed by the Atomki collaboration [12][13][14]. Dark-matter portals with a pseudoscalar mediator lighter than the Higgs boson can address the gamma-ray excess observed in the center of the galaxy by the Fermi-LAT collaboration, while avoiding constraints from direct detection and collider searches [15,16]. Axion-like particles (ALPs) have triggered interest way beyond their potential relevance in the context of the strong CP problem [17,18]. Pseudo Nambu-Goldstone bosons arise generically in models with spontaneous breaking of a global symmetry. Due to an (approximate) shift symmetry they can naturally be light with respect to the electroweak or even the QCD scale. Low-energy observables, cosmological constraints and ALP searches with helioscopes probe a significant region of the parameter space in terms of the mass of the ALP and its couplings to photons and electrons. Future helioscope experiments like the International Axion Observatory (IAXO) [19], and beam-dump experiments such as the facility to Search for Hidden Particles (SHIP) [20,21], will further improve these constraints for ALP masses below a GeV. Collider experiments have searched directly and indirectly for ALPs [22]. Besides ALP production in association with photons, jets and electroweak gauge bosons [23][24][25][26], searches for decays of the Z boson into a pseudoscalar a and a photon at LEP and the LHC provide limits for ALPs with up to electroweak scale masses [24,[27][28][29]. Constraints from flavor-violating couplings have recently been summarized in [30]. Utilizing Higgs decays to search for light pseudoscalars has been proposed in [31][32][33][34][35]. Several experimental searches looking for the decay h → aa have been performed, constraining various final states [36][37][38][39][40][41][42]. Surprisingly, no experimental analyses of the decay h → Za exist, even though analogous searches for heavy resonances decaying into a Z boson and a pseudoscalar a [43] as well as a search for a light Z boson in h → ZZ decays [44] have been performed. The reason is, perhaps, the suppression of the h → Za decay in the decoupling limit in two-Higgs-doublet models in general and supersymmetric models in particular [45]. In models featuring a gauge-singlet ALP, there is no dimension-5 operator mediating h → Za decay at tree level, and hence this mode has not received much theoretical attention either (see, however, a recent brief discussion in [26]).
In this paper we present a comprehensive analysis of the on-shell Higgs decay modes h → Za and h → aa as well as the on-shell Z-boson decay Z → γa starting from a general effective Lagrangian for a gauge-singlet ALP interacting with SM fields. We show that these decays can be used to probe the ALP couplings to SM particles in regions of parameters space inaccessible to any other searches. A first exposition of the main ideas of our approach has been presented in [46]. In the present paper we extend this JHEP12(2017)044 the possible two-particle decays of ALPs is presented in section 3, where we consistently include the tree-level contributions and one-loop corrections to the decay amplitudes. For ALP masses below a few GeV, we calculate the a → πππ decay rates and the effective ALP-photon couplings using a chiral Lagrangian. We also survey present constraints on the ALP-photon and ALP-electron couplings and point out that, under the assumption that the ALP couples at least to some SM particles with couplings larger than about (100 TeV) −1 , its mass must be above 1 MeV. In section 4 the preferred region of parameter space in which an ALP can explain the anomalous magnetic moment of the muon is derived. Section 5 is devoted to a detailed discussion of the exotic Higgs decays h → Za and h → aa. We discuss which regions of parameter space can be probed with 300 fb −1 of integrated luminosity in Run-2 of the LHC, and which regions can already be excluded using existing searches. In section 6 we extend this discussion to the exotic decay Z → γa, and we study Z-pole constraints from electroweak precision tests. We conclude in section 7. Technical details of our calculations are relegated to four appendices.

Effective Lagrangian for ALPs
We assume the existence of a new spin-0 resonance a, which is a gauge-singlet under the SM gauge group. Its mass m a is assumed to be smaller than the electroweak scale. A natural way to get such a light particle is by imposing a shift symmetry, a → a + c, where c is a constant. We will furthermore assume that the UV theory is CP invariant, and that CP is broken only by the SM Yukawa interactions. The particle a is supposed to be odd under CP. Then the most general effective Lagrangian including operators of dimension up to 5 (written in the unbroken phase of the electroweak symmetry) reads [52] L D≤5 where we have allowed for an explicit shift-symmetry breaking mass term m a,0 (see below). G A µν , W A µν and B µν are the field strength tensors of SU(3) c , SU(2) L and U(1) Y , and g s , g and g denote the corresponding coupling constants. The dual field strength tensors are defined asB µν = 1 2 µναβ B αβ etc. (with 0123 = 1). The advantage of factoring out the gauge couplings in the terms in the second line is that in this way the corresponding Wilson coefficients are scale invariant at one-loop order (see e.g. [53] for a recent discussion of the evolution equations beyond leading order). The sum in the first line extends over the chiral fermion multiplets F of the SM. The quantities C F are hermitian matrices in generation space. For the couplings of a to the U(1) Y and SU(2) L gauge fields, the additional terms arising from a constant shift a → a + c of the ALP field can be removed by field redefinitions. The coupling to QCD gauge fields is not invariant under a continuous shift transformation because of instanton effects, which however preserve a discrete version of the shift symmetry. Above we have indicated the suppression of the dimension-5 operators with a new-physics scale Λ, which is the characteristic scale of global symmetry breaking, JHEP12(2017)044 assumed to be above the weak scale. In the literature on axion phenomenology one often eliminates Λ in favor of the "axion decay constant" f a , defined such that Λ/|C GG | = 32π 2 f a . Note that at dimension-5 order there are no ALP couplings to the Higgs doublet φ. The only candidate for such an interaction is where c w ≡ cos θ w denotes the cosine of the weak mixing angle, and the last expression holds in unitary gauge. Despite appearance, this operator does not give rise to a tree-level h → Za matrix element; the resulting tree-level graphs precisely cancel each other [47,48]. Indeed, a term C Zh O Zh in the Lagrangian is redundant, because it can be reduced to the fermionic operators in (2.1) using the equations of motion for the Higgs doublet and the SM fermions [47,48]. The field redefinitions with ξ = C Zh /Λ, eliminate O Zh and shift the flavor matrices C F of the SU(2) L singlet fermions by 1 while the matrices C Q and C L of the SU(2) L doublets remain unchanged. There are no additional contributions to the operators in (2.1) involving the gauge fields, because the combination of axial-vector currents induced by the shifts in (2.4) is anomaly free. In this work we will be agnostic about the values of the Wilson coefficients. We will show that ALP searches at high-energy colliders are sensitive to couplings C i /Λ ranging from (1 TeV) −1 to (100 TeV) −1 . In weakly-coupled UV completions one expects that the operators describing ALP couplings to SM bosons have loop-suppressed couplings (see e.g. [54] for a recent discussion). This is in line with estimates based on naive dimensional analysis, which we briefly discuss in appendix A. Departures from these estimates can arise in models involving e.g. large multiplicities of new particles in loops. It is common practice in the ALP literature to absorb potential loop factors that may arise into the Wilson coefficients C i . As we will discuss in section 4, the puzzle of the anomalous magnetic moment of the muon can be resolved within our framework if C γγ /Λ = O(1/TeV). Probing this region at colliders is thus a particularly well motivated target [46]. We emphasize, though, that by using the search strategies developed here it will be possible to probe even loop-suppressed couplings as long as the new-physics scale Λ is in the TeV range.
The ALP can receive a mass by means of either an explicit soft breaking of the shift symmetry or through non-perturbative dynamics, like in the case of the QCD axion [3,4]. In the absence of an explicit breaking, QCD dynamics generates a mass term given by [6,55,56] m a, dyn ≈ 5.7 µeV 10 12 GeV f a ≈ 1.8 MeV |C GG | 1 TeV Λ . (2.5)

JHEP12(2017)044
When an explicit symmetry-breaking mass term m a,0 is included in the effective Lagrangian (2.1), the resulting mass squared m 2 a = m 2 a,0 + m 2 a, dyn becomes a free parameter. We will assume that m a v. At dimension-6 order and higher, several additional operators can arise. The ALP couplings to the Higgs field are those most relevant to our analysis. They are (2.6) The first two terms are the leading Higgs portal interactions, which give rise to the decay h → aa. Note that the second term, which explicitly violates the shift symmetry, is allowed only if the effective Lagrangian contains an explicit mass term for the ALP. Its effect is suppressed, relative to the first term, by a factor m 2 a,0 /m 2 h . The third term is the leading operator mediating the decay h → Za at tree level [47,48]. These decay modes will be of particular interest to our discussion in section 5.
After electroweak symmetry breaking (EWSB), the effective Lagrangian (2.1) contains couplings of the pseudoscalar a to γγ, γZ and ZZ. The relevant terms read where s w = sin θ w and c w = cos θ w , and we have defined The fermion mass terms resulting after EWSB are brought in diagonal form by means of field redefinitions, such that U † u Y u W u = diag(y u , y c , y t ) etc. Under these field redefinitions the matrices C F transform into new matrices (2.9) In any realistic model these couplings must have a hierarchical structure in order to be consistent with the strong constraints from flavor physics. We will discuss the structure of the flavor-changing ALP couplings in a companion paper [57]. For now we focus on the flavor-diagonal couplings. Using the fact that the flavor-diagonal vector currents are conserved, we can rewrite the relevant terms in the Lagrangian in the form where the sum runs over all fermion mass eigenstates, and we have defined (with i = 1, 2, 3)  Figure 1. Representative one-loop Feynman diagrams contributing to the decay a → γγ. The internal boson lines represent charged W bosons and the associated charged Goldstone fields. The last diagram contains the (gauge-dependent) self-energy Π γZ (0). One also needs to include the on-shell wave-function renormalization factors for the external photon fields.
shift-invariant coupling of an ALP to neutrino fields arises at dimension-8 order from an operator consisting of a times the Weinberg operator. Even in the most optimistic case, where no small coupling constant is associated with this operator, the resulting a → νν decay rate would be suppressed, relative to the a → γγ rate, by a factor of order m 2 a v 4 /Λ 6 . Alternatively, if Dirac neutrino mass terms are added to the SM, the corresponding couplings in (2.10) yield a a → νν decay rate proportional to m 2 ν . In either way, for Λ in the TeV range or higher, this decay rate is so strongly suppressed that if an ALP can only decay into neutrinos (e.g. since it is lighter than 2m e and its coupling to photons is exactly zero for some reason) it would be a long-lived particle for all practical purposes.

ALP decay rates into SM particles
The effective Lagrangian (2.1) governs the leading interactions (in powers of v/Λ) giving rise to ALP decays into pairs of SM gauge bosons and fermions, while the additional interactions in (2.6) are needed to parametrize the exotic decays of Higgs bosons into final states involving an ALP. In computing the various decay rates, we include the tree-level and one-loop contributions from the relevant operators. We find that fermion-loop corrections can be numerically important, and they can even be dominant in new-physics models where the coefficients C V V in (2.1) (with V = G, W, B) are loop suppressed.

ALP decay into photons
In many scenarios, the di-photon decay is the dominant decay mode of a light ALP. Because of its special importance, we have calculated the corresponding decay rate from the effective Lagrangian (2.1) including the complete set of one-loop corrections. The relevant Feynman diagrams are shown in figure 1. We define an effective coefficient C eff γγ such that To an excellent approximation (apart from a mild mass dependence in the loop corrections) the a → γγ decay rate scales with the third power of the ALP mass. For a very light ALP with m a < 2m e this is the only SM decay mode allowed, and with decreasing ALP mass the decay rate will eventually become so small that the ALP will leave the detector and appear as an invisible particle.

JHEP12(2017)044
The expression for C eff γγ depends on the ALP mass. If m a Λ QCD , then all loop corrections, including those involving colored particles, can be evaluated in perturbation theory. We obtain where τ i ≡ 4m 2 i /m 2 a for any SM particle, and N f c and Q f denote the color multiplicity and electric charge (in units of e) of the fermion f . The loop functions read The fermion loop function has the property that for heavy fermions (m f m a ). Thus, each electrically charged fermion lighter than the ALP adds a contribution of order c f f /(16π 2 ) to the effective Wilson coefficient C eff γγ , while fermions heavier than the ALP decouple. The calculation of the electroweak loop corrections to the decay rate is far more involved than that of the fermion loops. We have evaluated the relevant diagrams shown in figure 1 in a general R ξ gauge. After some intricate cancellations, the main result of these corrections is to renormalize the fine-structure constant α in the expression for the rate, which is to be evaluated at q 2 = 0, as appropriate for on-shell photons. As mentioned earlier, the Wilson coefficient C γγ is not renormalized at one-loop order. The remaining finite correction in (3.2) is strongly suppressed, since the loop function B 2 (τ W ) ≈ m 2 a 6m 2 W is proportional to the ALP mass squared. An interesting feature of our result for the effective ALP-photon coupling in (3.2) is that the loop-induced contributions from both fermions and W bosons vanish in the limit m a → 0. This is an advantage of our choice of operator basis.
It is interesting to compare our result for the fermionic contributions to the a → γγ decay rate with the corresponding effects on the di-photon decay rate of a CP-odd Higgs boson. In this case the Higgs boson couples to the pseudoscalar fermion current, and one finds an expression analogous to (3.2), but with the loop function [ [58]. The difference can be understood using the anomaly equation for the divergence of the axial-vector current, which allows us to rewrite the ALP-fermion coupling in (2.10) in the form where the dots represent similar terms involving gluons and weak gauge fields. The first term on the right-hand side is now of the same form as the coupling of a CP-odd Higgs boson to fermions, while the second term has the effect of subtracting "1" from the function B 1 (τ f ).

JHEP12(2017)044
At one-loop order, relation (3.2) involves all Wilson coefficients in the effective Lagrangian (2.1) except for C GG . Even if the original coefficient C γγ vanished for some reason, these loop contributions would induce an effective coefficient C eff γγ at one-loop order. The ALP-gluon coupling would first enter at two-loop order. Using results derived in the following section, its effect can be estimated as where for the light quarks q = u, d, s one should use a typical hadronic scale such as m π instead of m q in the argument of the logarithm. Numerically, this two-loop contribution can be sizable due to the large logarithm.
If the ALP mass is not in the perturbative regime, i.e. for m a 1 GeV, the hadronic loop corrections to the effective ALP-photon coupling can be calculated using an effective chiral Lagrangian. This is discussed in detail in appendix B. Including interactions up to linear order in the ALP field, and working at leading order in the chiral expansion, one obtains [52] Here f π ≈ 130 MeV is the pion decay constant, Σ = exp i √ 2 fπ τ A π A contains the pion fields and B 0 = m 2 π mu+m d is proportional to the chiral condensate. For simplicity we restrict ourselves to flavor SU(2) with just one generation of light quarks. The hermitian matrices m q = diag (m u , m d ) andĉ qq = diag (c qq + 32π 2 C GG κ q ) are diagonal in the quark mass basis. The parameters have been chosen such that there is no tree-level mass mixing of the ALP with the π 0 [52]. Note the unusual appearance of a "tree-level" contribution proportional to C GG to the coefficient of the ALP-photon coupling in (3.6). When higher-order corrections (including the effects of the strange quark) are taken into account, the coefficient of C GG inside the bracket is reduced by about 5% and one obtains [C γγ − (1.92 ± 0.04) C GG ] [59]. This large effect is a consequence of the axial-vector anomaly leading to enhanced π 0 , η, η couplings to two photons combined with a mass-mixing of the ALP with these mesons [60]. QCD dynamics generates a mass for the ALP given (at lowest order) by [6,55,56] (3.8) A possible explicit shift-symmetry breaking mass term m 2 a,0 would have to be added to this expression. The last term in (3.6) gives rise to a kinetic mixing between the ALP and the neutral pion. The physical states are obtained by bringing the kinetic terms into canonical form and rediagonalizing the mass matrix. This changes the mass eigenvalues for π 0 and a by tiny corrections of order f 2 π /Λ 2 relative to the leading terms. At the same time, the state π 0 receives a small admixture of the physical ALP state, such that Relation (3.9) holds as long as |m 2 π − m 2 a | 2 m π m a . In the opposite limit one would obtain π 0 = 1 √ 2 (π 0 phys + a phys ) + O( ), but such a large mixing requires a fine-tuning of the masses that is rather implausible. In the presence of the mixing in (3.9), the SM π 0 → γγ amplitude mediated by the axial-vector anomaly induces an additional contribution to the a → γγ amplitude. Combining all terms, we obtain (assuming m a = m π ) (3.11) The contribution from the coefficient c ss not shown here would be suppressed, for light ALPs, by a factor of order m 2 π /m 2 η relative to the contributions from c uu and c dd .

ALP decays into charged leptons
If the ALP mass is larger than 2m e ≈ 1.022 MeV, the leptonic decay a → e + e − or decays into heavier leptons (if kinematically allowed) can be the dominant ALP decay modes in some regions of parameter space. We have calculated the corresponding decay rates from the effective Lagrangian including the complete set of one-loop mixing contributions from the bosonic operators in (2.1) and (2.7), see figure 2. In analogy with (3.1), we write the result in the form (with = e, µ, τ ) Here Q = −1 is the electric charge of the charged lepton, and T 3 = − 1 2 is the weak isospin of its left-handed component. In the limit where m 2 is either much smaller or much larger than m 2 a , the loop function in the photon term is given by (3.14) The exact expression is given in appendix C. In (3.13) we have regularized the UV divergences of the various contributions using dimensional regularization in the MS scheme.
Only the sum of all contributions is scale independent, i.e. the scale dependence of c (µ) compensates the scale dependence of the other terms. We do not show the one-loop corrections proportional to the tree-level coefficient c itself. They contain IR divergences, which cancel in the sum of the decay rates for a → + − and a → + − γ soft with a soft photon in the final state. The scheme-dependent constant δ 1 in (3.13) arises from the treatment of the Levi-Civita symbol in d dimensions, as we also discuss in appendix C. We obtain δ 1 = − 11 3 . In a scheme where instead the Levi-Civita symbol is treated as a 4-dimensional object, one would have δ 1 = 0.
Relation (3.13) shows two important facts: first, at one-loop order ALP couplings to fermions are induced from operators in the effective Lagrangian coupling the ALP to gauge bosons; and second, it would be inconsistent to set c to zero in (2.1), since this scaledependent coefficient mixes with the coefficients of bosonic operators under renormalization. Hence it must contain µ-dependent terms, which cancel the explicit scale dependence in the above result. Because of the presence of such terms, the only information that can conclusively be extracted from the calculation of the low-energy contributions performed above are the coefficients of the large logarithms obtained by identifying the factorization scale µ with the UV cutoff Λ. The result for these logarithmic contributions simplifies when one adds up the various terms in (3.13), since they can be derived in the unbroken phase of the electroweak theory. We obtain 15) where the first two terms arise from the loops of SU(2) L and U(1) Y gauge bosons, for which tr(τ A τ A ) = 3 4 and Y L = − 1 2 , Y R = −1. The last term contains the finite large logarithm related to the long-distance photon contribution, with C γγ given in (2.8).

ALP decays into hadrons
At the partonic level, the pseudoscalar a can also decay into colored particles. At tree-level the relevant modes are a → gg and a → qq. In the hadronic world these decays are allowed only if m a > m π . However, below 1 GeV the number of possible hadronic decay channels is very limited, because the two-body decays a → ππ and a → π 0 γ are forbidden by CP invariance and angular momentum conservation, while the three-body modes a → ππγ, a → π 0 γγ and a → π 0 e + e − are strongly suppressed by phase space and powers of the fine-structure constant α [31]. The dominant decay modes in this region are a → 3π 0 and a → π + π − π 0 . As long as the ALP is sufficiently light, so that the energy of the final-state mesons is much less than 4πf π ≈ 1.6 GeV, the calculation of the decay rates for exclusive modes such as a → πππ can be performed using the effective chiral Lagrangian (3.6). ALP couplings to three pions arise from each of the three terms shown in the second line of this equation, where in the first two terms one must substitute relation (3.9) for the π 0 fields. Working consistently at leading order in the chiral expansion, we obtain where (with 0 ≤ r ≤ 1/9) (3.17) Both functions are normalized such that g ab (0) = 1, and they vanish at the threshold r = 1/9. If the ALP mass is in the perturbative regime (i.e., for m a Λ QCD ), its inclusive decay rate into hadrons can be calculated under the assumption of quark-hadron duality [62,63]. Setting the light quark masses to zero (since here by assumptions m a m q for all light quarks) and including the one-loop QCD corrections to the decay rate as calculated in [58], we obtain where n q = 3 is the number of light quark flavors. To good approximation this rate scales with the third power of the ALP mass. Decays into heavy quarks, if kinematically allowed, can be reconstructed by heavy-flavor tagging. The corresponding rates are (with Q = b or c) where at leading order in perturbation theory c eff QQ = c QQ .

JHEP12(2017)044
One-loop corrections to the ALP-quark couplings c qq for both light and heavy quarks can be calculated in analogy with those to the ALP-lepton couplings discussed in section 3.2. The obvious replacements to be applied to relation (3.13) are Q → Q q and T 3 → T q 3 . In addition, the W -boson contribution picks up a factor V ik V * jk or V * ki V kj (summed over k) for external up-type or down-type quarks with generation indices i and j, respectively. If the internal quark with index k is heavy, a non-trivial loop function arises. Note that these contributions can be off-diagonal in generation space. Finally, there is a new one-loop contribution involving the ALP-gluon coupling, whose form is with C F = 4/3. The perturbative calculation of this expression can be trusted as long as m a Λ QCD and m q Λ QCD . For the light quarks, the appropriate infrared scale is not the quark mass but a typical hadronic scale such as m π . We have derived the estimate (3.5) by using the above result for the gluon contribution to c qq in (3.2).

Summary of ALP decay modes
Above we have presented an overview of possible ALP decay modes into SM particles. The upper panel in figure 3 shows the various decay rates for a new-physics scale Λ = 1 TeV as a function of the ALP mass, under the assumption that the relevant coefficients |C eff γγ |, |C eff GG | and |c eff f f | are all equal to 1. For different values of these parameters, the rates need to be rescaled by factors (|C eff ii |/Λ) 2 . For example, in the lower panel we assume that the ALP-boson couplings are loop suppressed. If all Wilson coefficients are of the same magnitude and the ALP is lighter than the pion (or if it does not couple to colored particles at all), the dominant decay mode is a → γγ. The leptonic modes a → + − are only significant near the thresholds m a 2m , where they can be dominant. If the ALPboson couplings are loop suppressed, the leptonic decays can be dominant for ALP masses exceeding 2m e . The picture changes significantly if the ALP is heavy enough to decay hadronically, i.e. for m a > 3m π 0 ≈ 405 MeV. If the coupling to gluons is unsuppressed, the ALP then decays predominantly into hadronic final states. For m a > few GeV, the inclusive hadronic rate is approximately given by (3.18). If, on the other hand, the ALPgluon coupling is suppressed, there can be a potpourri of decay modes (a → hadrons, a → bb, a → cc, a → τ + τ − , a → γγ) with potentially similar rates. Which of these modes dominates depends on the details of the model.
If the total decay rate of the ALP is too small, the ALP leaves the detector before it decays. For example, a total rate of 10 −9 eV corresponds to a lifetime of 6.6 · 10 −7 s. If the ALP is produced in decays of heavier particles, the Lorentz boost can increase its lifetime significantly. It is also a possibility that the ALP decays invisibly into light particles of a hidden sector. In this case the decay products cannot be reconstructed, and hence the ALP signature would be that of missing energy and momentum.

Constraints on ALP couplings to photons and electrons
The couplings of ALPs to photons and electrons have been constrained over vast regions of parameter space using a variety of experiments in particle physics, astro-particle physics and cosmology. Since our work is motivated by the idea that ALPs could interact with SM particles with couplings of order (1 TeV) −1 to (100 TeV) −1 , such that these interactions can be probed at the LHC, we need to address the question of how the existing bounds can be satisfied. In figure 4 we show a compilation of existing exclusion regions for the ALP-photon and ALP-electron couplings. Before addressing these bounds in more detail, let us add an important remark concerning the ALP-lepton couplings. In the absence of a flavor symmetry, under which the three lepton flavors carry different charges (but which must be broken in order to explain neutrino oscillations), the matrices C L and C e entering the ALP-lepton couplings in (2.1) must, to an excellent approximation, be proportional to the unit matrix. Otherwise it is impossible to avoid flavor-changing neutral currents in the charged lepton sector, which are generated after electroweak symmetry breaking, see (2.9). The relevant couplings c eµ , c µτ and c eτ must satisfy very strong constraints from processes JHEP12(2017)044 such as µ → eγ and µ − → e − e + e − , and analogous ones involving heavier leptons (see [64] for a recent review). As a result, one expects that c ee c µµ c τ τ (3.21) to very good accuracy. Below we will sometimes make use of this relation.

Constraints on the ALP-photon coupling
Consider first the exclusion regions in the m a − |C eff γγ | plane shown in the left panel. The parameter space excluded from cosmological constraints is shaded gray. This includes constraints from measurements of the number of effective degrees of freedom, modifications to big-bang nucleosynthesis, distortions of the cosmic microwave-background spectrum and extragalactic background-light measurements [65,66]. Energy loss of stars through radiation of ALPs is constrained by the ratio of red giants to younger stars of the so-called horizontal branch (HB) [67][68][69] (shaded purple). Another strong constraint arises from the measurement of the length of the neutrino burst from Supernova SN1987a, which would have been shorter in the presence of an energy loss from ALP emission [70] (shaded yellow), as well as from the non-observation of a photon burst from SN1987a due to the decay of emitted ALPs [71] (shaded orange). These constraints require an extremely tight bound |C eff γγ |/Λ 10 −15 TeV −1 in the mass window between 150 eV and about 1 MeV. For smaller ALP masses the bounds are weaker, ranging from |C eff γγ |/Λ < 10 −9 TeV −1 for m a = 150 eV to |C eff γγ |/Λ < 3 · 10 −7 TeV −1 for m a < 4 eV. Below 4 eV the tightest bounds come from HB stars and axion helioscopes like the Tokyo Axion Helioscope (SUMICO) and the CERN Axion Solar Telescope (CAST), which search for ALPs produced in the Sun and exclude the blue parameter space [72][73][74]. Above the threshold m a = 2m e ≈ 1 MeV, decays of the ALPs into electron-positron pairs may affect the assumptions of some of these constraints in a non-trivial way. In the sub-eV mass range, light-shining-through-a-wall experiments (LSW) also provide interesting constraints.
Beam-dump searches are sensitive to ALPs radiated off photons, which are exchanged between the incoming beam and the target nuclei (Primakoff effect) and decay back to photons outside the target. The orange area is a compilation of different runs performed at SLAC [75,76]. Radiative decays Υ → γa of Upsilon mesons have been searched for at CLEO and BaBar [77,78], and yield the excluded area shaded light green. Bounds from collider searches for ALPs include searches for mono-photons with missing energy (e + e − → γa) at LEP (dark orange), tri-photon searches on and off the Z-pole (e + e − → 3γ) at LEP (light blue), and searches for the same final state at CDF (purple) and LHC (dark orange). A detailed discussion of these searches can be found in [23][24][25]. For ALP masses in the multi-GeV range, alternative searches for ALP production in ultra-peripheral heavyion collisions have the potential to improve the current bounds by up to two orders of magnitude, provided the a → γγ branching ratio is close to 100% [25]. First evidence for light-by-light scattering in 480 µb −1 of Pb-Pb collision data has recently been reported by ATLAS [79]. While the derivation of the precise bound on the ALP-photon coupling is beyond the scope of this work, the green area labeled "Pb" shows an estimate obtained based on a rescaling of the projected limit presented in [25] to the luminosity used in the JHEP12(2017)044 Figure 4. Existing constraints on the ALP-photon (left) and ALP-electron coupling (right) derived from a variety of particle physics, astro-particle physics and cosmological observations. Several of these bounds are model dependent. The BaBar constraint in the right-hand plot assumes c µµ ≈ c ee , see (3.21); otherwise, this is a bound on |c eff µµ |. See the text for more details.
ATLAS analysis. Beam-dump experiments and collider searches are directly sensitive to the presence of additional ALP couplings for masses m a > 2m e . The reach of beam-dump experiments, for example, would be strongly reduced if ALPs would decay into electrons before they leave the beam dump. The limits from collider searches and those derived fro heavy-ion collisions shown in the plot assume Br(a → γγ) = 1. The corresponding exclusion regions would move upwards if this assumption was relaxed. Also, in some cases specific assumptions about the relation between C γγ and C γZ were made, which have an influence on the results.
It follows from this discussion that the ALP-photon coupling is most severely constrained for all ALP masses below about 1 MeV. At tree-level, this requires that the combination C γγ − 1.92 C GG = C W W + C BB − 1.92 C GG of the Wilson coefficients of the operators in which the ALP couples to gauge fields in (2.1) must be extremely small, of order (10 −9 − 10 −7 ) (Λ/TeV) for m a < 150 eV, and less than 10 −15 (Λ/TeV) for 150 eV < m a < 1 MeV. If we assume that Λ lies within a few orders of magnitude of the TeV scale, these constraints would either require an extreme fine tuning or (better) a mechanism which enforces that C BB = −C W W and C GG = 0. (However, integrating out a single, complete electroweak multiplet will always generate contributions to C W W and C BB with same sign.) The assumption that such a cancellation can be engineered was made in the recent analysis in [26]. Moreover, relation (3.2) shows that even in this case an effective coupling C eff γγ = 0 will inevitably be generated at one-loop (and higher-loop) order as long as some couplings in the effective Lagrangian are set by the TeV scale. To JHEP12(2017)044 see this, consider the following numerical results in the relevant mass window: For ALP masses below 100 keV each loop contribution scales with m 2 a . We observe that reaching |C eff γγ |/Λ < 10 −15 TeV −1 requires a significant fine-tuning of essentially all Wilson coefficients in the effective Lagrangian (2.1). This includes the coefficient C W W , even though its one-loop contribution is very small. As we will show below, the one-loop radiative corrections to the ALP-electron coupling induce a contribution δc ee ≈ −0.8 · 10 −2 C W W independently of the ALP mass, which adds the terms 5 · 10 −5 C W W and 2 · 10 −7 C W W to the two values shown in (3.22). It follows that ALPs with masses in the range between 150 eV and 1 MeV are incompatible with the assumption of couplings to SM particles that could be probed at high-energy particle colliders. For masses below 150 eV, on the other hand, a mechanism which sets C γγ = 0 and C GG = 0 at tree level would be sufficient to satisfy the relevant constraints irrespective of the values of the remaining ALP couplings.
The left panel in figure 4 shows that above 30 MeV a window opens for |C eff γγ |/Λ ∼ 1 TeV −1 , and above 400 MeV the ALP-photon coupling is essentially unconstrained as long as it falls between (10 −6 − 10 +1 ) TeV −1 . The mass range m a > 30 MeV is thus the best motivated region to search for ALPs at high-energy particle colliders. It is interesting to study loop corrections also for this high-mass region. They can be sizable for all particles lighter than the ALP. For example, at m a = 10 GeV we find where we have included the two-loop estimate (3.5) for the contribution from C GG . In addition, we expect a two-loop contribution proportional to C W W of order 10 −2 inside the square bracket. Sizable loop corrections can be generated either if C GG = O(1), or if some of the fermion couplings are of O(10). For example, setting c ee = c µµ = c τ τ = 10 for the charged leptons would not lead to any tensions with perturbativity.

Constraints on the ALP-electron coupling
Exclusion limits on the ALP-electron coupling are shown on the right panel of figure 4. They include searches by the Edelweiss collaboration (shaded purple) [80] for ALPs produced in the Sun by the Compton process γe − → e − a, by bremsstrahlung e − X → Xa off JHEP12(2017)044 electrons or hydrogen and helium nuclei in the plasma, and by ALP radiation from excited ions. Even stronger limits for m a < 10 −5 GeV are derived from observations of Red Giants (shaded red). ALP radiation can lead to the cooling of the cores of these stars, which leads to delayed Helium ignition and modifies the brightness-temperature relation [69]. Axion radiation from electron beams is further constrained by beam-dump experiments performed at SLAC (shaded blue) [81]. The presence of a sizable ALP-photon coupling would reduce the reach of beam-dump experiments and could affect the astrophysical constraints in a non-trivial way. In particular, the Edelweiss bounds assume that ALPs produced in the Sun do not decay on their way to Earth, which would require that the ALP-photon coupling is tuned to zero with high precision, which is rather implausible in view of our discussion in the previous section. We note, however, that a viable scenario can be obtained by setting the tree-level ALP couplings to quarks and gauge bosons to zero. An ALP-photon coupling is then still induced at one-loop order, see (3.11), and for m a < m e it is to good approximation given by C eff Requiring that the average decay length of the ALP is larger than the Earth's distance to the Sun, we then obtain the bound The dashed and dotted lines intersecting the Edelweiss constraint in figure 4 indicate this bound for E a = 14 keV and 1 keV, respectively. Below these lines, ALPs with the corresponding minimum energies are sufficiently long-lived to travel from the Sun to the Earth before decaying. We also note that limits on the ALP-electron coupling in the mass range between 20 MeV and 10 GeV can be derived from dark-photon searches performed at MAMI [82] and BaBar [83]. While a proper conversion of these limits is non-trivial [84] and beyond the scope of this work, the bounds one obtains are typically rather weak, of order |c eff ee |/Λ 10 3 TeV −1 . Assuming the approximate universality of the ALP-lepton couplings shown in (3.21), a stronger constraint can be derived from a dark-photon search in the channel e + e − → µ + µ − Z performed by BaBar [85], which we will reanalyze in the context of our model in the next section. For C γγ = 0, this gives rise to the bound shaded in gray in figure 4.
Of the one-loop contributions to the effective ALP-electron coupling in (3.13), only the photon term shows a sizable sensitivity to the ALP mass, and only in the region where m a m e . We find (with µ = Λ = 1 TeV in the argument of the logarithms) To satisfy the model-independent bound |c eff ee |/Λ < 10 −6 TeV −1 in the mass range m a < 10 keV would require that |C γγ | and |C W W | (and hence both |C W W | and |C BB |) must be smaller than approximately 10 −4 (Λ/TeV) in this low-mass region.

Anomalous magnetic moment of the muon
The persistent deviation of the measured value of the muon anomalous magnetic moment a µ = (g − 2) µ /2 [86]  physics. The difference a exp µ − a SM µ = (29.3 ± 7.6) · 10 −10 , where we have taken an average of two recent determinations [87,88], differs from zero by about 4 standard deviations. It has been emphasized recently that this discrepancy can be accounted for by an ALP with an enhanced coupling to photons [11]. At one-loop order, the effective Lagrangian gives rise to the contributions to a µ shown in figure 5. The first graph, in which the ALP couples to the muon line, gives a contribution of the wrong sign [89,90]; however, its effect may be overcome by the second diagram, which involves the ALP coupling to photons (or to γZ), if the Wilson coefficient C γγ in (2.1) is sufficiently large [10,11]. Performing a complete one-loop analysis, we find that the effective ALP Lagrangian gives rise to the new-physics contribution (4.1) The loop functions read (with (4.2) They are positive and satisfy h 1,2 (0) = 1 as well as h 1 (x) ≈ (2/x)(ln x − 11 6 ) and h 2 (x) ≈ (ln x + 3 2 ) for x 1. The scheme-dependent constant δ 2 = −3 is again related to the treatment of the Levi-Civita symbol in d dimensions, see appendix C.
Note that in processes in which the ALP only appears in loops but not as an external particle, the scale dependence arising from the UV divergences of the ALP-induced loop contributions are canceled by the scale dependence of a Wilson coefficient in the D = 6 effective Lagrangian of the SM. In the present case the relevant term yielding a tree-level contribution to a µ reads (written in the broken phase of the electroweak theory) In order to calculate the Wilson coefficient K aµ one would need to consider a specific UV completion of the effective Lagrangian (2.1). The large logarithm in the term proportional to C γγ in (4.1) is, however, unaffected by this consideration. The coefficient we obtain for JHEP12(2017)044 Figure 6. Regions in ALP coupling space where the experimental value of (g − 2) µ is reproduced at 68% (red), 95% (orange) and 99% (yellow) confidence level (CL), for different values of m a . We assume K aµ (Λ) = 0 at Λ = 1 TeV and neglect the tiny contribution proportional to C γZ . For m a > 2m µ , the gray regions are excluded by a dark-photon search in the e + e − → µ + µ − + µ + µ − channel performed by BaBar [85]. this logarithm agrees with [11] (the remaining finite terms were not displayed in this reference). Two-loop light-by-light contributions proportional to (C γγ /Λ) 2 have been estimated in [11] and were found to be approximately given by For µ = Λ = 1 TeV this evaluates to δa µ | LbL ≈ 5.6 · 10 −12 C 2 γγ . In the region of parameter space we consider, where |C γγ |/Λ 2 TeV −1 (see below), the impact of this effect is tiny.
In our numerical analysis, we will assume that the contribution of K aµ (µ) is subleading at the high scale µ = Λ. If the Wilson coefficients c µµ and C γγ are of similar magnitude, the logarithmically enhanced contribution is the parametrically largest one-loop correction. It gives a positive shift of a µ provided the product c µµ C γγ is negative. The correction proportional to C γZ is suppressed by (1 − 4s 2 w ) and hence is numerically subdominant. Figure 7. Tree-level Feynman diagrams contributing to the process e + e − → µ + µ − a.

JHEP12(2017)044
Note also that the contribution proportional to (c µµ ) 2 is suppressed in the limit where m 2 a m 2 µ , while the remaining terms remain unsuppressed. Figure 6 shows the regions in the parameter space of the couplings c µµ and C γγ in which the experimental value of the muon anomalous magnetic moment can be explained in terms of the ALP-induced loop corrections shown in figure 5, without invoking a large contribution from the unknown short-distance coefficient K aµ (Λ). There is a weak dependence on the ALP mass, such that the allowed parameter space increases for m 2 a m 2 µ . Interestingly, we find that an explanation of the anomaly is possible without much tuning as long as one coefficients is of order Λ/TeV, while the other one can be of similar order or larger. Since c µµ enters observables always in combination with m µ , it is less constrained by perturbativity than C γγ .
An important constraint on the ALP-photon and ALP-muon couplings, C γγ and c µµ , can be derived from a search for light Z bosons performed by BaBar, which constrains the resonant production of muon pairs in the process The Feynman diagrams contributing to this process at tree level (and for m e = 0) are shown in figure 7. Neglecting the electron mass and averaging over the initial-state polarizations, we obtain for the cross section where r = m 2 a /s and = m 2 µ /s are dimensionless ratios, and √ s ≈ 10.58 GeV is the centerof-mass energy. Note that the contributions involving the ALP-muon coupling are chirally suppressed by a factor = m 2 µ /s and hence are numerically very small in the region where C γγ and c µµ take values of similar magnitude. The contributions involving the ALP-photon coupling are logarithmically divergent in the limit m µ → 0. Neglecting terms of O( ) and higher in the coefficient functions, which is an excellent approximation numerically, we find (4.6)

JHEP12(2017)044
In order to compute the resonant e + e − → µ + µ − a → µ + µ − µ + µ − cross section, we need to multiply expression (4.5) with the a → µ + µ − branching ratio. Assuming that only the Wilson coefficients C γγ and c µµ are non-zero, and that the ALP couplings to charged leptons are flavor universal, we obtain (for m a > 2m µ ) where the sum in the denominator extends over all lepton flavors with 2m < m a . If additional decay channels were present, the bounds derived below would become weaker. At one-loop order, the effective ALP-photon coupling receives contributions proportional to c µµ , which have been shown in (3.2) and (3.11). These loop-induced effects contribute to (4.5) at a level comparable to the chirally-suppressed tree-level contributions involving c µµ . In order to properly account for the full dependence on c µµ , one should thus use the effective ALP-photon coupling instead of C γγ in (4.5) and (4.7). For a given value of the ALP mass in the range 2m µ < m a < √ s − 2m µ the product is bounded from above by the values shown in figure 4 of [85]. In applying these bounds, we perform an average over the mass range [m a − 0.5 GeV, m a + 0.5 GeV] to smooth out the spiky structures seen in the figure. The resulting exclusion regions in the c µµ − C γγ plane arising at 90% CL are shown by the gray regions in figure 6. In the mass range just above the di-muon threshold, the exclusion region derived from the BaBar analysis lies close to the region where (g − 2) µ can be explained and indeed excludes a small portion of this region. On the other hand, for ALP masses below 2m µ no constraints arise, and for m a > 1.5 GeV the constraints quickly become rather weak. We emphasize, however, that ALP searches at the upcoming Belle II super flavor factory, both in the a → µ + µ − and a → γγ channels, have the potential to significantly tighten these constraints and exclude an ALP-based explanation of the muon anomaly in the mass range from 2m µ up to a few GeV.

Exotic decays of the Higgs boson into ALPs
The presence of ALP couplings to SM particles gives rise to the possibility of various exotic decay modes of the Higgs boson, which might be discoverable during the high-luminosity run of the LHC. The relevant decay modes are h → Za and h → aa. These offer a variety of interesting search channels for ALPs, depending on how the ALP and the Z boson JHEP12(2017)044 decay. In some regions of parameter space, the decay h → Za may be reconstructed in the h → Zγ search channel and appear as a new-physics contribution to this decay mode. The present experimental upper limits on the pp → h → Zγ rates reported by CMS [91] and ATLAS [92] (both at 95% confidence level (CL)) are 9 and 11 times above the SM value, respectively, thus leaving plenty of room for new-physics effects. A discovery of the h → Zγ decay mode and an accurate measurement of its rate are among the most pressing targets for the high-luminosity LHC run. Very importantly, we will show that ALP searches in the h → Za and h → aa channels with subsequent a → γγ or a → e + e − decays can potentially probe regions in the m a -C eff γγ and m a -c eff ee parameter spaces that are inaccessible to any other searches.
The lifetime of ALPs and their boost factor have important consequences for their detectability. For very light ALPs or very weak couplings, the decay length can become macroscopic and hence only a small fraction of ALPs decay inside the detector. Since to good approximation Higgs bosons at the LHC are produced along the beam direction, the average decay length of the ALP perpendicular to the beam axis is where θ is the angle of the ALP with respect to the beam axis, β a and γ a are the usual relativistic factors, and Γ a is the total decay width of the ALP. For the example of h → Za decay followed by a → γγ, the geometry is sketched in figure 8. Note that the quantity L ⊥ a (θ) (but not L a ) is invariant under longitudinal boosts along the beam axis, and we are thus free to define L a and the angle θ in the Higgs-boson rest frame. If the ALP is observed in the decay mode a → XX, we can express its total width in terms of the branching fraction and partial width for this decay, yielding

JHEP12(2017)044
We call f Za dec and f aa dec the fraction of all h → Za and h → aa events where the ALPs decay before they have traveled a perpendicular distance L det set by the relevant detector components needed for the reconstruction of the particles X (i.e., the electromagnetic calorimeter if X is a photon, and the inner tracker if X is an electron). Since two-body decays of the Higgs boson are isotropic in the Higgs rest frame, it follows that These integrals are discussed in more detail in appendix D. Both event fractions are exponentially close to 1 if L a L det . Numerically, one finds that f aa dec ≈ (f Za dec ) 2 to very good approximation, unless the ratio L det /L a 1. In the latter case one obtains We now define the effective branching ratios where Br(Z → + − ) = 0.0673 for = e, µ. If the decay length L a L det , the effective branching ratios are just the products of the relevant branching fractions for the individual decays. They depend on the squares of the Wilson coefficients C eff Zh and C eff ah , which govern the Higgs decay rates into ALPs, and on the branching ratio Br(a → XX) for the decay mode in which the ALP is reconstructed. In the opposite case, where the ALP decay length is larger than the detector scale L det , the dependence on the a → XX branching ratio drops out to good approximation, because the relevant product Br(a → XX)/L a ∝ Γ(a → XX) is governed by the a → XX partial decay rate. Via this rate enters a dependence on the Wilson coefficient C eff XX responsible for the decay a → XX. This behavior is illustrated in figure 9, which shows the effective branching ratio Br(h → Za → + − γγ) eff for different values of the a → γγ branching ratio and the relevant coefficient C eff γγ mediating the di-photon decay. We keep the h → Za branching fraction fixed at 10% for m a = 1 GeV. The two solid curves correspond to fixed Br(a → γγ) = 1 along with |C eff γγ |/Λ = 1/TeV (blue) and |C eff γγ |/Λ = 0.1/TeV (red). For sufficiently large ALP mass the same asymptotic value for the effective branching ratio is obtained, but the reach towards low masses depends sensitively on the value of C eff γγ . The two dotted lines are obtained in the same way, but with Br(a → γγ) = 0.1. In this case the asymptotic value for the effective branching ratio is reduced by a factor 10, but the behavior in the low-mass region is the same as before. As explained above, for low masses the effective branching ratio becomes independent of Br(a → γγ), while for large masses it becomes independent of C eff γγ . Figure 9. Effective h → Za → + − γγ branching ratio as a function of the ALP mass for a fixed value Br(h → Za) = 0.1 at m a = 1 GeV. The solid lines refer to a 100% a → γγ branching ratio along with |C eff γγ |/Λ = 1/TeV (blue) and |C eff γγ |/Λ = 0.1/TeV (red). The dotted lines are obtained by lowering the a → γγ branching ratio to 10%.

ALP searches in h → Za decay
The relevant Feynman diagrams contributing to the h → Za decay amplitude up to oneloop order are depicted in figure 10. The effective Lagrangian (2.1) does not contain a dimension-5 operator contributing to the h → Za decay amplitude at tree level. The only contribution arising at this order is due to fermion loop graphs. Because both the Higgs boson and the ALP couple to fermions proportional to the fermion mass, the only relevant effects comes from the top quark. The W -boson loop diagram shown in the second graph vanishes, since there are not enough 4-vectors available to saturate the indices of the Levi-Civita tensor in the aW W vertex. A tree-level contribution to the h → Za decay amplitude (third graph) arises first at dimension-7 order, from the third operator shown in (2.6). Evaluating all contributions, we obtain [47,48] where λ(x, y) = (1 − x − y) 2 − 4xy, and we have defined Zh .
Here y t and T t 3 = 1 2 are the top-quark Yukawa coupling and weak isospin, and C Zh = 0. The top-quark contribution involves the parameter integral Zh − 0.016 c tt + 0.030 C  Figure 10. Feynman diagrams contributing to the decay h → Za.
The left plot in figure 11 shows our predictions for the h → Za decay rate normalized to the SM rate Γ(h → Zγ) SM = 6.32 · 10 −6 GeV [93]. We set C Zh = 0 and display the rate ratio in the plane of the Wilson coefficients c tt and C Zh . Since only the relative sign of the two coefficients matters, we take C (7) Zh to be positive without loss of generality. We find that, in a large portion of parameter space, the exotic h → Za mode can naturally have a similar decay rate as the h → Zγ mode in the SM, especially if the top-quark contribution interferes constructively with the dimension-7 contribution proportional to C Zh . The argument for the absence of a tree-level dimension-5 contribution to the h → Za decay amplitude holds in all new-physics models, in which the operators in the effective Lagrangian arise from integrating out heavy particles whose mass remains large in the limit of unbroken electroweak symmetry [47,48]. However, this argument does not apply for the class of models featuring new heavy particles which receive all or most of their mass from electroweak symmetry breaking. Concrete examples of such models include little-Higgs models, in which fermionic top partners can have very large Higgs couplings [94,95], and triplet-doublet dark matter models with vector-like leptons [96,97], which are generalizations of the Wino-Higgsino dark matter scenario in the minimal supersymmetric standard model. The effective Lagrangian for such models generically contains operators which are non-polynomial in the Higgs field (see e.g. [98]). At dimension-5 order, there is a unique such operator relevant to the decay h → Za. It is given by [47,48] L non−pol Its contribution to the decay amplitude was already included in (5.7) and (5.9). The decay h → Za is unique in the sense that, at dimension-5 order, a tree-level hZa coupling can only arise in such special models. Note that the non-polynomial operator in (5.10) also arises at one-loop order in the SM. Integrating out the top-quark from the effective Lagrangian generates a contribution to C Zh given by the second term in (5.7) evaluated with F = 1. In the right plot in figure 11, we allow for non-zero C (5) Zh and display the rate ratio as a function of the effective Wilson coefficient C eff Zh defined in (5.9) for different ALP masses. In models where a tree-level dimension-5 contribution is present, one can naturally obtain h → Za rates exceeding the SM h → Zγ rate by orders of magnitude. For example, with |C eff Zh |/Λ = 0.3 TeV −1 and for a light ALP (m a < 1 GeV) one finds a ratio of about 60, corresponding to a 9% h → Za branching ratio. This would be a spectacular new-physics effect. We find that the decay rate is approximately independent of the ALP mass as long as m a is below a few GeV. The decay h → Za is kinematically allowed as long as m a < m h − m Z ≈ 33.9 GeV. Figure 11 shows that significant decay rates can be found even close to the kinematic limit.  Depending on the dominant branching ratio of the ALP, the decay h → Za can give rise to various interesting experimental signatures. ALP decays into photons can be searched for in the h → Za → + − γγ final state. No dedicated searches have been performed in this channel yet. However, for strongly boosted ALPs the two photons would be reconstructed as a single photon jet, and the decays h → Za would then lead to a modification of the observed pp → h → Zγ rate. Since there is no interference term, this rate would necessarily be enhanced in this case. From figure 11 it follows that this enhancement can easily be of O(1) and stronger. We estimate the mass below which a di-photon decay of the ALP will mimic a single photon in the detector to be about 47 MeV by following the analysis for h → aa decay of [103] and accounting for the different Lorentz boost factors (see the discussion in section 5.2). The current best limit on the cross section of JHEP12(2017)044  [91] then rules out the shaded area above the solid and dotted blue lines in the left panel of figure 12. The lines in this figure have the same meaning as in figure 9. Solid and dotted lines refer to Br(a → γγ) = 1 and Br(a → γγ) = 0.1, respectively. Blue lines are obtained with |C eff γγ |/Λ = 1/TeV, while red lines correspond to |C eff γγ |/Λ = 0.1/TeV. With present luminosity, only the former choice gives rise to non-trivial bounds. As explained above, for low ALP masses the constraints become independent of the a → γγ branching ratio. For very low ALP masses sensitivity is lost, because most of the ALPs decay outside the detector.
If the leptonic decay modes are relevant, ALPs can be searched for in h → Za → 4 decays. An analysis by ATLAS searching for new "dark" bosons Z d produced in Higgs decays h → ZZ d with subsequent decays ZZ d → 4 , where = e or µ, can be reinterpreted to constrain C eff Zh in the considered mass window m Z d = (15 − 35) GeV [44]. We show the excluded region in the left panel of figure 13, in which the solid and dotted contours correspond to Br(a → + − ) = 1 and 0.1, respectively. For these high ALP masses, the h → Za → 4 rate is essentially independent of the values of the Wilson coefficients |c eff |. We strongly encourage our experimental colleagues to extend these searches to lower masses and to separate the final-state lepton flavors. The expected asymmetry between electron, muon and tau final states from ALP decays would be a striking signature of a light pseudoscalar boson. The possibility to observe light new particles in Higgs decays with this final state has also been pointed out in [104]. A heavier ALP can also decay into heavy-quark pairs, which would provide spectacular signatures such as h → Za → + − bb, or into di-jets, i.e. h → Za → + − j(j), where a single jet would be observed in the case JHEP12(2017)044  of two strongly collimated jets. Very light or weakly coupled ALPs can remain stable on detector scales. In this case, a Higgs produced in vector-boson fusion or in association with a Z-boson or a top-quark pair can lead to interesting signatures of the type pp → hjj → Z + / E T + jj, pp → hZ → ZZ + / E T , or pp → htt → Z + / E T + tt.

ALP searches in h → aa decay
By means of the Higgs portal interactions in the dimension-6 effective Lagrangian (2.6), as well as by loop-mediated dimension-6 processes, a Higgs boson can decay into a pair of ALPs. We have calculated the h → aa decay rate including the tree-level Higgs-portal interactions as well as all one-loop corrections arising from two insertions of operators from the dimension-5 effective Lagrangian (2.1). The relevant diagrams are shown in figure 14.
Since both the Higgs boson and the ALP couple to fermions proportional to their mass, only the top-quark contribution needs to be retained in the second diagram. Keeping m a only in the phase space and neglecting it everywhere else, we find where the effective coupling is given by Note that the second Higgs-portal interaction in (2.6) does not contribute in this approximation, because its effect is suppressed by m 2 a /m 2 h . Numerically, we obtain for Λ = 1 TeV indicating that the top-quark contribution, in particular, can be sizable. Relation (5.13) shows that even if the portal coupling C ah vanishes at some scale, an effective coupling is induced at one-loop order if the ALP couples to at least one of the heavy SM particles (t, Z or W ). Also, because of the presence of UV divergences in the various terms, the coupling C ah (µ) must cancel the scale dependence of the various other terms, and hence it is not consistent to set it to zero in general. For a light ALP (m a < 1 GeV) a 10% h → aa branching ratio is obtained for |C eff ah |/Λ 2 = 0.62 TeV −2 . Note that a Wilson coefficient of this size could even be due to a loop-induced contribution from the top quark, if |c tt |/Λ ≈ 1.9 TeV −1 .

JHEP12(2017)044
Imposing the current upper limit Br(h → BSM) < 0.34 (at 95% CL) [99], we obtain (5. 16) More generally, if both coefficients are non-zero, the allowed values for C eff Zh and C eff ah are constrained to lie within the orange region in figure 15. At the end of LHC operation, with a projected integrated luminosity of 3000 fb −1 at √ s = 14 TeV, one expects the improved bound Br(h → BSM) < 0.1 [100], which would imply that the two coefficients must be inside the dashed black contour in the figure. The constraint on C eff ah alone would then be |C eff ah | < 0.62 (Λ/TeV) 2 . Invisible ALP decays would lead to invisible Higgs-boson decays, for which the bounds Br(h → invisible) < 0.23 from ATLAS [101] and Br(h → invisible) < 0.24 from CMS [102] imply the constraint |C eff ah | < 1.02 (Λ/TeV) 2 for Br(a → invisible) = 1. Depending on the pattern of ALP decay modes, promising signals arise from multiphoton and multi-lepton final states, but also from ALP decays into jets or b quarks. Very light ALPs can only decay into photons and are boosted along the beam direction with a boost factor γ a = m h /(2m a ) 1, for which the photons are highly collimated. For masses m a < 625 MeV, the opening angle between the final state photons ∆φ = arccos(1 − 2/γ 2 a ) ≈ 2/γ a is smaller than the angular resolution of the ATLAS and CMS electromagnetic calorimeters (ECALs) of ∼ 20 mrad, and hence the photons enter the same calorimeter cell [31,33,105]. However, shower-shape analyses allow one to differentiate between single and multiple photons even if the opening angle is below the angular resolution. To be conservative, and based on the analysis in [103], we therefore assume that ALP masses below 100 MeV cannot be distinguished from h → γγ decays. In this case, we can turn the limit on the signal strength parameter µ h→γγ exp = 1.14 + 0. 19 − 0.18 [99] into a constraint on the h → aa → γγ + γγ branching ratio, where the effective Higgs branching ratio Br(h → aa → γγ + γγ) eff is defined as in (5.5) and takes into account the lifetime of the ALPs. This constraint is shown by the contours in the low-mass region of the right panel of figure 12, where the meaning of the curves is the same as in figure 9. The solid and dotted curves correspond to Br(a → γγ) = 1 and Br(a → γγ) = 0.1, respectively, while the blue and red curves refer to |C eff γγ |/Λ = 1 TeV −1 and |C eff γγ |/Λ = 0.1 TeV −1 . ATLAS further provides limits on Br(h → aa → γγ + γγ) for the three mass values m a = 100 MeV, m a = 200 MeV and m a = 400 MeV, based on the √ s = 7 TeV dataset [103]. The corresponding limits are indicated by the three blue or red points in the figure. For ALP masses in the range m a = (10 − 62.5) GeV, ATLAS has performed a dedicated search for h → aa → 4γ [39]. We show the corresponding bounds in the right panel of figure 12. In this case the red contours overlap with the blue ones, since the value of C eff γγ becomes irrelevant as long as the a → γγ branching ratio takes a fixed value. It is apparent that the limits for very light ALP masses are independent of the choice of Br(a → γγ), while the limits for heavy ALPs are unchanged for smaller Wilson coefficients C eff γγ , as expected from the discussion of figure 9.

Probing the parameter space of ALPs
Given the rich phenomenology of ALP decays, there is a plethora of promising searches at the LHC for both h → Za and h → aa decays. If the a → γγ branching ratio is sufficiently large, these exotic Higgs decays with subsequent ALP decays into photons would give rise to very clean signatures, which can be used to discover or constrain the ALP-photon coupling in a vast region of so far unexplored parameter space [46]. Equally interesting are ALP decays into lepton pairs, which would also lead to clean final states. We now discuss the prospects for searches in these two channels and present projections for the reach of Run-2 of the LHC. ALP decays into hadronic final states, such as dijets or heavy QQ pairs, are experimentally more challenging and would require dedicated analyses. We emphasize that our focus in this work is on visibly decaying ALPs, which can be reconstructed in the detector. Searches for invisibly decaying ALPs can be performed using the missing-energy signature in mono-X final states such as pp → Z

Constraining the ALP-photon coupling
Present and future searches for h → γγ + γγ and h → + − + γγ decays at the LHC can probe a large range of ALP-photon couplings. In our estimates below we focus on Run-2 of the LHC, which will provide an integrated luminosity of 300 fb −1 at √ s = 13 TeV. We require 100 signal events in each search channel and require that the ALPs decay before the electromagnetic calorimeter, which is typically located at a distance of approximately 1.5 m from the beam axis. We assume that the Higgs bosons are produced in gluon fusion with a cross section of σ 13 TeV (gg → h) = 48.52 pb [108]. Projections for higher luminosity (3 ab −1 at √ s = 14 TeV) and for a 100 TeV proton-proton collider will be presented elsewhere [57].
In our analysis we consider very different experimental searches. Light ALPs can effectively enhance the h → γγ branching ratio, heavier ALPs produce clearly separated di-photon resonances in h → aa → γγ + γγ decays, and ALPs with very small couplings can lead to displaced vertices. Experimental strategies to isolate the signal and suppress the background differ significantly for these searches. We are not in a position to provide detailed estimates of detector and reconstruction efficiencies, or to perform solid background estimates. Nevertheless, we believe that our requirement of 100 signal events in the respective search channels is realistic. For comparison, we note that the current precision of the h → γγ rate measurements excludes more than 340 new-physics events in this channel [99], the upper limit on h → Zγ decay allows for 400 new-physics events [91], and the search for h → aa for heavy ALPs [39] is sensitive to 120-390 events depending on the ALP mass (all at 95% CL). 2 Note that in the present work we do not make use of displaced-vertex signatures, which will help to greatly reduce the background in the region of parameter space where only a small fraction of the ALPs decays inside the detector. We hope that our analysis will trigger sufficient interest in the experimental community that dedicated analysis strategies will be developed by the experimental collaborations themselves. We begin by presenting the projected reach of searches for the decay h → Za → + − + γγ, for which the effective branching ratio has been defined in the first line of (5.5). In this case we require that

JHEP12(2017)044
The green shaded regions in the left panel of figure 16 show the parameter space which can be probed in Run-2 for different values of the relevant Wilson ALP-Higgs coupling. The three lines limiting these regions correspond to |C eff Zh |/Λ = 0.72 TeV −1 (solid contour), 0.1 TeV −1 (dashed contour) and 0.015 TeV −1 (dotted contour), taking into account the model-independent upper bound from h → BSM derived in (5.11). Note that the dotted line roughly corresponds to a TeV-scale coupling suppressed by a loop factor. With JHEP12(2017)044 300 fb −1 of luminosity it is possible to extend the search to slightly smaller couplings, but reaching sensitivity to couplings smaller than |C eff Zh |/Λ < 0.01 TeV −1 would require a larger luminosity. To draw the contours in the figure we have assumed that Br(a → γγ) = 1; however, it is important to realize that their shape is essentially independent of the value of the a → γγ branching ratio as long as this quantity is larger than a certain critical value, which is set by the required number of signal events (and as long as the ALP mass is not too close to the kinematic limit). These limiting values are Br(a → γγ) > 3 · 10 −4 (solid), 0.011 (dashed) and 0.46 (solid). Importantly, it is thus possible to probe the ALP-photon coupling even if the ALP predominantly decays into other final states. The triangular shape of the region of the projected reach is a consequence of the fact that ALPs with either small masses or small couplings, which fall beyond the left boundary of the region of sensitivity, live long enough (on average) to leave the detector. As discussed in section 5, the line in the m a − |C eff γγ | plane where this happens only depends on the partial width Γ(a → γγ) ∝ m 3 a |C eff γγ | 2 /Λ 2 , but not on Br(a → γγ). This argument only breaks down near the kinematic boundary m a = m h − m Z , where the h → Za decay rate becomes sensitive to the ALP mass. This behavior can also be understood from figure 9. Note that the region in parameter space that can be probed using exotic Higgs decays into ALPs almost perfectly complements the regions covered by existing searches. This will also be true for the other search channels discussed below. Whereas existing searches probe signatures of long-lived ALPs, in our case the ALPs are so short lived that their decays can be reconstructed in the detector. The red band in figure 16 shows the parameter space in which the anomalous magnetic moment of the muon can be explained in terms of loop corrections involving a virtual ALP exchange, assuming |C eff γγ |/Λ ≤ |c µµ |/Λ ≤ 5 TeV −1 . The upper bound on |c µµ | ensures that there is a substantial a → γγ branching ratio everywhere inside the red band. Notice that almost this entire parameter space can be covered by searches for exotic Higgs decays, provided that the Higgs-ALP coupling C Zh is sufficiently large. In the right panel of figure 16 we present the parameter space already excluded by present analyses placing upper bounds on the h → Zγ branching ratio [91,92]. These bounds apply in the low-mass region, where the two photons produced in the decay of the ALP are seen as a single photon jet in the calorimeter. The excluded parameter space shaded in dark green is obtained assuming |C eff Zh |/Λ = 0.72 TeV −1 and Br(a → γγ) > 0.04. In figure 17 we present the projected reach of searches for the decay h → aa → γγ +γγ, for which the effective branching ratio has been defined in the second line of (5.5). As previously, we require that The lines limiting the green shaded regions in the left panel correspond to |C eff ah |/Λ 2 = 1 TeV −2 (solid), 0.1 TeV −2 (dashed) and 0.01 TeV −2 (dotted), where the last value corresponds to a TeV-scale coefficient times a loop factor. We have used Br(a → γγ) = 1 in the plot, but once again the contours are essentially independent of the a → γγ branching ratio except for ALP masses close to the kinematic limit m a = m h /2. The corresponding limiting a → γγ branching ratios are Br(a → γγ) > 0.006, 0.049 and 0.49, respectively. Figure 17. Constraints on the ALP mass and coupling to photons derived from various experiments (colored areas without boundaries, adapted from [24]) along with the parameter regions that can be probed using the Higgs decays h → aa → 4γ. The left panel shows the reach of LHC Run-2 with 300 fb −1 of integrated luminosity (shaded in light green). We require at least 100 signal events. The contours correspond to |C eff ah |/Λ 2 = 1 TeV −2 (solid), 0.1 TeV −2 (dashed) and 0.01 TeV −2 (dotted). The red band shows the preferred parameter space where the (g − 2) µ anomaly can be explained at 95% CL. The right panel shows the regions excluded by existing searches for h → γγ and h → 4γ (shaded in dark green), where we assume |C eff ah |/Λ 2 = 1 TeV −2 .

JHEP12(2017)044
With 300 fb −1 of luminosity it is possible to extend the search to slightly smaller couplings, but reaching sensitivity to couplings smaller than |C eff ah |/Λ 2 < 0.005 TeV −2 would require larger luminosity. In the right panel of figure 17 we show the exclusion regions derived from the experimental searches presented in the right panel of figure 12, now projected into the m a − |C eff γγ | plane. We assume |C eff ah |/Λ 2 = 1 TeV −2 . These bounds are valid for branching ratios Br(a → γγ) > 0.07, 0.57, and 0.04 for the cases of the low-mass region below 100 MeV, the mass range between 100 and 400 MeV, and the high-mass region, respectively. They are obtained from the absence of a significant enhancement of the h → γγ rate [99], the search for h → γγ + γγ for intermediate masses [103], and the corresponding search in the high-mass region [39]. The fact that the exclusion region obtained in the low-mass region with a luminosity of 25 fb −1 per experiment is not much weaker than our projection for 300 fb −1 shown by the solid line in the left panel indicates that our requirement of 100 signal events is not unreasonable.
While the graphical displays in figures 16 and 17 correctly represent the regions in the m a − |C eff γγ | parameter space which can be probed using exotic Higgs decays, it is important to emphasize that finding a signal in these search regions will require sufficiently large ALP-Higgs couplings, as indicated by the solid, dashed and dotted contour lines in the plots. Consequently, not finding a signal in any of these searches would not necessarily exclude the existence of an ALP in this parameter space. An alternative way to present our results, which makes this fact more explicit, is shown in figure 18 for h → Za (upper panel)

JHEP12(2017)044
m a = 10 GeV m a = 1 GeV m a = 100 MeV  some reason the couplings arise at tree level. The red line corresponds to a model in which C eff γγ , C eff Zh and C eff ah are generated from one-loop diagrams involving the three SM up-type quarks, which are assumed to have equal couplings c uu = c cc = c tt . The orange dashed line corresponds to a model in which only the top-quark coupling c tt is non-zero. This provides a concrete example of a scenario in which the loop-induced ALP-Higgs couplings can be sizable, while the induced ALP-photon coupling tends to be very small. In each case, the relevant coupling |c tt |/Λ is varied between 0.1 TeV −1 and 10 TeV −1 , as indicated by the labels along the line. The a → γγ branching ratios obtained in these scenarios are 7 · 10 −4 for m a = 10 GeV, 27% for m a = 1 GeV, and 100% for m a = 100 MeV. In the high-mass case (m a = 10 GeV), the di-jet final state a → 2 jets would provide for a more promising search channel.

Constraining the ALP-lepton couplings
The analysis of the previous section can be extended to any other decay mode of the ALP. As a second example we consider the decays a → + − , which are kinematically accessible if m a > 2m . We stress that analogous analyses to the ones presented here could (and should) be performed for all other possible ALP decay modes.
The a → e + e − decay mode is of particular interest, since in the sub-MeV region the ALP-electron coupling has been constrained using a variety of experimental searches, as discussed in section 3.5.2. Using exotic Higgs decays, it will be possible to probe the ALP-electron coupling in the largely unexplored region above 1 MeV. The decay chains h → Za → + 1 − 1 + e + e − and h → aa → e + e − + e + e − provide clean search channels in this parameter space. The corresponding projections are shown by the green shaded regions in figure 19, where we require that (with 1 = e, µ) (5.20) respectively. In contrast to ALP decay into photons, we now set L det = 2 cm, since the ALP decay into electrons should take place before the inner tracker. The region of sensitivity is limited by contours obtained for different values of the relevant ALP-Higgs couplings. As before, these values are |C eff Zh |/Λ = 0.72 TeV −1 (solid), 0.1 TeV −1 (dashed) and 0.015 TeV −1 (dotted) for h → Za → + 1 − 1 + e + e − , and |C eff ah |/Λ 2 = 1 TeV −2 (solid), 0.1 TeV −2 (dashed) and 0.01 TeV −2 (dotted) for h → aa → e + e − + e + e − . We have used Br(a → e + e − ) = 1 for the green-shaded region in the plot, but as previously the contours are essentially independent of the a → e + e − branching ratio unless this quantity falls below certain threshold values, which are the same as before. For h → Za, one needs Br(a → e + e − ) > 3 · 10 −4 (solid), 0.011 (dashed) and 0.46 (dotted). For h → aa, one needs instead Br(a → e + e − ) > 0.006 (solid), 0.049 (dashed) and 0.49 (dotted). Similar to the case of ALP decays into photons, searches for rare Higgs decays have the potential to probe so far unconstrained parameter space.
The orange and red regions overlaid in the plots show, for comparison, the corresponding parameter space that can be covered in searches for the decay modes a → µ + µ − and JHEP12(2017)044 Figure 19. Constraints on the ALP mass and coupling to leptons derived from various experiments (colored areas without boundaries, adapted from [80,81]) along with the parameter regions that can be probed using the Higgs decays h → Za → + 1 − 1 e + e − (left) and h → aa → e + e − e + e − (right). The areas shaded in light green show the reach of LHC Run-2 with 300 fb −1 of integrated luminosity. We require at least 100 signal events. The contours in the left panel correspond to |C eff Zh |/Λ = 0.72 TeV −1 (solid), 0.1 TeV −1 (dashed) and 0.015 TeV −1 (dotted), while those in the right panel refer to |C eff ah |/Λ 2 = 1 TeV −2 (solid), 0.1 TeV −2 (dashed) and 0.01 TeV −2 (dotted). The orange and red regions overlaid in the plots show the corresponding parameter space that can be covered in searches for the decay modes a → µ + µ − and a → τ + τ − (see text for more explanations).
a → τ + τ − . For the latter mode, we have adopted the τ reconstruction efficiencies from the h → aa → τ + τ − + τ + τ − search performed by CMS in [37]. For each ALP, they require one tau lepton to decay into a muon and the second one to decay hadronically (with 60% reconstruction efficiency), leading to a rate reduction by a factor 0.13 for each ALP. The exclusion contours have been computed assuming Br(a → + − ) = 1 for both cases, but as previously the contours are essentially independent of the branching ratio unless this quantity falls below certain threshold values. For a → µ + µ − these are the same as for the electron case. For a → τ + τ − the limiting branching fractions are larger, due to the lower reconstruction efficiency. For h → Za, one needs Br(a → τ + τ − ) > 2 · 10 −3 (solid) and 0.008 (dashed). For h → aa, one needs instead Br(a → τ + τ − ) > 0.041 (solid) and 0.36 (dashed). We observe that the ALP-muon and ALP-tau couplings which can be probed are significantly smaller than the ALP-electron couplings. This simply reflects that the relevant decay rates scale with the square of the charged-lepton mass.
So far we have discussed searches in the a → e + e − channel independently of other leptonic ALP decay modes. We emphasize, however, that in many new-physics models one would expect a strong correlation between these modes. Indeed, if the leptonic couplings c are approximately flavor universal, as shown in (3.21), then the orange and red areas labeled µ + µ − and τ + τ − in figure 19 can actually be interpreted as parameter regions in which one can probe the ALP-electron coupling. Indeed, if the ALP is heavy enough to decay into muons or taus, the branching ratios for decays into lighter leptons will be JHEP12(2017)044 tiny, and it will only be possible to reconstruct the decay in the heaviest lepton that is kinematically allowed. Note that the combination of the three different search regions nicely complements the region covered by beam-dump searches.
Once again, it is instructive to consider an alternative way of representing the information contained in figure 19. For three different values of the ALP mass, the green-shaded areas to the right of the solid or dashed contours in figure 20 show the regions in the parameter space of the relevant ALP-Higgs and ALP-lepton couplings which can be probed in the exotic Higgs decays h → Za → + fermions have equal couplings c f f . The relevant leptonic branching ratios is this model are Br(a → e + e − ) ≈ 98% for m a = 100 MeV, Br(a → µ + µ − ) ≈ 100% for m a = 1 GeV, and Br(a → τ + τ − ) ≈ 7.5% for m a = 10 GeV.

Constraints from Z-pole measurements
The ALP couplings to electroweak gauge bosons can also be probed through precision measurements of the properties of Z bosons. As a concrete example, consider the production of a photon in association with an ALP in e + e − collisions. The relevant Born-level diagrams are shown in figure 21. Neglecting the electron mass, we find the cross section where √ s is the center-of-mass energy and θ denotes the scattering angle of the photon relative to the beam axis. ALP emission from the initial-state leptons vanishes in the limit m e = 0 and is otherwise strongly suppressed. The vector and axial-vector form factors are given by where Γ Z is the total width of the Z boson. If one makes the ad hoc assumption that the ALP only couples to photons, while C γZ = 0, then measurements of this cross section at LEP can be used to constrain the coupling C γγ [24]. However, in view of the general relations (2.8) this assumptions seems very artificial. Let us instead analyze the general structure of the cross section in more detail. At low energy (s m 2 Z ) the photon contribution dominates and produces a cross section (after integration over angles) At high energy (s m 2 Z ) one finds to good approximation where we have used that (1−4s 2 w ) ≈ 0 in the first term in the expression for V (s) in (6.2). By combining measurements of the cross sections at high and low energies it is thus possible to constraint the two coefficients C γγ and C γZ in a model-independent way. A much enhanced sensitivity to the aγZ coupling is obtained on the Z pole, where the cross section is given by subtracted by performing a scan about the peak position. In this way one obtains access to C γZ directly. This example nicely illustrates the main idea of our approach. By using on-shell decays of narrow, heavy SM particles into ALPs rather than the production of ALPs via an off-shell particle we obtain a much better sensitivity to the ALP couplings.
For the case of on-shell Higgs decays studied in [46] and in section 5 of the present work, the relevant enhancement factor is (m h /Γ h ) 2 ≈ 9.4 · 10 8 (assuming a SM Higgs width).
It has been pointed out in [26] that the Drell-Yan process pp → (γ/Z) * → γa at the LHC already provides better constraints on the ALP couplings than the corresponding process e + e − → (γ/Z) * → γa at LEP, which we have discussed above. An analogous statement applies for the on-shell decay, which we discuss in sections 6.1 and 6.2. Zpole measurements are also interesting in view of electroweak precision observables placing constraints on the Wilson coefficients C γγ and C γZ (or alternatively C W W and C BB ). These constraints are derived in section 6.3. Ultra-high precision studies of rare Z-boson decays could be performed at a future e + e − collider operating on the Z pole, which could provide samples of almost 10 12 Z bosons per year [109]. Projections for ALP searches at such a facility will be presented elsewhere [57].

ALP searches in Z → γa decay
The second operator in (2.7) induces the exotic Z-boson decay Z → γa at tree level. Including also the one-loop contributions from fermion loops, we obtain the decay rate where the effective Wilson coefficient C eff γZ is given by Here v f = 1 2 T f 3 − s 2 w Q f is the Z-boson vector coupling to fermion f , and we have defined the mass ratios τ f = 4m 2 f /m 2 a and τ f /Z = 4m 2 f /m 2 Z . The relevant loop function reads It obeys B 3 (τ f , τ f /Z ) ≈ 1 for all light fermions other than the top quark, for which B 3 (τ t , τ t/Z ) ≈ B 1 (τ t/Z ) ≈ −0.024 is very small. As in the case of the a → γγ decay discussed in section 3.1, the main effect of electroweak radiative corrections would be to JHEP12(2017)044 renormalize the gauge couplings. In the present case the coupling α associated with the photon is evaluated at q 2 = 0, while the coupling α(m Z )/(s 2 w c 2 w ) associated with the Z boson should be evaluated at q 2 = m 2 Z as indicated. The Z → γa branching fraction is obtained by dividing this partial decay rate by the Z-boson total width Γ Z . This yields Br(Z → γa) = 8.17 · 10 −4 C eff By requiring the Z-boson total width to agree with the direct measurement Γ Z = (2.495 ± 0.0023) GeV performed at LEP [110], an upper bound on the Wilson coefficient |C eff γZ | can be extracted. At 95% CL we find Br(Z → BSM) < 0.0018 and This bound is obtained by neglecting the ALP mass and gets weaker when m a approaches the kinematic threshold at m a = m Z .
To analyze the reach of this decay mode in probing the ALP-γZ, ALP-photon and ALP-electron couplings, we follow a similar strategy as discussed for Higgs decays in section 5. As before, the lifetime of the ALP is taken into account by defining the average decay length of the ALP perpendicular to the beam axis, L ⊥ a (θ) given in (5.1), where the relevant boost factor in the Z-boson rest frame is now β a γ a = (m 2 Z − m 2 a )/(2m a m Z ). The fraction f γa dec of all Z → γa events in which a decays before traveling a characteristic distance L det is given by the same expression as in the first line of (5.3). In analogy with (5.5), we define the effective branching ratio Br(Z → γa → γXX) eff = Br(Z → γa) Br(a → XX) f γa dec . (6.11)

JHEP12(2017)044
The ALP branching ratios determine which final states are the most interesting ones. ALP decays into photons lead to the experimental signature Z → γa → γγγ. Bounds on this branching ratio can be derived from precision studies of Z-boson decays performed at LEP, the Tevatron and the LHC [39,[111][112][113]. The most stringent constraint is set by a recent ATLAS analysis finding Br(Z → γγγ) < 2.2 · 10 −6 at 95% CL [39]. Assuming Br(a → γγ) = 1 or 0.1, this constraint sets bounds on the Wilson coefficient |C eff γZ |, which are depicted by the red solid and dashed lines in the left panel of figure 22. The photons have to pass an isolation cut of 4 GeV in transverse energy. However, to be conservative we take the lower bound at 10 GeV as in the h → γγγγ search presented in the same paper. The constraint Br(Z → γγ) < 1.46 · 10 −5 obtained at 95% CL by CDF [113] becomes relevant below m a < 73 MeV, where the two photons are too collimated to be distinguished in the detector. It implies the exclusion regions shown in violet, which has been derived assuming |C eff γγ |/Λ = 1 TeV −1 . ALP decays into lepton pairs give rise to the final states Z → γa → γ + − . OPAL sets the most stringent constraints on these processes, namely Br(Z → γe + e − ) < 5.2·10 −4 , Br(Z → γµ + µ − ) < 5.6·10 −4 and Br(Z → γτ + τ − ) < 7.3·10 −4 at 95% CL [114]. The limits on |C eff γZ | derived from these searches are shown in the right panel of figure 22, assuming Br(a → + − ) = 1.

Probing the ALP-photon and ALP-lepton couplings
Future LHC searches for Z → γa → γγγ decays can probe a large region in the m a − |C eff γγ | parameter space. The green contours in the left panel in figure 23 depict the region where at least 100 signal events are expected at the LHC with √ s = 13 TeV and 300 fb −1 of integrated luminosity. The Z-boson production cross section is σ(pp → Z) = 58.9 nb [115]. The solid, dashed and dotted blue contours correspond to |C eff γZ |/Λ = 1 TeV −1 , 0.1 TeV −1 and 0.01 TeV −1 , respectively. As before, the triangular shape is explained by the fact that ALPs with small masses and couplings are more likely to escape detection. We use Br(a → γγ) = 1 in the plot, but lowering this branching ratio does not change the contours significantly until a critical value is reached, where less than 100 events are produced for all masses and values of C eff γZ . These limiting values are Br(a → γγ) > 7 · 10 −6 (solid), 7 · 10 −4 (dashed) and 0.07 (dotted). To reach couplings smaller than |C eff γZ |/Λ = 0.0026 TeV −1 would require more luminosity. The parameter space shaded in dark green is excluded by present data from CDF [113] and ATLAS [39] (see the left panel of figure 22) under the assumption that |C eff γZ |/Λ = 1 TeV −1 as well as Br(a → γγ) > 0.065 (low-mass region) and Br(a → γγ) > 0.015 (high-mass region).
Comparing the left panel in figure 23 with the corresponding plots in figures 16 and 17 seems to indicate that ALP searches in on-shell Z → γa decays offer the highest sensitivity to the ALP-photon coupling. This is not necessarily true. The point is that, unlike the case of the Higgs-boson decays considered earlier, in the present case the ALP production process Z → γa and the ALP decay process a → γγ are governed by Wilson coefficients C γγ and C γZ , which are correlated via the relations in (2.8), since both couplings originate from the gauge-invariant operators with Wilson coefficients C W W and C BB in (2.1). It is thus very unlikely that |C eff γγ | can take a value that is much smaller than |C eff γZ |. In particular, we note that integrating out a single, complete electroweak multiplet will always generate JHEP12(2017)044 Figure 23. Constraints on the ALP mass and coupling to photons derived from various experiments (colored areas without boundaries, adapted from [24]) along with the parameter regions that can be probed in LHC Run-2 with 300 fb −1 of integrated luminosity using the decay Z → γa → γγγ. We require at least 100 signal events. Left: regions that can be probed are shaded in light green. The contours correspond to |C eff γZ |/Λ = 1 TeV −1 (solid), 0.1 TeV −1 (dashed) and 0.01 TeV −1 (dotted). The dark green regions are excluded by existing measurements assuming that |C eff γZ |/Λ = 1 TeV −1 . The red band shows the preferred parameter space where the (g − 2) µ anomaly can be explained at 95% CL. Right: regions that can be probed in scenarios where the ALP couples only to hypercharge gauge fields (solid blue) or only to SU(2) L gauge fields (solid orange). This plot refers to |C eff γZ |/Λ = 1 TeV −1 .
contributions to C W W and C BB with the same sign. If this is the case, then |C γZ | ≤ c 2 w |C γγ | , (single electroweak multiplet) (6.12) and to very good approximation the same inequality holds for the effective Wilson coefficients including loop corrections. Since |C eff γZ |/Λ > 0.0026 TeV −1 is required to obtain at least 100 signal events, in the presence of the bound (6.12) one cannot probe smaller values of |C eff γγ |. To illustrate this point, we show in the right panel of figure 23 the sensitivity regions obtained for the two cases where the ALP coupling to photons originates only from a coupling to hypercharge (blue line) or only from a coupling to SU(2) L gauge bosons (orange line). In the first case C γZ = −s 2 w C γγ , while in the second one C γZ = c 2 w C γγ . In both cases we have assumed Br(a → γγ) = 1, but the contours are essentially independent of this branching ratio as long as Br(a → γγ) > 1.3 · 10 −4 for U(1) Y and Br(a → γγ) > 1.2 · 10 −5 for SU(2) L . The sensitivity regions are now significantly reduced, but they still cover the parameter space relevant for an explanation of (g − 2) µ .
In the leptonic decay channels, future LHC analyses can search for Z → γa → γ + − decays with = e, µ, τ . Figure 24 shows the regions where at least 100 events are expected in the electron (green), muon (orange) and tau (red) channels (red) for |C eff γZ |/Λ = 1 TeV −1 (solid), 0.1 TeV −1 (dashed) and 0.01 TeV −1 (dotted). We have used Br(a → + − ) = 1 in each case, but as previously the contours are essentially independent of the a → e + e − JHEP12(2017)044 Figure 24. Constraints on the ALP mass and coupling to leptons derived from various experiments (colored areas without boundaries, adapted from [80,81]) along with the parameter region that can be probed using the decay Z → γa → γe + e − . The areas shaded in light green show the reach of LHC Run-2 with 300 fb −1 of integrated luminosity. We require at least 100 signal events. The contours correspond to |C eff γZ |/Λ = 1 TeV −1 (solid), 0.1 TeV −1 (dashed) and 0.01 TeV −1 (dotted). The orange and red regions overlaid in the plots show the corresponding parameter space that can be covered in searches for the decay modes a → µ + µ − and a → τ + τ − (see text for more explanations).

Electroweak precision tests
Since we consider ALPs whose mass is significantly lighter than the electroweak scale, loop corrections to electroweak precision observables can in general not simply be described in terms of the usual oblique parameters S, T and U . Instead, one needs to evaluate the relevant electroweak observables at one-loop order explicitly. Following Peskin and Takeuchi [116], we thus consider the ALP-induced one-loop corrections to three different definitions of the sine squared of the weak mixing angle s 2 w , namely s 2 * defined in terms of the neutral-current couplings ∼ (T f 3 − Q f s 2 * ) of the Z boson to fermions on the Z pole, s 2 W = 1 − m 2 W /m 2 Z defined in terms of the W -and Z-boson masses, and s 2 0 defined via We also consider the ρ * parameter defined by the low-energy ratio of charged-to neutral-current amplitudes. In terms of vacuum-polarization functions defined by the decomposition Π µν AB (q) = Π AB (q 2 ) g µν + O(q µ q ν ), and working to one-loop order, these quantities can be expressed as In the correction terms the lowest-order expressions s 2 w = g 2 /(g 2 + g 2 ) and c 2 w = g 2 /(g 2 + g 2 ) can be used. Note that our relation for s 2 0 differs from a corresponding relation in [117], where the polarization function Π γγ (m 2 Z ) in the first term has been expanded about q 2 = 0. In a new-physics model containing light new particles, such as ours, such an expansion is not legitimate. We find that, at dimension-6 order, the ALP-induced contributions to the vacuum-polarization functions derived from the effective Lagrangian (2.1) involve intermediate (aV ) states with V = γ, Z, W , see the first graph in figure 25. These contributions vanish at q 2 = 0, and hence they do not give a contribution to the ρ * parameter. The individual Π AB (q 2 ) functions are quadratically divergent, however these divergences cancel if we consider the differences between the various definitions of s 2 w . In the class of new-physics models in which the non-polynomial operator (5.10) is present, there is an additional contribution to Π ZZ (q 2 ) shown in the second graph in figure 25, which does not vanish at q 2 = 0, and hence a contribution to the ρ * parameter arises in these models. Setting the ALP mass to zero for simplicity, we obtain 14) and  Figure 26. Allowed regions in the parameters space of the Wilson coefficients C W W − C BB (left) and C γγ − C γZ (right) obtained from a global two-parameter electroweak fit [118] with C (5) Zh = 0 at 68% CL (red), 95% CL (orange) and 99% CL (yellow). We assume that contributions from dimension-6 operators not containing the ALP field can be neglected at Λ = 1 TeV.
where δ 2 = −3, and we have defined The imaginary parts in the above expressions arise from loop graphs containing a photon and an ALP and reflect the existence of the on-shell decay Z → γa considered in section 6.1.
In cross sections these imaginary parts only enter at two-loop order and thus can be omitted here. We can then match the above results with the S, T , U parameters defined in terms of ρ * and the quantities given in (6.14) and (6.15) [116]. This leads to where we have set µ = Λ. The coupling α in the T parameter should be evaluated at q 2 = 0. The presence of UV divergences in these expressions signals that additional short-distance contributions from dimension-6 operators not containing the pseudoscalar a are required in order to cancel the scale dependence. Like in section 4, we will assume that these are small at the new physics scale, since they are not enhanced by the large logarithm ln(Λ 2 /m 2 Z ). Figure 26 shows the allowed parameter space for C (5) Zh = 0 in the plane of the Wilson coefficients C γγ − C γZ (left) and C W W − C BB (right) obtained from the global electroweak fit [118]. The various coefficients are related by (2.8). We observe that the coefficients JHEP12(2017)044 Figure 27. Allowed regions in the parameters space of the Wilson coefficients C W W −C BB obtained from a global three-parameter electroweak fit [118] at 68% CL (red), 95% CL (orange) and 99% CL (yellow). The plots show projections onto the planes where C Zh /Λ = 0 (left), 0.36 TeV −1 (center) and 0.72 TeV −1 (right). We assume that contributions from dimension-6 operators not containing the ALP field can be neglected at Λ = 1 TeV.
C γγ and C BB are largely unconstrained, while C γZ and C W W are restricted to relatively narrow ranges. At 99% CL, we obtain to good approximation |C γZ |/Λ < 6 TeV −1 and |C W W |/Λ < 8 TeV −1 . The flat directions arise because for C W W = 0 (corresponding to C γZ = −s 2 w C γγ ) the contributions to S and U in (6.17) become independent of C BB and C γγ ). We have also performed a global fit for three degrees of freedom including the effect of C (5) Zh . Its contribution to the T -parameter is negative and thus creates a slight tension with the current best fit. To lie within one or two standard deviations of the current best fit point requires |C (5) Zh |/Λ < 0.53 TeV −1 and |C (5) Zh |/Λ < 1.39 TeV −1 respectively. Given the model-independent bound (5.11), the tension is therefore very minor. Figure 27 depicts the results of this fit projected onto the planes where C  Another precision test can be performed by considering the running of the electromagnetic coupling constant from q 2 = 0 to q 2 = m 2 Z . In our model we obtain where the vacuum-polarization functions now contain the ALP contribution only. Dropping again a small imaginary part and setting µ = Λ, we find A measurement of α(m Z ) has been performed by the OPAL collaboration at a center of mass energy of 193 GeV [119]. The precision of this measurement is at the percent level, which is still compatible with values of C W W and C BB of O(30/TeV). A significant improvement on the precision is expected from a future circular e + e − collider FCC-ee [120], which will be able to measure Z-pole observables with unprece-JHEP12(2017)044 Figure 28. Allowed regions in the parameters space of the Wilson coefficients C W W − C BB (left) and C γγ − C γZ (right) obtained from projections for the two-parameter global electroweak fit at a future FCC-ee machine [121] at 68% CL (red), 95% CL (orange) and 99% CL (yellow), setting C Zh = 0. We assume that contributions from dimension-6 operators not containing the ALP field can be neglected at Λ = 1 TeV. For the parameter space within the dashed black contour, a FCC-ee measurement of α(m Z ) is within its projected errors at 95% CL [120]. dented precision. In particular, α(m Z ) can be determined with an uncertainty of about 10 −5 . In figure 28, we show projections for the two-parameter electroweak fit based on the data obtained at such a machine [121], assuming that the central values of C W W and C BB vanish. In the same figure, we superimpose the expected 95% CL bound derived from the measurement of α(m Z ) (dashed contours), assuming that the theoretical error on this quantity will have decreased below the experimental uncertainty by the time the measurement can be performed. Combining these measurements can constrain |C W W |/Λ < 2 TeV −1 and |C BB |/Λ < 3 TeV −1 , or equivalently |C γγ |/Λ < 2.5 TeV −1 and |C γZ |/Λ < 1 TeV −1 (at 95% CL).

Conclusions
Pseudoscalar particles with an approximate shift symmetry, so-called axion-like particles (ALPs), appear as pseudo Nambu-Goldstone bosons in any theory in which a global symmetry is spontaneously broken. If the mass scale of new physics is high, a light pseudo Nambu-Goldstone boson could be a harbinger of a new UV sector, which cannot otherwise be probed directly. The discovery of an ALP would not only confirm the existence of a UV theory beyond the SM, but by measuring its couplings important information on the properties of this theory can be derived.
Based on the most general effective Lagrangian for a pseudoscalar with an approximate shift symmetry (softly broken only by an explicit mass term), we have computed the partial decay widths for ALPs into pairs of photons, leptons, jets and heavy quarks at one-loop order. Since the decay a → γπ is not allowed, relevant hadronic decay modes only open JHEP12(2017)044 up when the ALP is heavy enough to decay into three pions. We have calculated the a → πππ partial widths of the ALP in terms of a chiral Lagrangian for the first time. We have emphasized that even loop-suppressed Wilson coefficients can lead to non-negligible branching ratios of ALPs decaying into photons or leptons. The assumption of stable ALPs therefore becomes unrealistic above a certain mass, if sizable couplings C ii /Λ of order (0.01 − 1) TeV −1 to any SM fields exist. For the same reason, an ALP with such couplings cannot be lighter than about 1 MeV, since in the presence of loop corrections it is impossible to satisfy the very strong cosmological bounds on the ALP-photon coupling without excessive fine tuning.
Significant insights can be gained by considering the exotic, on-shell decays h → Za, h → aa and Z → γa. These three decays offer complementary information on a possible UV sector beyond the SM. While Z decays are induced by dimension-5 operators coupling the ALP to electroweak gauge bosons, Higgs decays probe the dimension-6 Higgs portal in the case of h → aa and dimension-5 or 7 operators in the case of h → Za. The non-polynomial dimension-5 operator in (5.10) inducing h → Za decay at Born level only arises if the heavy particles in the UV theory obtain a dominant fraction of their mass from electroweak symmetry breaking. Discovering an ALP in any or a combination of these exotic decays therefore allows us to extract non-trivial details about the underlying UV theory. To see the important role of Higgs decays, consider as a concrete example a scenario in which the only tree-level ALP couplings to SM fields are flavor-universal couplings to the up-type quarks, c uu = c cc = c tt = Λ/TeV. ALP couplings to other SM particles are induced only by means of quark loops. Assuming m a = 1 GeV, one then finds the Higgs branching ratios Br(h → Za) = 2.5 × 10 −4 and Br(h → aa) = 8.5 × 10 −3 , which are 0.15 and 5.5 times the SM h → γZ branching ratio, respectively. The loop-induced couplings to electroweak gauge bosons are rather small, and correspondingly Br(Z → γa) = 4.8 × 10 −9 is absolutely negligible. On the other hand, the loop-induced ALP-photon coupling C eff γγ ≈ 0.008 lies in the range of sensitivity of our approach. In this scenario it would be easy to discover the ALP in h → aa decay, challenging to probe its couplings in h → Za decay, and hopeless to see any hints of ALPs in Z → γa decay.
We have presented a comprehensive discussion of the LHC reach in searches for ALPs in exotic Higgs-and Z-boson decays. Taking into account constraints from existing searches, model-independent bounds on non-SM Higgs or Z decays and finite-lifetime effects of light ALPs or ALPs with small couplings, we found LHC searches for the decays h → Za, h → aa and Z → γa to be sensitive to new-physics scales as high as 100 TeV for ALP masses in the GeV range. Depending on the decay mode of the ALP, several striking signatures can be observed. We have especially considered subsequent ALP decays into photons and charged leptons, taking into account the possibility that light boosted ALPs decay into collimated photon jets, which cannot be distinguished from a single photon experimentally. Cosmological bounds, ALP searches with helioscopes, beam-dump experiments and searches for ALPs at lepton and hadron colliders significantly constrain the parameter space for ALPs decaying into di-photons or e + e − pairs. Intriguingly, we project the best sensitivity for ALP searches in on-shell Higgs-and Z-boson decays at the LHC for ALP masses in the range above approximately 10 MeV and up to about 90 GeV, a region of parameter space JHEP12(2017)044 mostly unconstrained by existing bounds once we assume that the relevant ALP couplings are of order 1/TeV or smaller. For ALPs in the GeV mass range, this reach extends many orders of magnitude beyond current bounds, without the need to assume any large Wilson coefficients. Even with loop-suppressed ALP-Higgs couplings, the bounds on the ALP-photon coupling can be improved by up to five orders of magnitude using searches for the decays h → Za → + − + γγ and h → aa → 4γ. Improvements by several orders of magnitude can also be obtained from a search for the decay Z → γa → 3γ; however, the reach in this case depends on the correlation of the aγγ and aγZ couplings, which depends on the underlying UV model. Importantly, these bounds can be derived even if the a → γγ branching ratio is significantly less than 1. In the leptonic decay channels a → + − , completely uncharted territory in parameter space can be probed, extending down to ALP-lepton couplings as small as (10 6 TeV) −1 .
We have further computed the parameter space for which the long-standing (g − 2) µ anomaly can be explained by ALPs coupling to muons and photons. A possible resolution by a loop contribution from ALPs is largely independent of its mass and requires a sizable coupling to photons and a coupling of similar size (and the correct sign) to muons. For example, a good fit can be found for m a = 1 GeV and C γγ ≈ −c µµ ≈ 1.5 (Λ/TeV). We have translated the bound from a Babar search for a new Z boson in the e + e − → µ + µ − + µ + µ − channel into a constraint on the c µµ − C γγ plane, thereby directly constraining a possible explanation of (g − 2) µ by ALP exchange. We find that future searches for e + e − → µ + µ − + µ + µ − as well as e + e − → µ + µ − + γγ at Belle II have the potential to discover or exclude an ALP explanation of the anomaly for 2m µ < m a 2 GeV. Remarkably, the complete unconstrained parameter space for which an ALP can explain the muon anomaly can be probed by the exotic Higgs-and Z-boson decays studied in this paper. Barring for scenarios in which the aγγ coupling is very large, whereas the aZh, aah and aγZ couplings are all more than one-loop suppressed, searches for ALPs in exotic decays of on-shell Higgs and Z bosons at the LHC can therefore exclude or confirm an ALP explanation of (g − 2) µ . Electroweak precision tests constrain the ALP couplings to electroweak gauge bosons and to Zh. These coefficients control the Z → γa and h → Za decay rates, respectively. We have computed the one-loop corrections to the oblique parameters and to α(m Z ) and derived the corresponding bounds on the Wilson coefficients from the global electroweak fit, finding that they are rather weak. We have also presented projections for a future FCC-ee machine, where it will be possible to probe ALP couplings to electroweak gauge bosons of order 1/TeV. The LHC has an unprecedented reach in searching for ALPs in exotic, on-shell decays of Higgs and Z bosons. We strongly encourage experimental searches in the full mass range and in all three channels discussed in this paper. A UFO file for the ALP model discussed in the present work is available from the authors upon request.
Note added. We would like to thank the referee for valuable comments on the Edelweiss and BaBar bounds and for encouraging us to include figures 18 and 20. After the submission of this paper, two new analyses discussing first experimental results in heavyion collisions [135] and bounds from beam-dump searches [136] have been submitted to JHEP12(2017)044 the arXiv. These results supersede some of the constraints on the ALP-photon coupling shown in our work, but they do not change any of the conclusions derived in this paper. Physics, Fundamental Interactions and Structure of Matter (PRISMA -EXC 1098), and grant 05H12UME of the German Federal Ministry for Education and Research (BMBF).

A Naive dimensional analysis estimates
Here we collect order-of-magnitude estimates for the Wilson coefficients in the effective Lagrangians (2.1) and (2.6) based on naive dimensional analysis. Using the counting rules derived in [122][123][124], one obtains Zh = (4π) 3C (7) Zh , (A.1) where the subscript V is the second relation can be G, W or B. The barred coefficients on the right-hand sides of these relations can naturally be of O(1) in strongly coupled theories. When the effective Lagrangians are rewritten in terms of a parameter f defined such that 4πf ≡ Λ (this parameter is related to the ALP decay constant f a by f = −2C GG f a ), one obtains expressions analogous to (2.1) and (2.6), in which the Wilson coefficients are replaced by the barred Wilson coefficients and Λ is replaced by f . The only exception are the ALP-gauge-boson couplings, which are given byC V V /(4π) 2 . It would therefore have been more natural to introduce a loop factor 1/(4π) 2 in the three terms shown in the second line of (2.1). 3 Following a standard practice in the ALP literature, we have refrained from doing so. In light of these remarks, it becomes evident that an explanation of the (g−2) µ anomaly requires a somewhat unnaturally large value of the ALP-photon coupling. From figure 6 we see that we typically need |C γγ |/Λ 0.5/TeV, corresponding to |C γγ |/Λ 6/TeV. Generating such a large coefficient may require to have a large multiplicity of new TeVscale particles in a loop or lowering Λ below the TeV scale, but it does not appear to be impossible.
Extensions of the SM in which the electroweak symmetry is realized non-linearly provide an explicit example of strongly coupled models, in which operators of higher dimension in the effective Lagrangian are suppressed by powers of 1/f rather than 1/Λ [125]. In realistic composite Higgs scenarios the ratio ξ = v 2 /f 2 is tightly constrained by electroweak JHEP12(2017)044 precisions tests, implying ξ < 0.05 [126], and Higgs phenomenology, yielding ξ < 0.1 [101], both at 95% CL. As a result, it is unlikely that f can be significantly below the TeV scale in these models [127].

B Couplings of a light ALP to hadrons
At energies below a few GeV, the effective Lagrangian (2.1) supplemented by the QCD Lagrangian gives rise to the terms where m a,0 denotes a possible ALP mass term resulting from an explicit breaking of the shift symmetry. We will for simplicity only consider the two light u and d quarks. We use a compact matrix notation, where in the mass basis m q = diag(m u , m d ) and c qq = diag(c uu , c dd ) are diagonal hermitian matrices. Before mapping this expression onto an effective chiral Lagrangian, it is convenient to remove the ALP-gluon coupling by means of the chiral rotation q → exp iκ q a 2f a γ 5 q , where κ q is a diagonal matrix satisfying tr κ q = 1, and f a is referred to as the ALP decay constant. Under the chiral rotation the measure of the path integral is not invariant [128,129], and this generates extra terms adding to the anomalous couplings in (B.1). In order to remove the ALP-gluon coupling we need to require that wherê m q (a) = exp iκ q a 2f a γ 5 m q exp iκ q a 2f a γ 5 ,ĉ qq = c qq + 32π 2 κ q C GG . (B.5) Matching the above effective Lagrangian onto a chiral Lagrangian, one obtains [52,60,130] where Σ containing the pion fields has been defined after (3.6). Note that nowm q (a) is evaluated by replacing γ 5 in (B.5) by its eigenvalue +1. The covariant derivative is defined as D µ Σ = ∂ µ Σ − ieA µ [Q, Σ], where Q contains the quark electric charges in units of e. Even if the explicit mass term m a,0 is absent, QCD dynamics generates a mass for the ALP [6,55,56], thereby breaking the continuous shift symmetry. Expanding the terms in the first line to quadratic order in the pion and ALP fields, one finds the mass eigenvalues where we have adopted the choice (3.7) for the κ q parameters, which eliminates the mass mixing of the ALP with the neutral pion. This choice leads to the effective chiral Lagrangian given in (3.6). The coefficient in front of the ALP-photon coupling now takes the form where E/N = C γγ /C GG and we have used that m u /m d ≈ 0.48 (see e.g. [131,132] for two recent lattice determination of this ratio). The term proportional to the explicit isospin breaking caused by the mass difference between up and down quarks results from the coupling of the neutral pion to G A µνG µν,A . The corresponding matrix element has been evaluated in [133] and is found to be The pion then decays into two photons via the axial anomaly. The contribution 5/3 arises from an analogous coupling the flavor-singlet meson ϕ 0 (the analogue of η 1 in flavor SU(3)) [134]. Next-to-leading order corrections to the result (B.8) have been worked out in [59]. They lead to a coefficient [E/N − (1.92 ± 0.04)], which we use in our numerical analysis.

C Technical details of the loop calculations
The loop function g(τ ) entering the expression for the effective ALP-lepton coupling in (3.13) is given by the parameter integral g(τ ) = 5 + 4 3 where τ = 4m 2 /m 2 a − i0. The asymptotic expansions for small and large values of τ have been shown in (3.14).

D Effect of a finite ALP lifetime
The two event fractions defined in (5.3) obey the exact relations where the function F (x) is given by Using the first result, we have obtained the asymptotic relations given in (5.4).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.