Leptophilic ALPs in Laboratory Experiments

Abstract

We study the collider phenomenology of leptophilic axion-like particles (ALPs), i.e. pseudoscalar particles that couple only to charged leptons. Loops of charged leptons induce effective interactions of the ALPs with photons, which depend on the momenta of the interacting particles and differ between pseudoscalar and derivative lepton couplings. We systematically discuss the form of the interaction with photons for general external momenta and identify the regimes when it can be safely approximated by an effective coupling constant. We use these results to derive novel constraints from LEP and calculate state-of-the-art limits from E137 and NA64 for four different scenarios, in which the ALPs couple either to a single lepton generation or universally to all, for both pseudoscalar and derivative lepton couplings. We collect complementary bounds from astrophysics, flavour, and other laboratory experiments to chart the allowed parameter space of leptophilic ALPs in the MeV-GeV mass range.

1 Introduction

Some of the simplest modifications of the Standard Model (SM) involve a single axion-like particle (ALP), which is a pseudo-Goldstone boson of a suitable, spontaneously broken global symmetry. In contrast to the QCD axion Wilczek:1977pj ; Weinberg:1977ma ; Peccei:1977hh ; Peccei:1977ur , this symmetry is subject to a large explicit breaking, which induces an ALP mass that is simply taken as a free parameter. Of particular interest are masses in the MeV–GeV range, which lead to a very rich phenomenology at particle colliders Bauer:2017ris , depending on the precise form of ALP interactions with SM particles. Though the explanation of the dark matter abundance is one of the main motivations of these scenarios (see, e.g. refs. Nomura:2008ru ; Freytsis:2010ne ; Dolan:2014ska ; Berlin:2015wwa ; Fan:2015sza ; No:2015xqa ; Bauer:2017ota ; Baek:2017vzd ; Kamada:2017tsq ; Kaneta:2017wfh ; Banerjee:2017wxi ; Arcadi:2017wqi ; Hochberg:2018rjs ; Berlin:2018bsc ; deNiverville:2019xsx ; Buttazzo:2020vfs ; Darme:2020sjf ; Ge:2021cjz ; Gola:2021abm ; Domcke:2021yuz ; Zhevlakov:2022vio ; Bauer:2022rwf ; Panci:2022wlc ; Bharucha:2022lty ; Fitzpatrick:2023xks ; Ghosh:2023tyz ; Dror:2023fyd ; Capozzi:2023ffu ; Armando:2023zwz ), it is important to map current limits from laboratory experiments, astrophysics and cosmology into the general parameter space of the ALP Lagrangian.

While this task has been performed in much detail for ALP couplings to photons, couplings to leptons have received far less attention (notable exceptions include refs. Altmannshofer:2022ckw ; Bauer:2021mvw ; Ferreira:2022xlw ; Armando:2023zwz ; Cornella:2019uxs ; Calibbi:2020jvd ; Bertuzzo:2022fcm ). In the present work, we aim to systematically study the model-independent laboratory constraints on “leptophilic ALPs”, which, by definition, couple dominantly to charged leptons. As a result, we provide an exhaustive overview of the relevant constraints from beam-dump, flavour and collider experiments on leptophilic ALPs, which complement limits from astrophysics Raffelt:2006cw ; Carenza:2021pcm ; Ferreira:2022xlw ; Fiorillo:2025sln .

We restrict ourselves to flavour-diagonal couplings and consider four simple scenarios, where either the coupling to a single lepton family dominates ($e,\mu,\tau$) or the ALP couples in a flavour-universal way to leptons. In contrast to previous analyses, we explicitly consider two cases of possible leptonic ALP interactions, which take the form of either pseudoscalar or derivative couplings. It is well-known that these interactions are largely equivalent for tree-level processes involving a single ALP but differ for multi-ALP interactions and, in particular, by an effective ALP-photon coupling induced by the chiral anomaly. Therefore, the two leptophilic interactions lead to different predictions for observables that involve an effective ALP-photon-photon interaction, which is induced by a leptonic loop. The associated vertex function depends on the external momenta and, thus, on the typical energy scales of the relevant process (apart from a constant part depending on whether pseudoscalar or derivative couplings are considered). This energy dependence implies that the phenomenology of the leptophilic ALP cannot be simply obtained from existing limits on the ALP-photon coupling but has to be re-calculated using the full momentum dependence for all processes involving ALP-photon interactions Bharucha:2022lty ; Ferreira:2022xlw . The complete analysis of the resulting limits from E137 Bjorken:1988as , NA64 Gninenko:2719646 ; NA64:2019auh and LEP Jaeckel:2015jla represents the main result of this work.

Although such loop-induced processes are suppressed with respect to tree-level diagrams involving only leptons, we will show that beam dump experiments like E137 can give highly relevant constraints on leptophilic ALP, even when both the production and decay of the ALPs proceed exclusively via photon couplings. When the relevant energy scales are much larger than the ALP and lepton masses, the loop function becomes suppressed. The resulting vertex inherits this suppression for pseudoscalar interactions, while the derivative vertex is essentially given by the constant part induced by the chiral anomaly. This implies that, e.g. limits from high-energy colliders like LEP that have been derived for effective photon couplings Jaeckel:2015jla are relevant only for leptophilic ALPs with derivative couplings, and similar conclusions hold for highly energetic leptonic beam dumps. More generally, a large separation of scales often permits the approximation of the fully off-shell vertex function by the on-shell coupling, which largely simplifies calculations, as we discuss in detail. We also point out the great potential of muon beam-dump experiments, such as NA64$\mu$ NA64:2024klw , to constrain models of leptophilic ALPs Li:2025yzb . Indeed, NA64$\mu$ already sets relevant constraints on several of the scenarios that we consider, which will substantially improve in coming years with larger data sets.

The rest of this work is organised as follows. In Section 2, we define the basic setup of the leptophilic ALP models and discuss the difference between the pseudoscalar and derivative basis and the resulting ALP-photon-photon vertex function. In Section 3, we use these results in order to calculate ALP production at leptonic beam dumps in all relevant channels, focusing on E137 and NA64. Constraints on leptophilic ALPs from other experiments are discussed in Section 4, and in particular, we derive novel limits by recasting LEP1 data on multi-photon events (largely following ref. Jaeckel:2015jla ). Here, we also collect existing bounds in the literature from $B$ factories, $W^{+}$ boson and rare meson decays, SN 1987A, leptonic anomalous magnetic moments and cosmology. Our results are presented in Section 5, which mainly consist of novel limits on leptophilic ALPs from E137 and NA64 that have been carefully derived including the track-length distribution of the beam electrons as well as secondary electrons and positrons, which we combine with other existing limits in order to provide a concise overview of all relevant constraints on leptophilic ALPs. We summarise our conclusions in Section 6 and include various appendices with details on the experimental setup in appendix A and the cross-section calculations in appendix B. We also compare our results for the exclusion limit from E137 with previous analyses in appendix C.

2 Setup

In this section, we introduce the effective Lagrangians for leptophilic ALP models, discuss in detail the phenomenological and theoretical differences between pseudoscalar and derivative couplings and study the general induced one-loop ALP couplings to photons.

2.1 ALP Couplings to Fermions

The ALP is assumed to be a pseudoscalar that is not charged under the SM gauge group, arising as a pseudo-Goldstone boson from the spontaneous breaking of a global Peccei-Quinn (PQ) symmetry at some UV scale $\Lambda$. At scales much smaller than $\Lambda$, the effective Lagrangian is governed by the non-linearly realised PQ symmetry under which the ALP shifts, $a\to a+{\rm const}$. While for the QCD axion, this symmetry is only broken by QCD instantons, for the ALP that we consider here, there is also a soft breaking by an explicit ALP mass term. Apart from this mass term, the most general Lagrangian then involves shift-invariant couplings to fermions and anomalous couplings to gauge bosons, just as the QCD axion Georgi:1986df .

Here, we focus on leptophilic ALPs, which are described by an effective Lagrangian only involving charged leptons. Up to higher-dimensional operators and restricting to flavour-diagonal couplings, the general ALP Lagrangian reads

$$\mathcal{L}_{\rm D} =\frac{1}{2}\partial_{\mu}a\,\partial^{\mu}a-\frac{1}{2}m_{a}^{2}a^{2}+\frac{\partial_{\mu}a}{2\Lambda}\sum_{\ell=e,\mu,\tau}C_{\ell}\,\overline{\psi}_{\ell}\gamma^{\mu}\gamma_{5}\psi_{\ell}\,,$$ (1)

with real couplings $C_{\ell}$. Equivalently, one can change the field basis and write this Lagrangian in terms of pseudoscalar couplings to fermions. Under the chiral ALP-dependent rotation of the lepton fields

$$\psi_{\ell} \to e^{i\frac{C_{\ell}}{2\Lambda}a\gamma_{5}}\psi_{\ell}\,, \overline{\psi}_{\ell} \to\overline{\psi}_{\ell}\,e^{i\frac{C_{\ell}}{2\Lambda}a\gamma_{5}}\,,$$ (2)

the derivative couplings above are cancelled when plugging the transformed fields into the lepton kinetic terms, and the ALP appears in the lepton mass term. Since this transformation is anomalous, i.e. it does not leave the fermion path integral measure invariant, it also gives rise to a contribution to the ALP couplings to photons, so that in the new field basis, the ALP Lagrangian reads (see e.g. refs. Bauer:2017ris ; Armando:2023zwz for details)

$$\mathcal{L}_{\rm D}=\frac{1}{2}\partial_{\mu}a\,\partial^{\mu}a-\frac{1}{2}m_{a}^{2}a^{2}-ia\sum_{\ell=e,\mu,\tau}\frac{m_{\ell}}{\Lambda}C_{\ell}\overline{\psi}_{\ell}\gamma_{5}\psi_{\ell}+\frac{\alpha}{4\pi}\frac{a}{\Lambda}F^{\mu\nu}\tilde{F}_{\mu\nu}\sum_{\ell=e,\mu,\tau}C_{\ell}\,,$$ (3)

up to higher-order terms in $a/\Lambda$, and $\tilde{F}^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}F_{\rho\sigma}/2$ with $\epsilon^{0123}=1$. Note that the truncation explicitly breaks the shift symmetry of lepton couplings, which is restored by higher-order terms where the shift symmetry just corresponds to a re-phasing of the leptons.

It is well-known Chisholm:1961tha ; Kamefuchi:1961sb that all physical observables remain unchanged under non-linear field redefinitions of the form in eq. (2), which is just a change of variables in the path integral formulation of QFT. In the context of effective field theories, it has been shown in ref. Arzt:1993gz that such field redefinitions are largely equivalent to using the equations of motion for the involved fields (even when off-shell), although care has to be taken beyond the leading order Criado:2018sdb . Working consistently at leading order in $1/\Lambda$, i.e. only using vertices involving a single ALP, the Lagrangians in eq. (1) and eq. (3) therefore make identical predictions, even though the latter features ALP couplings to both leptons and photons.

This equivalence implies that a leptophilic Lagrangian in the pseudoscalar basis defined by

$$\mathcal{L}_{\rm P}=\frac{1}{2}\partial_{\mu}a\,\partial^{\mu}a-\frac{1}{2}m_{a}^{2}a^{2}-ia\sum_{\ell=e,\mu,\tau}\frac{m_{\ell}}{\Lambda}C_{\ell}\overline{\psi}_{\ell}\gamma_{5}\psi_{\ell}\,,$$ (4)

does not give the same results as (1) when the corresponding process involves effective vertices to photons. For this reason, there are two “leptophilic” ALP models with distinct phenomenology, which in the following we refer to as the “derivative” (eq. (1)) and “pseudoscalar” (eq. (4)) model. To follow the conventional notation in the literature, we also introduce the dimensionless coupling $g_{a\ell}$ via the relation

$$\frac{C_{\ell}}{\Lambda}=\frac{g_{a\ell}}{m_{\ell}}\,,$$ (5)

such that the pseudoscalar interactions can be written as $\mathcal{L}_{\text{P}}\supset-ig_{a\ell}a\overline{\psi}_{\ell}\gamma_{5}\psi_{\ell}$. In order to clearly distinguish between derivative and pseudoscalar leptophilic models, we use the notation $C_{\ell}/\Lambda$ exclusively for the derivative interaction and the coupling $g_{a\ell}/m_{\ell}$ exclusively for the pseudoscalar interaction. Moreover, we will discuss below four different scenarios, depending on the values of $C_{\ell}$ for the three families. Either the ALP couples to a single lepton species, e.g. $C_{e}\neq 0$, while $C_{\mu}=C_{\tau}=0$, or it has flavour-universal couplings in the derivative basis, i.e. $C_{e}=C_{\mu}=C_{\tau}$, corresponding to universal PQ charges.

Finally, we briefly discuss the question of which type of leptophilic interaction (eq. (1) or eq. (4)) is best motivated from a UV perspective. Simple UV completions seem to suggest the pseudoscalar form, which can arise from a renormalisable Yukawa coupling with a non-linearly realised Goldstone boson in a straightforward way. For example, one can imagine a standard QCD axion scenario like the DFSZ model Zhitnitsky:1980tq ; Dine:1981rt , in addition to some small explicitly PQ breaking mass term, for example induced by $M_{\rm Planck}$-suppressed operators that generally spoil the solution to the Strong CP Problem in QCD axion models (the so-called PQ Quality Problem Kamionkowski:1992mf ). In DFSZ scenarios, SM fermions are charged under the PQ symmetry, which is spontaneously broken by the Higgs fields (at least two) and a SM singlet scalar. The axion is then a linear combination of the pseudoscalar components in the Higgs fields and the singlet and couples to fermions through the Yukawa interactions (see e.g. refs. DiLuzio:2020wdo ; Badziak:2021apn for the explicit construction of generalised DFSZ models)

$$\mathcal{L}_{\textrm{DFSZ}}\supset-\bar{u}_{\textrm{L}}m_{U}u_{\textrm{R}}e^{i\chi_{H_{u}}\frac{a}{v_{\rm PQ}}}-\bar{d}_{\textrm{L}}m_{D}d_{\textrm{R}}e^{i\chi_{H_{d}}\frac{a}{v_{\rm PQ}}}-\bar{e}_{\textrm{L}}m_{E}e_{\textrm{R}}e^{i\chi_{H_{d}}\frac{a}{v_{\rm PQ}}}+\text{h.c.}\,,$$ (6)

where $\chi_{H_{u,d}}$ are the PQ charges of the Higgs doublets $H_{u,d}$ and $v_{\rm PQ}\propto\Lambda$ is the PQ breaking scale.

Expanding the exponential functions, one obtains mass terms for the fermions as well as pseudoscalar ALP-fermion interaction terms, while a tree-level coupling between ALPs and photons is absent in this basis. A chiral redefinition of the fermion fields would trade the pseudoscalar couplings of eq. (6) for derivative couplings and generally non-vanishing couplings to gauge bosons, unless the PQ charges are such that the electromagnetic (EM) anomaly coefficient vanishes. In the four leptophilic scenarios discussed above, this anomaly is indeed non-vanishing since, by construction, only leptons contribute, and the leptonic PQ charge matrix is not traceless. Thus, in order to obtain one of these scenarios with derivative lepton couplings and a vanishing photon coupling, one would need to extend the simple DFSZ UV model with new heavy fermions that have suitable PQ charges to cancel the EM anomaly, along the lines of common KSVZ models Kim:1979if ; Shifman:1979if .

On the other hand, since the derivative coupling is the only dimension-5 operator that preserves the shift-symmetry of the ALP, other shift-symmetric couplings (such as those to gauge bosons, e.g. in a KSVZ model) will generally induce this form of coupling through RGE-evolution Chala:2020wvs . An example is the photophobic ALP discussed in ref. Craig:2018kne , which, however, certainly requires more model-building effort than the above-mentioned DFSZ-like constructions. Moreover, as discussed in ref. Craig:2018kne , it comes with additional interactions with the electroweak gauge bosons (which we do not consider here).

In the remainder of this article, we choose to remain agnostic about the nature of the UV completions of both bases and simply derive constraints for both scenarios in eq. (1) and eq. (4), which will deviate when probing ALP couplings to photons. These couplings are effectively generated from lepton loops and differ for the two interactions by a constant, as we discuss in more detail now.

2.2 Effective ALP Coupling to Photons

Figure 1: Triangle diagram inducing the effective ALP-photon interaction.

Even though leptophilic ALPs do not couple to photons at tree level, an effective ALP-photon interaction is induced by the triangle diagram shown in figure 1, which corresponds to an off-shell ALP-photon-photon vertex function $\Gamma_{a^{*}\gamma^{*}\gamma^{*}}^{\mu\nu}(q_{1},q_{2})$ that depends on the outgoing photon momenta $q_{1,2}$, with incoming ALP momentum $q_{a}=q_{1}+q_{2}$. It is convenient to introduce an effective coupling as a scalar form factor by

$$\Gamma_{a^{*}\gamma^{*}\gamma^{*}}^{\mu\nu}(q_{1},q_{2})=ig_{a^{*}\gamma^{*}\gamma^{*}}^{\rm eff}(q_{a}^{2},q_{1}^{2},q_{2}^{2})q_{1\alpha}q_{2\beta}\epsilon^{\mu\nu\alpha\beta}\,,$$ (7)

such that $g_{a^{*}\gamma^{*}\gamma^{*}}^{\rm eff}(q_{a}^{2},q_{1}^{2},q_{2}^{2})$ is given by a scalar one-loop integral with three external momenta $q_{a},q_{1},q_{2}$, which can be expressed in terms of Passarino-Veltman functions. For ALPs with a derivative coupling to a single lepton species of the form in eq. (1), the effective off-shell coupling reads Ferreira:2022xlw (the superscript “D” denotes the derivative basis)

$$g_{a^{*}\gamma^{*}\gamma^{*}}^{\rm eff,D}(q_{a}^{2},q_{1}^{2},q_{2}^{2}) =\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\quantity[1+2m_{\ell}^{2}C_{0}\quantity(q_{1}^{2},q_{2}^{2},q_{a}^{2},m_{\ell}^{2},m_{\ell}^{2},m_{\ell}^{2})]\,,$$ (8)

where $C_{0}$ is the Passarino-Veldman three-point function defined as

$$\begin{split}C_{0}&\quantity(q_{1}^{2},q_{2}^{2},(q_{1}+q_{2})^{2},m_{1}^{2},m_{2}^{2},m_{3}^{2})\equiv\int\frac{\differential^{4}k}{i\pi^{2}}\frac{1}{\quantity(k^{2}-m_{1}^{2})\quantity(k_{1}^{2}-m_{2}^{2})\quantity(k_{12}^{2}-m_{3}^{2})}\,,\end{split}$$ (9)

with $k_{1}\equiv k-q_{1}$ and $k_{12}=k-q_{1}-q_{2}$, following the conventions in ref. Shtabovenko:2023idz . If the ALP couples to several leptons, the effective coupling is given by the sum over all contributions.

Since the Lagrangian in eq. (1) is equivalent to the one in eq. (3) for single ALP vertices, the scalar vertex function for the pseudoscalar coupling in eq. (4) differs from eq. (8) by a constant, which is given by the coefficient of the $\tfrac{1}{4}aF\tilde{F}$ operator in eq. (3), $\alpha C_{\ell}/(\pi\Lambda)$. Indeed, for the effective coupling in the pseudoscalar basis (denoted by the superscript “P”), one has

$$g_{a^{*}\gamma^{*}\gamma^{*}}^{\rm eff,P}(q_{a}^{2},q_{1}^{2},q_{2}^{2}) =\frac{\alpha}{\pi}\frac{g_{a\ell}}{m_{\ell}}\quantity[2m_{\ell}^{2}C_{0}\quantity(q_{1}^{2},q_{2}^{2},q_{a}^{2},m_{\ell}^{2},m_{\ell}^{2},m_{\ell}^{2})]\,,$$ (10)

where we used $g_{a\ell}/m_{\ell}=C_{\ell}/\Lambda$.

It is instructive to evaluate the fully off-shell couplings for the simplifying cases where one or more external particles are on-shell. If a single photon is on-shell, we find (in agreement with ref. Bauer:2021mvw )

$$g_{a^{*}\gamma^{*}\gamma}^{\rm eff,D}(q_{a}^{2},q_{1}^{2}) =\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\left(1-\int_{0}^{1}dx\int_{0}^{1}dy\frac{m_{\ell}^{2}}{m_{\ell}^{2}-x(1-x)(q_{a}^{2}-2yq_{a}q_{1})}\right)\,,$$ (11)

which is valid in the derivative basis, i.e. using the Lagrangian in Eq. (1). For the analogous result in the pseudoscalar basis, one simply needs to drop the constant part.

Taking the ALP on-shell, $q_{a}^{2}=m_{a}^{2}$, the vertex function depends only on a single off-shell photon squared momentum $q_{1}^{2}\equiv t$, and reads

$$g^{\rm eff,D}_{a\gamma^{*}\gamma}(t) =\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\left(1-\int_{0}^{1}dx\int_{0}^{1}dy\frac{m_{\ell}^{2}}{m_{\ell}^{2}-x(1-x)(m_{a}^{2}-y(m_{a}^{2}-t))}\right)\,.$$ (12)

Finally, if all particles are on-shell, i.e. if $t=0$, the vertex function becomes a constant

$$g^{\rm eff,D}_{a\gamma\gamma} =\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\left(1-\int_{0}^{1}dx\int_{0}^{1}dy\frac{m_{\ell}^{2}}{m_{\ell}^{2}-x(1-x)(1-y)m_{a}^{2}}\right)=\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\left(1-\tau_{\ell}f(\tau_{\ell})^{2}\right)\,,$$ (13)

with $\tau_{\ell}=4m_{\ell}^{2}/m_{a}^{2}$ and the function

$$f(\tau)=\begin{cases}\arcsin(\frac{1}{\sqrt{\tau}}),&\mbox{\quad if\quad}\tau\geq 1\\ \frac{\pi}{2}+\frac{i}{2}\log(\frac{1+\sqrt{1-\tau}}{1-\sqrt{1-\tau}}),&\mbox{\quad if\quad}0<\tau<1\end{cases}\,,$$ (14)

consistent with the expression derived in refs. Ferreira:2022xlw ; Bauer:2017ris . This constant can be identified with the effective Lagrangian coupling

$$\mathcal{L}_{\gamma} =\frac{g_{a\gamma\gamma}^{\rm eff}}{4}aF^{\mu\nu}\tilde{F}_{\mu\nu}\,.$$ (15)

and controls the ALP decay into two photons. Since, asymptotically,

$$\tau f^{2}(\tau)=\begin{cases}1+\frac{1}{3\tau}&\tau\gg 1\\ \frac{\tau}{4}\left(\pi+i\ln\frac{4}{\tau}\right)^{2}&\tau\ll 1\end{cases}\,,$$ (16)

the effective photon coupling defined in eq. (15) for a derivatively coupled ALP in eq. (1) has the limits

$$g_{a\gamma\gamma}^{\rm eff,D} \propto\ 1-\tau_{\ell}f(\tau_{\ell})^{2}\propto\begin{cases}-\tfrac{1}{12}\,(m_{a}/m_{\ell})^{2}&\mbox{\quad if\quad}m_{a}\ll m_{\ell}\\ 1&\mbox{\quad if\quad}m_{a}\gg m_{\ell}\end{cases}\,.$$ (17)

Instead, the effective photon coupling for an ALP coupled to the pseudoscalar current in eq. (4) has the asymptotic behaviour

$$g_{a\gamma\gamma}^{\rm eff,P} \propto-\tau_{\ell}f(\tau_{\ell})^{2}\propto\begin{cases}-1-\tfrac{1}{12}\,(m_{a}/m_{\ell})^{2}&m_{a}\ll m_{\ell}\\ \order{\log[2](m_{a}/m_{\ell})}\left(m_{\ell}/m_{a}\right)^{2}&m_{a}\gg m_{\ell}\end{cases}\,.$$ (18)

in agreement with ref. Ferreira:2022xlw .

Therefore, for couplings to leptons significantly lighter than the ALP ($\tau_{\ell}\ll 1$), the contribution to the effective photon coupling coming from $C_{0}$ itself is negligible. In the case of derivative couplings, this leaves only a constant that results in an effective coupling of $g_{a\gamma\gamma}^{\rm eff,D}\approx\alpha C_{\ell}/(\pi\Lambda)$. Instead, for pseudoscalar couplings, the ALP-photon coupling is strongly suppressed in this limit. Conversely, in the case of couplings to leptons significantly heavier than the ALP ($\tau_{\ell}\gg 1$), the limits are effectively reversed. The derivative coupling leads to a strong suppression in $g_{a\gamma\gamma}^{\rm eff,D}$ due to cancellation between the constant and $C_{0}$, while the pseudoscalar coupling leads to a similar result as the tree-level coupling $g_{a\gamma\gamma}^{\rm eff,P}\approx-\alpha g_{a\ell}/(\pi/m_{\ell})$, albeit with the opposite sign. The ALP decay rate to two photons will feature the same limiting behaviour, see section 2.3. To summarise, the effective on-shell photon couplings are given by

$$g_{a\gamma\gamma}^{\rm eff,D} =\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\left(1-\tau_{\ell}f(\tau_{\ell})^{2}\right)\,, g_{a\gamma\gamma}^{\rm eff,P} =-\frac{\alpha}{\pi}\frac{g_{a\ell}}{m_{\ell}}\tau_{\ell}f(\tau_{\ell})^{2}\,,$$ (19)

for derivative and pseudoscalar couplings, respectively.

The loop function can also be calculated analytically for the case in which one of the photons is off-shell, i.e. eq. (12). In this case, one finds

$$2m_{\ell}^{2}C_{0}(0,t,m_{a}^{2},m_{\ell}^{2},m_{\ell}^{2},m_{\ell}^{2})=B\quantity(\frac{4m_{\ell}^{2}}{m_{a}^{2}},\frac{4m_{\ell}^{2}}{t})\,,$$ (20)

where the function $B(\tau_{1},\tau_{2})$ is defined as

$$B(\tau_{1},\tau_{2})=\frac{\tau_{1}\tau_{2}}{\tau_{1}-\tau_{2}}\quantity(f(\tau_{1})^{2}-f(\tau_{2})^{2})\,.$$ (21)

Since $t$ can also be negative, we extend the definition of $f$ to

$$f(\tau)=\begin{cases}\arcsin(\frac{1}{\sqrt{\tau}}),&\mbox{\quad if\quad}\tau\geq 1\\ \frac{\pi}{2}+\frac{i}{2}\log(\frac{1+\sqrt{1-\tau}}{1-\sqrt{1-\tau}}),&\mbox{\quad if\quad}0<\tau<1\\ \frac{i}{2}\log(\frac{\sqrt{1-\tau}+1}{\sqrt{1-\tau}-1}),&\mbox{\quad if\quad}\tau<0\end{cases}\,.$$ (22)

Note that eq. (20) is consistent with eq. (13) since in the limit of $t\to 0$

$$B(\tau_{1},\tau_{2})\xrightarrow{\tau_{2}\to\infty}-\tau_{1}f(\tau_{1})^{2}\,,$$ (23)

as $\tau_{2}f(\tau_{2})^{2}$ becomes constant. Instead, in the limit of $t\gg m_{\ell}^{2}$, one has

$$B(\tau_{1},\tau_{2})\xrightarrow{\tau_{2}\to 0}-\tau_{2}f(\tau_{2})^{2}\propto\tau_{2}\ln^{2}{\tau_{2}}\,,$$ (24)

which is particularly relevant for ALP searches with large momentum transfer.

The limit for $C_{0}$ in eq. (20) is sufficient for our phenomenological discussion since at least one of the photons will be on-shell, or at least approximately so, see section 3.2. To summarise the resulting effective coupling, one has

$$g^{\rm eff,D}_{a\gamma^{*}\gamma}(t) =\frac{\alpha}{\pi}\frac{C_{\ell}}{\Lambda}\quantity[1+B\quantity(\frac{4m_{\ell}^{2}}{m_{a}^{2}},\frac{4m_{\ell}^{2}}{t})]\,, g^{\rm eff,P}_{a\gamma^{*}\gamma}(t) =\frac{\alpha}{\pi}\frac{g_{a\ell}}{m_{\ell}}B\quantity(\frac{4m_{\ell}^{2}}{m_{a}^{2}},\frac{4m_{\ell}^{2}}{t})\,,$$ (25)

for derivative and pseudoscalar coupling, respectively. For $t\gg m_{\ell}^{2}$, the vertex function for the pseudoscalar couplings is suppressed, while it approaches a constant for the derivative coupling. In the following, we will often suppress the argument, as it is clear from the subscript whether we refer to the on- or off-shell coupling.

2.3 ALP decay rates

The decay rates for the decay of an ALP into a pair of photons or leptons are Bauer:2017ris

$$\Gamma_{a\rightarrow\gamma\gamma} =\frac{|g_{a\gamma\gamma}^{\textrm{eff}}|^{2}m_{a}^{3}}{64\pi}\,,$$ (26)

$$\Gamma_{a\rightarrow\ell\ell} =\frac{g_{a\ell}^{2}m_{a}}{8\pi}\sqrt{1-\frac{4m_{\ell}^{2}}{m_{a}^{2}}}=\frac{C_{\ell}^{2}m_{\ell}^{2}m_{a}}{8\pi\Lambda^{2}}\sqrt{1-\frac{4m_{\ell}^{2}}{m_{a}^{2}}}\,,$$ (27)

where the effective on-shell photon coupling for derivative and pseudoscalar couplings is given in eq. (19) with the limiting values in eqs. (17) and (18), respectively. For leptophilic ALPs with arbitrary couplings to the three lepton species, the total decay rate is given by

$$\Gamma_{a} =\Gamma_{a\rightarrow\gamma\gamma}+\sum_{l=e,\mu,\tau}\Gamma_{a\rightarrow\ell\ell}\,.$$ (28)

3 ALP Production at Lepton Beam Dumps

Figure 2: Schematic illustration of a beam dump. A BSM particle $X$ is produced at the target and passes through the dump. It can be detected if it decays into SM particles in the decay volume.

Beam-dump experiments are a valuable asset in the search for physics beyond the Standard Model (BSM). The general setup of a beam-dump experiment is depicted schematically in figure 2. A high-intensity beam, usually consisting of protons, electrons or muons, impinges upon a target material. Inside this target, countless Standard Model processes take place, creating particles such as pions or kaons Ertas:2021mkt . However, a small fraction of the interactions might produce a weakly interacting BSM particle $X$, for example, an ALP. In order to shield the detectors from the SM background, the beam dump is positioned downstream of the target. The dump absorbs the shower of particles, and only weakly interacting particles are able to pass. Once past the dump, the particle $X$ may decay back into Standard Model particles, which leads to a signature in the detectors placed downstream of the dump Cesarotti:2023sje .

Beam-dump experiments are well-suited to searches for particles in the mass range $\order{\text{MeV}}-\order{\text{GeV}}$ with tiny couplings to Standard Model particles. These feebly interacting particles have long lifetimes and decay lengths and are consequently difficult to detect in collider experiments. In contrast, the long decay volumes of beam-dump experiments allow for greater sensitivity to these types of particles since they are able to decay inside the experiment and be detected Cesarotti:2023sje ; Kim:2024vxg .

The most sensitive beam-dump experiments for leptophilic ALPs are the ones employing lepton beams. Concretely, we consider the E137 beam-dump experiment performed at SLAC Bjorken:1988as to search for neutral metastable particles, as well as the NA64 beam-dump experiment running at the CERN Super Proton Synchrotron. The latter performs three different types of searches:

• The standard invisible mode of the NA64 experiment differs from E137 and figure 2 in that the experiment employs an active beam dump. Instead of the detector being downstream of the target and detecting potential decay products, the target itself is the detector, and the experiment searches for missing energy Gninenko:2719646 .

• In contrast to the invisible setup, NA64 run in visible mode uses a similar search strategy to the one employed by E137 NA64:2019auh .

• The NA64$\mu$ experiment NA64:2024klw employs a muon beam to search for new particles that couple primarily to muons, such as the gauge boson of a new $\text{U(1)}_{L_{\mu}-L_{\tau}}$ gauge symmetry. Clearly, this muon mode also possesses great sensitivity to muonphilic ALP models.

The details of the different experiments are summarised in appendix A.

In the remainder of this section, we explain how the constraints for these experiments are calculated. Complementary constraints from other types of experiments will be discussed in section 4.

3.1 ALP Yield in E137

For electron beam-dump experiments like E137 or NA64 run in visible mode, the general formula for the expected number of ALPs that can be detected by the experiment is given by the expression Liu:2023bby

$$N_{a}=\frac{N_{e}X}{M_{\textrm{target}}/N_{\textrm{avo}}}\int_{E_{\textrm{min}}}^{E_{\textrm{max}}}\differential E\int_{x_{\textrm{min}}}^{x_{\textrm{max}}}\differential x\int_{0}^{T}\differential t\,I_{e}(E_{0},E,t)\derivative{\sigma}{x}e^{-L_{\textrm{sh}}\quantity(\frac{1}{l_{a}}+\frac{1}{l_{\lambda}})}\quantity(1-e^{-\frac{L_{\textrm{dec}}}{l_{a}}})\,,$$ (29)

with the number of electrons on target (EOT) $N_{e}$, the unit length of the radiation target $X$, the molar mass of the target atom $M_{\textrm{target}}$, and the Avogadro constant $N_{\textrm{avo}}$. Furthermore, $L_{\textrm{sh}}$ denotes the length of the beam dump, while $L_{\textrm{dec}}$ is the length of the decay volume, i.e. the distance from the end of the beam dump to the detector. A summary of the experimental parameters of E137 and NA64 is given in appendix A.3.

The variable $x=E_{a}/E$ defines the ratio of the energy transferred from an electron (or muon for muon beam dumps) of energy $E$ to an ALP with energy $E_{a}$. The lower integration bound for $x$ can be determined from the requirement that the ALP energy must be greater than the energy threshold of the detector $E_{\text{cut}}$, as well as the ALP rest mass. Assuming the ALP and electron/muon initial and final states are collinear Liu:2017htz , the upper bound can be approximated by requiring energy conservation in the case that the outgoing electron/muon is at rest, with the entire energy (excluding the lepton’s rest mass) being transferred to the ALP. This means the bounds for $x$ are

$$x_{\textrm{min}}=\frac{\max(m_{a}\,,\,E_{\textrm{cut}})}{E}\,,\quad x_{\textrm{max}}=\frac{E-m_{\ell}}{E}=1-\frac{m_{\ell}}{E}\,.$$ (30)

In this work, we only consider ALPs with mass below $1\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$, i.e. below the energy threshold, such that $x_{\text{min}}=E_{\text{cut}}/E$.

The upper integration boundary for the electron energy is given by the incident beam energy $E_{\textrm{max}}=E_{0}$. The energy of the electrons decreases after impacting the beam dump due to shower creation and scattering processes Tsai:1966js . This is not the case for muon beams (see section 3.4 below). The lower energy bound is given for the electron energy $E_{\textrm{min}}$ at which the bounds $x_{\textrm{min,max}}$ are equal. This integration range for the energy $E$ is, therefore, in the interval between

$$E_{\textrm{min}}=E_{\textrm{cut}}+m_{\ell}\,,\quad E_{\textrm{max}}=E_{0}\,.$$ (31)

The function $I_{e}(E_{0},E,t)$ denotes the energy distribution function for electrons with an energy $E$ at a depth $t$ inside the beam dump. The maximum depth for $t$ is given by $T=\rho L_{\textrm{tar}}/X$, where $\rho$ is the density of the target and $L_{\textrm{tar}}$ is the thickness of the target.

Figure 3: Integrated electron/positron energy distributions for E137.

In this work, we consider the primary shower of electrons $I_{e}^{(1)}$, as well as the secondary shower of electrons and positrons $I_{e,p}^{(2)}$, which are created through pair production from photons in the shower Tsai:1966js . We use the expressions for $I_{e}^{(1,2)}$ derived in ref. Tsai:1966js . However, we do not use the simplified expressions, which assume that $E_{0}-E\ll E_{0}$, but instead use the full expressions for the distributions. Since the electron energy distribution is the only function of $t$ in the integral (29), we define the integrated energy distribution

$$I_{\textrm{int}}(E_{0},E,T)\equiv\int_{0}^{T}\differential t\,I_{e}(E_{0},E,t)\,.$$ (32)

The integrated distributions for first and second-generation electrons/positrons in E137 are shown in figure 3. In addition, we indicate the distribution of first-generation electrons obtained when using the approximation $E_{0}-E\ll E_{0}$. The approximate expression differs substantially from the exact first-generation distribution. However, it coincidentally matches the more complete result, which includes second-generation effects, quite well.

The absorption of ALPs by electrons in the target is described by $l_{\lambda}$, which leads to exponential attenuation. This is relevant for thick targets but can be neglected here Liu:2023bby . The ALP decay length $l_{a}$ depends on the couplings as well as the energy and mass of the ALP. It can be written as

$$l_{a}=\frac{p_{a}}{m_{a}}\frac{1}{\Gamma_{a}}=\frac{\sqrt{E_{a}^{2}-m_{a}^{2}}}{m_{a}}\frac{1}{\Gamma_{a}}=\frac{\sqrt{x^{2}E^{2}-m_{a}^{2}}}{m_{a}}\frac{1}{\Gamma_{a}}\,.$$ (33)

The first exponential function in eq. (29) can be understood as the probability that the ALP only decays after leaving the beam dump, such that the decay products do not get absorbed. Conversely, the last term describes the probability that the ALP does not decay behind the detector, i.e. that it decays inside the decay volume and can be detected.

With the simplifications discussed above, the general formula we use for the expected number of detected ALPs is

$$N_{a}=\frac{N_{e}X}{M_{\textrm{target}}/N_{\textrm{avo}}}\int_{E_{\textrm{min}}}^{E_{0}}\differential E\int_{x_{\textrm{min}}}^{x_{\textrm{max}}}\differential x\,I_{\textrm{int}}(E_{0},E,T)\derivative{\sigma}{x}e^{-\frac{L_{\textrm{sh}}}{l_{a}}}\quantity(1-e^{-\frac{L_{\textrm{dec}}}{l_{a}}})\,.$$ (34)

The only missing ingredient is the differential cross section $\mathrm{d}\sigma/\mathrm{d}x$ for the production of ALPs, which will be discussed next.

3.2 Dark Bremsstrahlung and Primakoff Production

The main ALP production mechanisms in beam-dump experiments are Dark Bremsstrahlung (DB) and Primakoff production, shown in the top row of figure 4.

Figure 4: Top row: Feynman diagrams contributing to the process $eN\rightarrow eNa$. The diagrams correspond from left to right to the DB $s$-channel process, the DB $u$-channel process and the Primakoff process. Bottom row: Corresponding Feynman diagrams contributing to the process $e\gamma\rightarrow ea$. Note that the lepton loop has been replaced with the effective Primakoff coupling.

In the DB diagrams, the leptophilic ALP couples to the beam electron/muon at tree level, while for the Primakoff process, it only couples to the photons via the effective loop coupling.

Due to the high energy of the beam compared to the ALP and lepton masses, the photon that is exchanged between the nucleus and the electron/muon is approximately on-shell Kim:1973he . As a result, the nucleus can be replaced by an effective flux of photons, which is described by the Improved Weizsäcker-Williams (IWW) approximation. The IWW approximation allows the cross section to be calculated using the reduced $2\to 2$ scattering amplitude $\overline{\left|\mathcal{M}\right|}^{2}$. Thus, for a lepton scattering with a nucleus

$$\ell^{-}(p)+N(P_{i})\rightarrow\ell^{-}(p^{\prime})+N(P_{f})+a(k)\,,$$ (35)

the simplified photoproduction processes are

$$\ell^{-}(p)+\gamma(q)\rightarrow\ell^{-}(p^{\prime})+a(k)\,.$$ (36)

The simplified DB and Primakoff diagrams are shown in the bottom row of figure 4.

The differential cross section for the tree-level DB process can be calculated analytically (see appendix B.1) and is found to be

$$\begin{split}\derivative{\sigma_{\textrm{DB}}}{x}&=\frac{\alpha}{(4\pi)^{2}}\sqrt{x^{2}-\frac{m_{a}^{2}}{E^{2}}}\chi\frac{1-x}{x}\\ &\times\frac{e^{2}E^{2}g_{a\ell}^{2}\theta_{a,\text{max}}^{2}x^{5}}{3(1-x)\left(m_{a}^{2}(1-x)+m_{\ell}^{2}x^{2}\right)^{2}\left(m_{a}^{2}(1-x)+x^{2}\left(E^{2}\theta_{a,\text{max}}^{2}+m_{\ell}^{2}\right)\right)^{3}}\\ &\times\Bigl{[}3m_{a}^{4}(1-x)^{2}x^{2}\left(E^{2}\theta_{a,\text{max}}^{2}+3m_{\ell}^{2}\right)+3x^{6}\left(E^{2}\theta_{a,\text{max}}^{2}m_{\ell}+m_{\ell}^{3}\right)^{2}\Bigr{.}\\ \Bigl{.}&\quad+m_{a}^{2}(1-x)x^{4}\left(2E^{4}\theta_{a,\text{max}}^{4}+9E^{2}\theta_{a,\text{max}}^{2}m_{\ell}^{2}+9m_{\ell}^{4}\right)+3m_{a}^{6}(1-x)^{3}\Bigr{]}\,,\end{split}$$ (37)

where $\theta_{a,\text{max}}$ denotes the largest angle relative to the beam axis that the ALP can have in order to be detected. The Primakoff process and the interference term involve the momentum-dependent effective coupling to photons, so we do not obtain a simple analytical expression for the corresponding differential cross sections. Instead, we evaluate them numerically.

In order to analyse the mass and energy dependence of the cross sections, we plot the contributions to $\sigma$ (normalised by the coupling) coming from DB, Primakoff, and interference processes in E137 in figure 5. We assume that the ALP couples exclusively to electrons and compare the normalised cross sections for two masses: $m_{a}=$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ and $m_{a}=$450\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$.

Figure 5: Normalised cross sections for non-vanishing ALP coupling to electrons as a function of the electron energy $E$. The ALP mass is $m_{a}=$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ on the left, $m_{a}=$450\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ on the right. The solid lines represent derivative couplings, while the dashed lines represent pseudoscalar couplings.

For $m_{a}=$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, the DB cross section is clearly the main contributor to the total cross section. This is to be expected, considering these interactions are mediated by tree-level interactions between the ALPs and electrons. However, it would be wrong to conclude that the Primakoff process does not give important contributions elsewhere, as can be seen for $m_{a}=$450\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. For this mass, Primakoff production outweighs DB production, assuming the ALP exhibits a derivative coupling to electrons (while for pseudoscalar coupling the Primakoff contribution is again negligible). It should be noted that comparing $\sigma$ is only illustrative, as the actual integration over $x$ weighs $\differential\sigma/\differential x$ with the exponential functions given at the end of eq. (34). However, it is clear that while the Primakoff production of leptophilic ALPs occurs only via an effective loop coupling, it can still have an impact on constraints that can be set on the ALP mass and couplings.

Finally, there are still the contributions arising from secondary positrons created in the shower. Both the DB and Primakoff processes give exactly the same contributions when the electron line in the Feynman diagrams is replaced with a positron line. The only difference between the contributions are the electron and positron energy distribution functions $I_{e,p}$. This means that instead of considering the electron and positron production processes separately, we summarise them together by using the total electron-positron distribution function

$$I_{\textrm{int}}(E_{0},E,T) =\int_{0}^{T}\differential t\quantity(I_{e}^{(1)}(E_{0},E,t)+I_{e}^{(2)}(E_{0},E,t))+\int_{0}^{T}\differential tI_{p}^{(2)}(E_{0},E,t)$$ (38)

$$=\int_{0}^{T}\differential t\quantity(I_{e}^{(1)}(E_{0},E,t)+2I_{e}^{(2)}(E_{0},E,t))\,.$$ (39)

For muonphilic ALPs, there can also be ALP bremsstrahlung from secondary muons produced in the electromagnetic shower Marsicano:2018vin . Comparing our results to the ones in ref. Marsicano:2018vin , we find that the loop-induced Primakoff process dominates ALP production for muonphilic ALPs. Therefore, we do not include the contributions of secondary muons in the present work.

3.3 Annihilation Processes

In addition to the contributions from DB and Primakoff production, we also consider ALP production from processes in which a positron created in a secondary shower annihilates with an electron inside the beam dump. Since only positrons from the shower contribute to these processes, the first important change compared to DB and Primakoff production is that the integrated energy distribution appearing in eq. (34) only contains the contributions from positrons, i.e.

$$I_{\textrm{int}}(E_{0},E,T)=\int_{0}^{T}\differential t\,I_{p}^{(2)}(E_{0},E,t)\,.$$ (40)

If the annihilation of the positron and electron has a centre-of-mass energy that is exactly $\sqrt{s}=m_{a}$, an ALP can be produced resonantly:

$$e^{-}(p_{1})+e^{+}(p_{2})\rightarrow a(k)\,.$$ (41)

The Feynman diagram for this process is shown in figure 6.

Figure 6: Resonant ALP production through the process $e^{-}e^{+}\rightarrow a$.

The matrix element for this process is simply

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{R}}}}^{2} =\frac{g_{ae}^{2}s}{2}\,,$$ (42)

and the differential cross section is given by

$$\derivative{\sigma_{\textrm{R}}}{x} =\frac{\pi g_{ae}^{2}}{4m_{e}\sqrt{1-\frac{4m_{e}^{2}}{m_{a}^{2}}}}\,\delta\quantity(E-\quantity(\frac{m_{a}^{2}}{2m_{e}}-m_{e}))\delta\quantity(x-\quantity(1+\frac{m_{e}}{E}))\,.$$ (43)

The first delta distribution in eq. (43) is equivalent to the condition $s=m_{a}^{2}$. After inserting this cross section into eq. (34), the two delta distributions cancel out the integrations, leaving behind

$$N_{a,\text{res}}=\frac{N_{e}X}{M_{\textrm{target}}/N_{\textrm{avo}}}\frac{\pi g_{ae}^{2}}{4m_{e}\sqrt{1-\frac{4m_{e}^{2}}{m_{a}^{2}}}}I_{\textrm{int}}(E_{0},E_{*},T)e^{-\frac{L_{\textrm{sh}}}{{l_{a}}_{*}}}\quantity(1-e^{-\frac{L_{\textrm{dec}}}{{l_{a}}_{*}}})\,.$$ (44)

Here, $E_{*}$ and ${l_{a}}_{*}$ are the positron energy and decay length, respectively, after applying the conditions given by the two delta distributions in eq. (43). The mass range in which resonant production can be detected is bounded from below by the fact that the ALP must have an energy larger than the cutoff energy of the detector and from above by the maximum possible positron energy, i.e.

$$E_{a}=E+m_{e}=\frac{m_{a}^{2}}{2m_{e}} >E_{\textrm{cut}}\,,$$ (45)

$$E=\frac{m_{a}^{2}}{2m_{e}}-m_{e} <E_{0}\,.$$ (46)

This means resonantly produced ALPs can only be detected when the ALP mass is in the interval $\sqrt{2m_{e}E_{\textrm{cut}}}<m_{a}<\sqrt{2m_{e}(E_{0}+m_{e})}$.

In addition to the resonant production discussed above, electron-positron annihilation can also proceed non-resonantly. The corresponding cross sections are discussed in appendix B.2.

3.4 ALP Yields in NA64

While eq. (34) holds for E137 and the visible mode setup of NA64, there are some relevant changes when calculating the expected yield of ALPs for the invisible setup as well as the muon setup of NA64.

NA64 in invisible mode

Since the invisible setup of NA64 is an active beam dump that measures the missing energy of an event instead of detecting particles downstream of the dump, one can consider the decay volume to be infinitely large, i.e. one can set $L_{\textrm{dec}}\rightarrow\infty$. This simplifies the equation for $N_{a}$ to

$$N_{a}=\frac{N_{e}X}{M_{\textrm{target}}/N_{\textrm{avo}}}\int_{E_{\textrm{min}}}^{E_{0}}\differential E\int_{x_{\textrm{min}}}^{x_{\textrm{max}}}\differential xI_{\textrm{int}}(E_{0},E,T)\derivative{\sigma}{x}e^{-\frac{L_{\textrm{sh}}}{l_{a}}}\,,$$ (47)

where the various quantities that enter this equation are given in appendix A.3. In principle, it would also be possible to simplify the expression for $\mathrm{d}\sigma/\mathrm{d}x$, since there is no longer a restriction on the emission angle $\theta_{a}$ of the ALP. Since the main contribution to the cross section comes from small angles Dusaev:2020gxi , one can safely extend the integration to large angles of order unity. Concretely, in analogy to E137, we introduce an arbitrary cutoff for the angle based on the geometry of the detector given by

$$\theta_{a,\text{max}}=\frac{R_{\textrm{sh}}}{L_{\textrm{sh}}}\,,$$ (48)

where $R_{\textrm{sh}}\approx$60\text{\,}\mathrm{cm}$$ is the approximate radius of the detector, in this case of the HCALs.

NA64$\mu$

The general setup of the NA64$\mu$ experiment is similar to that of NA64 run in invisible mode, with it being an active beam dump. This means that the formula (47) applies to NA64$\mu$ as well, with the caveat that the electron lines in the DB and Primakoff diagrams in figure 4 need to be replaced with muon lines, and the EOT $N_{e}$ should be replaced with the muons on target (MOT) $N_{\mu}$. Since there are no muons for the anti-muons in the beam to annihilate with, the contribution discussed in section 3.3 is absent. Finally, an important difference is the energy distribution function of the muons compared to the electrons. Since muons are much more massive than electrons, they lose significantly less energy due to bremsstrahlung compared to electrons, allowing them to easily penetrate the beam dump NA64:2024klw . This means the muon’s energy distribution function is independent of the depth $t$. Therefore, we can use the thin-target approximation to calculate the distribution function for muons, which is given by Tsai:1966js

$$I_{\mu}(E_{0},E,t)=\delta(E_{0}-E)\quad\Rightarrow\quad I_{\textrm{int}}(E_{0},E,T)=T\delta(E_{0}-E)\,.$$ (49)

This means the expected yield of ALPs simplifies to

$$N_{a}=\frac{N_{\mu}X}{M_{\textrm{target}}/N_{\textrm{avo}}}T\int_{x_{\textrm{min}}}^{x_{\textrm{max}}}\differential x\derivative{\sigma}{x}e^{-\frac{L_{\textrm{sh}}}{l_{a}}}\,,$$ (50)

where it is implied that $E=E_{0}$ in the expressions of $x_{\textrm{min,max}}$, $l_{a}$, and the cross section. It is important to note that this makes $N_{a}$ very sensitive to the length $T$ of the beam dump. This is in contrast with the other electron beam dumps, in which the electron shower attenuates after several radiation lengths.

4 Other Constraints

In this section we review other experimental and theoretical constraints of relevance for leptophilic ALPs.

4.1 Perturbative unitarity

A simple but important constraint on the ALP coupling to leptons can be obtained by considering the scattering process $\ell^{+}\ell^{-}\to\ell^{+}\ell^{-}$ via $s$-channel ALP exchange. Even though the ALP interaction is non-renormalisable, the matrix element for this process becomes independent of the centre-of-mass energy for $\sqrt{s}\gg m_{a}$. The requirement of unitarity of the $J=0$ partial wave $a_{0}$, i.e. $|\text{Re}\ a_{0}|<1/2$, implies Cornella:2019uxs

$$\frac{C_{\ell}}{\Lambda}=\frac{g_{a\ell}}{m_{\ell}}<\frac{1}{m_{\ell}}\sqrt{\frac{8\pi}{3}}\,.$$ (51)

4.2 LEP bounds on exotic $Z$ boson decays

Ref. Jaeckel:2015jla studied constraints on ALP couplings to photons from LEP measurements on the $Z$ pole. Here, we apply this analysis to the present setup, where photon couplings are induced via lepton loops. As we will see, these loops are strongly suppressed by $m_{\ell}^{2}/M_{Z}^{2}$ for the case of pseudoscalar ALP couplings, meaning relevant constraints from LEP only arise for derivative couplings.

There are two processes of interest: $e^{+}e^{-}\to Z\to a\gamma$ via on-shell production of a $Z$ boson, and $e^{+}e^{-}\to\gamma^{*}\to a\gamma$ via an off-shell photon at the $Z$ peak. For light ALPs with mass below $m_{\pi}$, the photons from the radiative ALP decay are very collimated, so that both processes look like the forbidden SM decay $Z\to 2\gamma$ Chala:2015cev . For $Z$-boson decays we follow ref. Jaeckel:2015jla and use the bound

$${\rm BR}(Z\to a\gamma)\times{\rm BR}(a\to\gamma\gamma)\times P(l_{a})\leq 5.2\times 10^{-5}\,,$$ (52)

which is inferred from the 95% confidence level LEP limit ${\rm BR}(Z\to\pi^{0}\gamma)\leq 5.2\times 10^{-5}$ L3:1995nbq . Here, we conservatively include only the fraction ${\rm BR}(a\to\gamma\gamma)$ of ALPs decaying into two photons since it is unclear whether a leptonically decaying ALP could mimic a $\pi^{0}$ decay when strongly boosted. Moreover, the factor

$$P(l_{a})=1-e^{-\tfrac{L}{l_{a}}}\,,$$ (53)

gives the probability that the ALP decays inside the decay volume and can be detected. Here, $l_{a}$ is the ALP decay length in the laboratory frame (see section 3.1), and we take $L=10\,\mathrm{cm}$. The above limits are used for $m_{a}\leq m_{\pi^{0}}=135{\,\rm MeV}$. For larger ALP masses up to 1 GeV we use the limits obtained in ref. Jaeckel:2015jla from simulating angular distributions for the three-photon decay ($e^{+}e^{-}\to Z\to a+\gamma\to 3\gamma$), which were compared bin by bin to the distributions given in ref. L3:1995nbq . The resulting bounds are weaker than those for $m_{a}\leq m_{\pi^{0}}$ by about 30%.

These constraints can be directly interpreted as a bound on lepton couplings, using the decay rates for derivative (D) and pseudoscalar (P) couplings Bauer:2021mvw

$$\Gamma(Z\to a\gamma)_{\rm D} =\frac{m_{Z}^{3}}{96\pi^{3}}\frac{\alpha^{2}}{s_{\rm w}^{2}c_{\rm w}^{2}}\frac{|C_{\ell}|^{2}}{\Lambda^{2}}\left(s_{\rm w}^{2}-1/4\right)^{2}\left|1+B(\tau_{\ell},\tau_{Z})\right|^{2}\left(1-\frac{m_{a}^{2}}{m_{Z}^{2}}\right)^{3}\,,$$

$$\Gamma(Z\to a\gamma)_{\rm P} =\frac{m_{Z}^{3}}{96\pi^{3}}\frac{\alpha^{2}}{s_{\rm w}^{2}c_{\rm w}^{2}}\frac{|g_{a\ell}|^{2}}{m_{\ell}^{2}}\left(s_{\rm w}^{2}-1/4\right)^{2}\left|B(\tau_{\ell},\tau_{Z})\right|^{2}\left(1-\frac{m_{a}^{2}}{m_{Z}^{2}}\right)^{3}\,,$$ (54)

with $s_{\text{w}}$ denoting the sine of the Weinberg angle, $\tau_{Z}=4m_{\ell}^{2}/M_{Z}^{2}$ and $B(\tau_{1},\tau_{2})$ given in eq. (21). Since here we consider $m_{a}\ll m_{Z}$, one has $\tau_{Z}\ll\tau_{\ell}$, so that to good approximation

$$B(\tau_{\ell},\tau_{Z}) =\tau_{Z}\left(f^{2}(\tau_{\ell})-f^{2}(\tau_{Z})\right)\ll 1\,.$$ (55)

Constraints on the pseudoscalar scenario are thus much weaker than those in the derivative scenario. In fact, they are weaker than the perturbative unitarity bounds for all leptons and, therefore, irrelevant.

Figure 7: Constraints on ALP derivative couplings from LEP searches for $Z\to 3\gamma$ from direct ALP production either via on-shell or off-shell photons (purple), see text for details. Also shown are the limits from perturbative unitarity (gray), which always dominate over analogous constraints on pseudoscalar couplings. Note that for $m_{a}\leq m_{\pi^{0}}$, the bounds are slightly stronger than above, due to highly collimated photons looking like the forbidden $Z\to\gamma\gamma$ decay, see ref. Jaeckel:2015jla for more details.

For the derivative scenario, the constraints are more relevant but still rather weak due to the small $Z$ boson couplings to SM leptons resulting from the accidental cancellation $s_{w}^{2}\approx 1/4$. As we will see below, the LEP limits from $Z$ decays are coincidentally of the same magnitude as the limits from ALP production via off-shell photons.

Indeed, the same signature as $Z\to a\gamma$ is generated from ALP production via off-shell photons $e^{+}e^{-}\to\gamma^{*}\to a\gamma$ with $s=M_{Z}^{2}$. Therefore, the cross-section is limited by

$$\sigma({e^{+}e^{-}\to\gamma^{*}\to a\gamma})|_{s=M_{Z}^{2}} \leq\sigma_{\rm limit}({e^{+}e^{-}\to Z\to a\gamma})|_{s=M_{Z}^{2}}$$

$$=\frac{12\pi}{m_{Z}^{2}}{\rm BR}(Z\to e^{+}e^{-}){\rm BR}_{\rm limit}(Z\to a\gamma)$$

$$=\frac{7.9\times 10^{-9}}{{\,\rm GeV}^{2}}\times\frac{1}{{\rm BR}(a\to\gamma\gamma)}\times\frac{1}{P(l_{a})}\,,$$ (56)

where we have used the narrow-width approximation in the second line and eq. (52) in the third line. The cross section at the $Z$ pole can be calculated using the general vertex function with a single off-shell photon in the derivative and pseudoscalar scenario, respectively (cf. eq. (25))

$$\Gamma^{\mu\nu,D}_{a\gamma^{*}\gamma} =i\frac{\alpha C_{\ell}}{\pi\Lambda}\epsilon^{\mu\nu\alpha\beta}q_{\alpha}k_{\beta}\left[1+B(\tau_{\ell},\tau_{Z})\right]\,,$$

$$\Gamma^{\mu\nu,P}_{a\gamma^{*}\gamma} =i\frac{\alpha g_{a\ell}}{\pi m_{\ell}}\epsilon^{\mu\nu\alpha\beta}q_{\alpha}k_{\beta}B(\tau_{\ell},\tau_{Z})\,,$$ (57)

where the external on-shell photon has outgoing momentum $q$ and polarisation index $\mu$, $k$ is the outgoing momentum of the axion, and $p=k+q$ is the inflowing momentum of the off-shell photon with polarisation index $\nu$, with $p^{2}=s=M_{Z}^{2}$. The total cross-section is given by

$$\sigma({e^{+}e^{-}\to\gamma^{*}\to a\gamma})_{\rm D} =\frac{\alpha^{3}}{24\pi^{2}}\frac{|C_{\ell}|^{2}}{\Lambda^{2}}\frac{\left(1-\frac{m_{a}^{2}}{m_{Z}^{2}}\right)^{3}\left(1+\frac{2m_{\ell}^{2}}{m_{Z}^{2}}\right)}{\sqrt{1-4m_{\ell}^{2}/M_{Z}^{2}}}\left|1+B(\tau_{\ell},\tau_{Z})\right|^{2}\,,$$

$$\sigma({e^{+}e^{-}\to\gamma^{*}\to a\gamma})_{\rm P} =\frac{\alpha^{3}}{24\pi^{2}}\frac{|g_{a\ell}|^{2}}{m_{\ell}^{2}}\frac{\left(1-\frac{m_{a}^{2}}{m_{Z}^{2}}\right)^{3}\left(1+\frac{2m_{\ell}^{2}}{m_{Z}^{2}}\right)}{\sqrt{1-4m_{\ell}^{2}/m_{Z}^{2}}}\left|B(\tau_{\ell},\tau_{Z})\right|^{2}\,,$$ (58)

Again, constraints on the lepton couplings in the pseudoscalar scenario are weaker than those from perturbative unitarity, such that constraints from ALP production at LEP are relevant only in the derivative scenario. They are accidentally of the same size as the constraints from $Z$ boson decays, since

$$\frac{\Gamma_{Z}}{12\pi M_{Z}{\rm BR}(Z\to e^{+}e^{-})}\frac{e^{2}s_{\rm w}^{2}c_{\rm w}^{2}}{(s_{\rm w}^{2}-1/4)^{2}}\simeq 0.94\,.$$ (59)

The resulting constraints on the lepton couplings in the derivative scenario are shown in the ($m_{a}-C_{\ell}/\Lambda$ plane) in figure 7 for all three individual lepton couplings as well as the universal case, with the limits from perturbative unitarity shown in grey.

4.3 B factories

In ref. BaBar , the BaBar collaboration reported results from a search for dark photons $A^{\prime}$ produced in the reaction $e^{+}e^{-}\to\gamma A^{\prime}$, $A^{\prime}\to\ell^{+}\ell^{-}$ with $\ell=e,\mu$ in terms of a 90% CL upper bound on the kinetic mixing parameter $\epsilon$ as a function of the dark photon mass $m_{A^{\prime}}$. We reinterpret these results in terms of the couplings of leptophilic ALPs by comparing the cross section for dark photon production with the one for non-resonant production of ALPs, see appendix B.2.

Since the BaBar search focuses on dark photon decays to light charged leptons, the bound does not apply to models of ALPs coupled exclusively to tau leptons. For ALPs coupled exclusively to electrons, we only reinterpret the bounds up to the muon production threshold, $m_{a}<2m_{\mu}$. Similarly, for ALPs coupled exclusively to muons, we only reinterpret the BaBar bound for $m_{a}>2m_{\mu}$. Note that even though there is no tree-level coupling to electrons in this case, the ALP can still be produced through the s-channel process $e^{-}e^{+}\to\gamma^{*}\to a\gamma$ via the effective photon coupling induced through the muon loop.

In ref. Belle-II:2024wtd , the Belle II collaboration reported results from their search for a muonphilic scalar $S$ via the process $e^{-}e^{+}\to\mu^{-}\mu^{+}S$ with $S\to\mu^{-}\mu^{+}$. The Belle II collaboration assumed an interaction Lagrangian of the form $\mathcal{L}_{\text{int}}\supset-ig_{a\mu}a\overline{\psi}_{\mu}\psi_{\mu}$, which allows us to directly use their results, neglecting the small differences of scalars compared to pseudoscalars. Moreover, this constraint applies to both the derivative and pseudoscalar coupling cases since it involves a tree-level coupling between the muon and the ALP.

4.4 Meson and $W^{+}$ boson decays

Ref. Altmannshofer:2022ckw pointed out the possibility of searching for leptophilic ALPs using rare decays such as $\pi^{+}\to e^{+}\nu a$ or $K^{+}\to e^{+}\nu a$ or $W^{+}$ boson decays such as $W^{+}\to e^{+}\nu a$. The resulting constraints can be directly applied to our setup, noting that the derivative scenario corresponds to the weak-violating ALP in ref. Altmannshofer:2022ckw and the pseudoscalar scenario to the weak-conserving ALP.

4.5 Anomalous cooling of SN 1987A

A burst of ALPs can be produced by the electrons, muons or photons that are present in large numbers in the hot and dense plasmas of core-collapse supernovae. In 1987, one of these explosions, named SN 1987A, was observed in the Large Magellanic Cloud, around 50 kpc away from Earth. This is the only supernova to date for which the accompanying neutrino emission was measured by experiments here on Earth. The duration of this neutrino burst is in overall agreement with simulations Fiorillo:2023frv , which can be used to constrain ALP models. If any exotic particle would be produced in the central neutrinosphere of the nascent proto-neutron star in the supernova core and was coupled so weakly that it escaped without re-depositing its energy, it would provide an additional cooling channel that is not present in simulations. Following approximate, order-of-magnitude physical arguments supported by early simulations of the explosion, the additional cooling power cannot be larger than the neutrino luminosity — otherwise the duration of the neutrino burst would be significantly shortened, contrary to observations Raffelt:1996wa ; Raffelt:2006cw .

For derivative electron couplings, ref. Ferreira:2022xlw takes into account the loop-induced ALP-photon interactions, while for pseudoscalar couplings between electrons and ALPs, we show the tree-level-only result from ref. Ferreira:2022xlw as a gray dashed line in the right panel of figure 10. However, omitting the loop-induced photon coupling in the latter case is not a strictly conservative assumption since ALPs are not only produced through interactions with photons but can also be reabsorbed before leaving the neutrinosphere, leading to a weaker constraint for comparably large couplings. Hence, these constraints should be viewed as an approximation indicating the order-of-magnitude at which SN constraints become relevant.

For interactions with muons, the cooling constraint on the derivative coupling was presented in ref. Croon:2020lrf . The constraint is, in principle, different for pseudoscalar couplings since the differing loop-induced photon coupling significantly affects the upper end of the bound for $5{\,\rm MeV}\lesssim m_{a}\lesssim 2m_{\mu}\simeq 211{\,\rm MeV}$. For this case, ref. Caputo:2021rux gives a rough estimate of the cooling constraint in the “axion-sphere” approximation. They also derive an estimated bound from the ALPs transferring energy into the mantle of the progenitor star, which would contribute to the explosion energy of the circumstellar material. For typical supernovae such as SN 1987A, this explosion energy is observed to be around $10^{51}$ erg. For $m_{a}\lesssim 10$ MeV, ref. Caputo:2021rux estimates that ALPs decaying into photons in the progenitor mantle yield a contribution larger than the observed value unless $g_{a\mu}/m_{\mu}\gtrsim 0.91\cdot 10^{-3}~{}{\rm GeV}^{-1}$. We show this constraint in the right panel of figure 11.

For ALPs coupled only to tau leptons, one can safely integrate them out, given that the tau mass is much larger than typical temperatures reached in the proto-neutron star. This leaves only a constant on-shell coupling to photons for the pseudoscalar case, which is constrained by the usual analyses, while for the derivative case, the coupling is suppressed by $m_{a}^{2}/m_{\tau}^{2}$. In the universal coupling case, production from electrons, muons, and the loop-induced photon interaction (also receiving a contribution from $g_{a\tau}$) have to be taken into account, and the same holds true for reabsorption and decay events. Thus, neither of the bounds on electron or muon couplings are accurate estimates for this case.

4.6 Proton beam-dump experiments and forward-physics detectors

Proton beam dump experiments provide sensitive probes of the effective ALP-photon coupling  Dobrich:2015jyk ; Dobrich:2019dxc . A detailed reinterpretation of these constraints in the context of leptophilic ALPs is not straightforward since the typical momentum transfer $t$ from the target nucleus to the ALP can be comparable to $m_{\ell}^{2}$, such that the full momentum dependence of the lepton loop needs to be included for an accurate estimate.

To obtain a first estimate of the expected constraints, we implemented the constant piece of the loop-induced photon coupling $g_{a\gamma\gamma}^{\text{eff}}=\alpha C_{\ell}/(\pi\Lambda)$ in the public code ALPINIST Jerhot:2022chi , which allows for a reinterpretation of proton-beam dump data for ALPs with different coupling structures. The resulting constraints are found to be weaker than the ones from E137 and NA64. Given the results from ref. Rella:2022len for muonphilic scalars, we anticipate the same conclusion to hold also for muonphilic ALPs, which can be produced via Bremsstrahlung from secondary muons. Since a more accurate calculation is not expected to give significantly stronger bounds, we refrain from a detailed study and omit these constraints from our analysis. For the same reason, we do not attempt to reinterpret the recent search for ALPs coupled to photons at FASER FASER:2024bbl .

4.7 Muon and Electron $g-2$

Figure 8: One-loop contributions (left) and two-loop contributions (right) of muonphilic ALPs to $(g-2)_{\mu}$.

Electrophilic or muonphilic ALPs contribute to the anomalous magnetic moment of the electron or the muon, $a_{\ell}\equiv(g-2)_{\ell}/2$ with $\ell=e,\mu$, respectively. The main contributions for the muonic case only are shown in figure 8.

At the one-loop level, there is no difference between derivative and pseudoscalar couplings due to the direct coupling of the ALP to the external charged lepton lines. Conversely, the two-loop contribution, which features an effective ALP-photon coupling, distinguishes the two coupling types. The expression for the one-loop contribution is given in ref. Bauer:2017ris , while the expressions for the two-loop contributions can be found in refs. Buttazzo:2020vfs ; Armando:2023zwz and Buen-Abad:2021fwq for pseudoscalar and derivative couplings, respectively. Note that, in the latter case, the result depends on the UV cutoff scale $\Lambda$, which we will set here to $\Lambda\sim$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$.

The theoretical situation regarding the muon anomalous magnetic moment is currently under scrutiny, see e.g. ref. Davies:2025pmx and references therein. Indeed, up to a few years ago, the theoretical consensus was a value Aoyama:2020ynm that currently shows a negative 5.2$\sigma$ discrepancy with the current experimental world average Muong-2:2023cdq , a discrepancy that would only be exacerbated by ALPs contributions. However, the latest lattice determinations of the so-called leading order hadronic vacuum polarisation (LOHVP) contribution increase the theoretical prediction of $a_{\mu}$, effectively resulting in a compatibility at the 0.9$\sigma$ level Borsanyi:2020mff ; Boccaletti:2024guq . Yet these determinations are incompatible with the previous LOHVP determinations at the 4$\sigma$ level. While waiting for the resolution of this theoretical controversy, we decided to adopt a conservative approach here and include in our results only a line corresponding to the ALP parameter values such that its presence would decrease $a_{\mu}^{\textrm{theory}}$ by $\delta a_{\mu}=-10^{-9}$, i.e. corresponding to a negative $2\sigma$ shift from the current experimental value, when current theoretical and experimental errors are combined. Nevertheless, this result must not be considered as a binding constraint, as nothing prevents the presence of additional new processes contributing to $a_{\mu}$ and counteracting the ALP effect.

A similar conservative approach is required as well for the electron anomalous magnetic moment, albeit for different reasons. In this case, the discrepant determinations of the LOHVP do not play a relevant role because this contribution is negligible due to its electron mass suppression. However, at the current level of accuracy of $10^{-13}$, the theoretical prediction is sensitive to different measurements of the fine-structure constant employed in its determination, where a 5.4$\sigma$ difference has been observed in its two latest experimental determinations Parker_2018 ; Morel:2020 . The latest experimental determination of $a_{e}$, measured with a precision similar to the theoretical uncertainty, sits between the two theoretical determinations Fan:2022eto . Therefore, we adopt a similar approach to the muon case, showing in our results the line corresponding to ALP parameters values such that they would decrease $a_{e}^{\textrm{theory}}$ by $10^{-13}$. Once again, this constraint does not have to be considered binding due to possible additional contributions to this channel.

4.8 ALP lifetime

Axion-like particles decaying into leptons and photons in the early universe may change its expansion rate, the temperature of the electron-photon bath relative to the one of neutrinos and the primordial element abundances Depta:2020zbh . While a detailed calculation of these constraints is challenging, in particular for the muonic decay mode, it is expected that ALP decays are harmless for cosmology if the ALP lifetime is shorter than $0.1\,\mathrm{s}$. To give the reader a feeling for how relevant this consideration is, we indicate a line labelled “$\tau_{a}=$0.1\text{\,}\mathrm{s}$$” in our plots. This line should, however, not be interpreted as a hard constraint since additional effects can change the cosmological evolution, and the resulting constraints may be weaker.

5 Results

In this section, we present and discuss the bounds on the leptophilic ALP models obtained from the beam-dump experiments E137 and NA64 discussed in section 3 as well as the other constraints discussed in section 4. In order to remain agnostic about the UV completion of this scenario, we analyse different coupling structures independently:

• Scenario E: Non-vanishing electron coupling $\hat{g}_{ae}\neq 0$  ,

• Scenario M: Non-vanishing muon coupling $\hat{g}_{a\mu}\neq 0$  ,

• Scenario T: Non-vanishing tau coupling $\hat{g}_{a\tau}\neq 0$  ,

• Scenario U: Universal lepton coupling $\hat{g}_{ae}=\hat{g}_{a\mu}=\hat{g}_{a\tau}$  ,

where we have introduced the definitions $\hat{g}_{a\ell}\equiv g_{a\ell}/m_{\ell}$ or $\hat{g}_{a\ell}\equiv C_{\ell}/\Lambda$ for the pseudoscalar or the derivative case, respectively. For each scenario, we consider both derivative and pseudoscalar couplings and set all other ALP couplings to zero. All bounds are calculated at $90\%$ confidence level, except the ones from LEP and $W^{+}\to e^{+}\nu a$, which are at 95% confidence level.

5.1 Effect of the coupling type on the E137 bound

Before presenting our full results, we compare the effect of the choice of coupling type on the bounds derived from the E137 experiment. The comparisons of the bounds for all scenarios are shown in figure 9.

In this figure, aside from comparing the derivative coupling given in eq. (1) and the pseudoscalar coupling in eq. (4), we also show a comparison with a simplifying approximation, in which the effective off-shell coupling of the ALP to photons is set to the constant on-shell value, in order to judge the relevance of the momentum dependence of the effective photon coupling. We approximate the couplings according to the limits of the derivative and pseudoscalar couplings given in eqs. (17) and (18), respectively. Since the ALP mass is larger than that of the electron in the most relevant regions of parameter space, we assume $\tau_{e}=4m_{e}^{2}/m_{a}^{2}\ll 1$. For the muon and tau contribution, on the other hand, we can approximate $\tau_{\mu,\tau}\gg 1$ in the parameter regions probed by E137. In the on-shell approximation, we thus approximate the derivative coupling by the constant values

$$g_{a\gamma^{*}\gamma,}^{\textrm{eff,D}}\to g_{a\gamma\gamma}^{\textrm{eff,D}}\approx g_{a\gamma\gamma,\text{const}}^{\textrm{eff,D}}=\begin{cases}\frac{\alpha C_{e}}{\pi\Lambda}\,,&\mbox{for\quad}\text{Scenario E}\,,\\ -\frac{\alpha C_{\mu}}{\pi\Lambda}\frac{m_{a}^{2}}{12m_{\mu}^{2}}\,,&\mbox{for\quad}\text{Scenario M}\,,\\ -\frac{\alpha C_{\tau}}{\pi\Lambda}\frac{m_{a}^{2}}{12m_{\tau}^{2}}\,,&\mbox{for\quad}\text{Scenario T}\,,\\ \frac{\alpha C_{\ell}}{\pi\Lambda}\quantity(1-\frac{m_{a}^{2}}{12m_{\mu}^{2}}-\frac{m_{a}^{2}}{12m_{\tau}^{2}})\,,&\mbox{for\quad}\text{Scenario U}\,.\end{cases}$$ (60)

Similarly, the approximate pseudoscalar couplings read

$$g_{a\gamma^{*}\gamma,}^{\textrm{eff,P}}\to g_{a\gamma\gamma}^{\textrm{eff,P}}\approx g_{a\gamma\gamma,\text{const}}^{\textrm{eff,P}}=\begin{cases}0\,,&\mbox{for\quad}\text{Scenario E}\,,\\ \frac{\alpha g_{a\mu}}{\pi m_{\mu}}\quantity(-1-\frac{m_{a}^{2}}{12m_{\mu}^{2}})\,,&\mbox{for\quad}\text{Scenario M}\,,\\ \frac{\alpha g_{a\tau}}{\pi m_{\tau}}\quantity(-1-\frac{m_{a}^{2}}{12m_{\tau}^{2}})\,,&\mbox{for\quad}\text{Scenario T}\,,\\ \frac{\alpha g_{a\ell}}{\pi m_{\ell}}\quantity(-1-\frac{m_{a}^{2}}{12m_{\mu}^{2}}-1-\frac{m_{a}^{2}}{12m_{\tau}^{2}})\,,&\mbox{for\quad}\text{Scenario U}\,.\end{cases}$$ (61)

Specifically, in the approximation for the case of pseudoscalar couplings to electrons only, we set the loop-induced coupling to photons to zero and consider only the tree-level contribution from electrons. For Scenario U, we also consider a simplified tree-level model of ALPs, where we set all effective loop-induced photon couplings to zero in order to estimate their relevance for the final exclusion bounds.

Figure 9: Bounds derived from E137 using different types of couplings to electrons (top-left), muons (top-right), tau leptons (bottom-left) and universal couplings to all leptons (bottom-right).

Beginning with Scenario E in the top-left panel of figure 9, one can see the bounds on the derivative coupling, the pseudoscalar one and their respective constant approximations. For clarity, only the contributions coming from DB and Primakoff production are shown, but not the positron annihilation processes. For ALP masses $m_{a}\lesssim$100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, the agreement of the bounds between all couplings is close to perfect. However, for larger masses, the differences between the couplings become evident. While the lower bounds on the derivative coupling remain relatively constant until reaching the tip, the pseudoscalar bounds are weaker for the higher ALP mass range.

The derivative and pseudoscalar couplings perfectly match their respective constant approximations. This is unsurprising, especially at the tip of the plot, as the ALP mass is much larger than the electron mass, making eq. (60) an excellent approximation for the effective photon coupling in the derivative basis. Similarly, the fact that the bound for the pseudoscalar coupling can be reproduced by setting the loop-induced photon coupling to zero agrees with the expectation from eq. (18), as the effective photon coupling in the pseudoscalar basis is strongly suppressed by the electron mass.

We note that this excellent agreement between the approximated couplings and the exact couplings is still present at the very left side of the plot, where $m_{a}\approx$10\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$ and the approximation $\tau\ll 1$ starts to become invalid. This agreement can be explained by the fact that the Primakoff production mechanism, which is the contribution sensitive to the effective coupling, is dwarfed by the tree-level DB production channel, as can be seen in figure 5. No matter which coupling type is used, their differing effects are irrelevant compared to DB production, hence the agreement is good for low masses as well. At the tip of the exclusion, on the other hand, the Primakoff production actually outweighs the DB production (see figure 5), such that the two coupling scenarios give slightly different results.

This difference becomes much more pronounced in the top-right panel of figure 9, which shows the comparison of the bounds derived using Scenario M. Since in this scenario only the muon coupling is present at tree level, the only production mechanism in E137 is Primakoff production. Below the threshold $m_{a}=2m_{\mu}\approx$211\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, the ALP can only decay into photons, whereas above the muon threshold, the decay into muons proceeds so fast that the ALPs can never leave the dump and reach the detector. As a result, the bounds on this model are very sensitive to the type of coupling used. Indeed, since the effective coupling in the derivative basis is smaller than in the pseudoscalar basis for large lepton masses, the bounds for the derivative coupling are shifted to higher $\hat{g}_{a\mu}$ to compensate.

In addition to the derivative and pseudoscalar couplings, we also show the bounds derived using constant couplings. For the constant approximation of the pseudoscalar coupling, the bounds are generally in very good agreement with the bounds derived from the full pseudoscalar coupling, except close to the threshold $m_{a}=2m_{\mu}$ where the approximation $\tau\gg 1$ becomes invalid. For the constant approximation of the derivative coupling, the bounds of the full and approximate couplings start to deviate significantly for low masses. In other words, it is essential to consider the full momentum dependence of the effective photon coupling to accurately describe the bounds on the ALP-muon coupling for the derivative coupling case.

One can observe the same effect in Scenario T, which is shown in the bottom-left panel of figure 9. Since the ALP masses are far below the tau threshold, the ALPs can only decay into photons via the effective coupling. As in Scenario M, the ALPs can only be produced via Primakoff production. As a result, the qualitative discussion for Scenarios M and T are almost identical. The only major difference is that, since the relevant ALP mass range remains well below $2m_{\tau}$, the constant coupling approximation works very well for $m_{a}\gtrsim$100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$. For lower masses, however, the momentum dependence of the couplings in the derivative case again becomes apparent.

Finally, the bottom-right panel of figure 9 considers Scenario U, i.e. universal couplings to all leptons. The general shape of the bounds is similar to Scenario E for $m_{a}<2m_{\mu}$, but features the cutoff at $m_{a}=2m_{\mu}$ seen in Scenario M. Since there is now once again a tree-level coupling to the beam electrons, the DB production of ALPs is the dominant process for lower ALP masses, such that pseudoscalar and derivative couplings agree in this mass range. The difference between the two coupling types only becomes apparent at masses around $m_{a}\gtrsim$100\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, when the loop-induced decays into photons start to compete with the tree-level decays into electrons Dolan:2014ska .

Similarly to Scenario E, both the derivative coupling and the pseudoscalar one are described quite accurately by their relative constant coupling approximation in most of the mass range. However, the additional muon and tau couplings in Scenario U lead to a non-trivial momentum dependence in the effective Primakoff coupling near the threshold $m_{a}=2m_{\mu}$, particularly in the pseudoscalar case. By coincidence, the derivative case is well approximated by considering only tree-level couplings, i.e. setting the effective ALP-photon coupling to zero and neglecting both Primakoff production and decays into photons.

5.2 Combined Constraints

Having discussed the pertinent differences between pseudoscalar and derivative couplings for the bounds from E137, as well as the relevance of the momentum dependence in the effective couplings, we now turn to the various other constraints on leptophilic ALPs. Our main results are shown in figures 10 to 13, corresponding to the four scenarios introduced above. In each case, we show in the left (right) panel the constraints for derivative (pseudoscalar) couplings. The bounds that have been (re)computed by us, namely the ones from sections 3 and 4.2, are shown in the foreground, each highlighted with a differently coloured area. Conversely, the constraints from sections 4.1, 4.3, 4.4 and 4.5, which have been taken directly from the literature, are shown in the background as a unified, gray-shaded area. Non-binding constraints, namely the ones coming from anomalous magnetic moments described in section 4.7 and the ALP lifetime discussed in section 4.8, are indicated for reference as dashed lines. Finally, for the NA64$\mu$ experiment, we also indicate projected bounds for a total of $10^{14}$ MOT, compared to the current $1.98\cdot 10^{10}$ MOT.

Figure 10: Bounds on an ALP with derivative (left) and pseudoscalar (right) coupling to electrons.

Figure 11: Bounds on an ALP with derivative (left) and pseudoscalar (right) coupling to muons.

Figure 12: Bounds on an ALP with derivative (left) and pseudoscalar (right) coupling to tau leptons.

Figure 13: Bounds on an ALP with universal derivative (left) and pseudoscalar (right) coupling to all charged leptons.

Let us start our discussion from Scenario E, for which bounds are given in figure 10. Contrary to the previous section, where only contributions from DB and Primakoff production were considered for estimating E137 constraints, an additional bump is now present because we include the resonant ALP production yields as well. We find that E137, which ran over 40 years ago, is still able to set some of the most stringent constraints on the electrophilic ALP models, especially for high ALP masses. However, the bounds set by NA64, both in its visible and invisible modes, complement those by E137, particularly for higher values of the coupling. This is particularly relevant in the pseudoscalar case, where these additional constraints cover all the otherwise unconstrained regions for ALP masses below 10 MeV and couplings above $10^{-2}\,\mathrm{GeV^{-1}}$. While most of the parameter space under consideration is now covered for the derivative case, this is not the case in the pseudoscalar one. In particular, the absence of constraints coming from LEP and $W^{+}$ decays leave the large coupling region for an ALP mass above $\sim$ 200 MeV unconstrained. To determine the viable parameter regions of smaller couplings, it would be desirable to perform a detailed calculation of the constraints from SNe for pseudoscalar couplings, including the effect of loop-induced couplings to photons.

Moving onto Scenario M, for which bounds are shown in figure 11, we observe, as expected, that the strongest constraints come from NA64$\mu$, which is exceptionally well suited to probing muonphilic ALPs. Contrary to the previous scenario, E137 is only relevant for couplings between about $10^{-3}\,\mathrm{GeV}^{-1}$ and $10^{-4}\,\mathrm{GeV}^{-1}$, and masses above 100 (10) MeV, for the derivative (pseudoscalar) case. However, if NA64$\mu$ extends its programme to $10^{14}$ MOT, it is expected to push its constraints down to $10^{-5}\,\mathrm{GeV}^{-1}$, hence totally encompassing E137 ones. Complementary constraints also come from Belle II and from LEP (only in the derivative case) and SNe. We therefore conclude that, as for Scenario E, the parameter space is more constrained in the derivative coupling case than in the pseudoscalar coupling one.

Concerning Scenario T, from figure 12, it becomes clear that currently not many constraints are stronger than the unitarity bounds. We also observe that, in this scenario, the main additional bound comes from E137, which is stronger in the pseudoscalar coupling case than in the derivative one. As a result, in this scenario, the pseudoscalar coupling case is more tightly constrained than the derivative one. Current bounds from NA64 have marginal relevance in the former case, namely only in the ALP mass range $10\text{--}30$ MeV, while in the latter one, they improve the unitarity bound by a factor of 2. The final NA64$\mu$ results will be able to provide much stronger constraints, improving the current bound by about two orders of magnitudes in the coupling and encompassing almost completely the ones from E137 in the derivative case.

We conclude with Scenario U, for which constraints are reported in figure 13. As expected, if the ALP couples to all charged leptons, the different constraints from the three previous scenarios become relevant simultaneously. We therefore find that the combination of E137, NA64$\mu$, BaBar and Belle II measurements exclude ALP couplings above about $10^{-2}\,\mathrm{GeV}^{-1}$. The strongest constraint comes again from E137, which pushes the bound down to almost $10^{-5}\,\mathrm{GeV}^{-1}$ in the $1\text{--}200$ MeV ALP mass range, while for smaller ALP masses NA64 excludes couplings larger than $10^{-3}\,\mathrm{GeV}^{-1}$. With $10^{14}$ MOT, NA64$\mu$ is expected to push the coupling bounds down to about $10^{-5}\,\mathrm{GeV}^{-1}$ for ALP masses up to the muon threshold, going beyond the E137 reach and covering the vast majority of the currently unconstrained parameter space studied here. For the case of universal couplings to charged leptons, there is no relevant difference between the derivative and the pseudoscalar coupling cases, except a small unconstrained window for ALP masses in the $30\text{--}200$ MeV mass range for the latter case.

6 Conclusion

Axion-like particles (ALPs) arising as pseudo-Goldstone bosons from a spontaneously broken approximate global symmetry are one of the most well-motivated ways to extend the Standard Model with a new state below the GeV scale coupled feebly to known particles. While couplings of ALPs to SM gauge bosons and SM quarks have been studied in great detail, the case of leptophilic ALPs, which couple dominantly to charged leptons, has received much less attention. In the present work, we have filled this gap by providing an exhaustive overview of the relevant constraints from beam-dump, flavour and collider experiments, as well as from theoretical considerations and astrophysics.

First of all, we have pointed out that the precise definition of leptophilic ALPs is rather subtle since the assumed absence of tree-level couplings to photons leads to a distinction between ALPs with pseudoscalar couplings and ALPs with derivative couplings. We take a complete model-independent perspective and discuss both interactions in parallel, highlighting their differences.

In both cases, effective interactions of ALPs to photons are generated through lepton loops, which depend on the momenta of the ALP and the two photons. If all particles are on-shell, one reproduces the standard photon coupling, but for off-shell particles, the loop function depends non-trivially on the momenta, which are specific to the relevant experimental setup. We have shown that this effect plays an important role in the analysis of electron beam-dump experiments such as E137 and collider constraints from, e.g. LEP and have derived suitable approximations that can be used to simplify the analysis.

In our analysis, we have then considered separately the three cases where ALPs couple only to one generation of charged leptons, as well as the case of universal couplings to all charged leptons. In all cases, we find that E137 gives highly relevant constraints, even if both the production and decay of the ALPs proceed exclusively via loop-induced photon couplings. Further relevant constraints come in particular from searches for missing energy and displaced decays at NA64, in particular when using a muon beam. Indeed, future runs of NA64$\mu$ constitute a highly promising way to further probe the parameter space of leptophilic ALPs.

For cases in which the ALPs couple to muons, beam dump experiments are not sensitive to ALP masses above the threshold $m_{a}>2m_{\mu}$, where ALPs rapidly decay back into muons. In this parameter region, future searches for di-muon resonances at Belle II Belle-II:2024wtd are particularly promising. Below the threshold, relevant sensitivity may come from the planned proton beam-dump experiment SHiP Albanese:2878604 as well as from the proposed forward-physics experiment FASER2 Adhikary:2024nlv , both of which can exploit the presence of secondary muons in hadron showers to search for leptophilic ALPs. A detailed study of this production mode, including Monte Carlo simulations of the secondary muon distributions, is left for future work. Finally, an improved understanding of the bounds from SN 1987A for all the different coupling scenarios is needed in order to identify the most interesting regions of parameter space at very small couplings. Combining all of these different strategies will enable us to further explore the parameter space of leptophilic ALPs.

Appendix A Experimental details

A.1 E137

An electron beam with an energy of $20\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ was directed at an aluminium target. Over the course of the experiment, a total of $30\text{\,}\mathrm{C}$ (Coulombs) of electrons were deposited in the dump, corresponding to $N_{e}=1.87\cdot 10^{20}$ electrons on target (EOT). The target was placed in front of a hill $179\text{\,}\mathrm{m}$ in length, which absorbed all SM particles except neutrinos. The detector was placed behind the hill, providing a decay length of $204\text{\,}\mathrm{m}$. The detector consisted of an 8-radiation-length shower calorimeter, which could detect photons and electrons. For the first phase of the experiment, the detector was arranged in a $2\times 3$ mosaic of $$1\text{\,}\mathrm{m}$\times$1\text{\,}\mathrm{m}$$ chambers, which was later changed to $$3\text{\,}\mathrm{m}$\times$3\text{\,}\mathrm{m}$$ chambers. The cutoff energy for the detector, i.e. the minimum ALP energy to be considered a candidate, was set to $2\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. No candidate events were observed.

A.2 NA64

Invisible mode

The detector consists of one electromagnetic calorimeter (ECAL) and three hadronic calorimeters (HCAL). The ECAL consists of lead and scintillator plates, having a total length of $45\text{\,}\mathrm{c}\mathrm{m}$ and a width of $$23\text{\,}\mathrm{c}\mathrm{m}$\times$23\text{\,}\mathrm{c}\mathrm{m}$$. The ECAL is the active target and detects the energy deposition of electrons. Each HCAL is made of iron and scintillator plates intended to detect charged and neutral hadrons. The three HCALs have a combined length of approximately $5\text{\,}\mathrm{m}$, and a width of $$60\text{\,}\mathrm{c}\mathrm{m}$\times$60\text{\,}\mathrm{c}\mathrm{m}$$. Using this setup, the calorimeters are able to determine the energy deposited inside the beam dump. If an ALP or other weakly interacting particle were to be produced inside the dump and was long-lived enough, it would escape the detector and carry away a fraction of the original energy of the beam. This missing energy can be measured and registered as an event Gninenko:2719646 .

The main beam has an energy of $100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. The energy cutoff for the detector is set at $50\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. In the runs from 2016–2022, NA64 had a total of $N_{e}=9.37\cdot 10^{11}$ EOT. In the 2016 run, a fourth HCAL was placed in line with the other three HCALs, which was moved downstream to be a zero-degree HCAL to help with background detection in later runs. A detailed description of the layout of the experiment can be found in ref. MolinaBueno:2868332 . No candidate events were observed NA64:2023wbi .

Visible mode

The setup used for the visible mode is described in ref. NA64:2019auh . The beam is directed toward an electromagnetic calorimeter composed of tungsten and scintillator plates (WCAL). The WCAL serves as an active beam dump, with the first layers being the main regions in which BSM particles would be produced. The rest of the WCAL serves as a dump for the SM particles created in other processes, while weakly interacting particles would pass through the dump. After a decay length of approximately $3.5\text{\,}\mathrm{m}$, the WCAL is followed by the ECAL and HCALs discussed above. If an ALP or any other particle was produced and subsequently decayed inside the decay volume, the ECAL would register the energy of the decay products. If the energy that is deposited in the ECAL is the same as the missing energy in the WCAL, i.e. $E_{\textrm{ECAL}}=E_{\textrm{beam}}-E_{\textrm{WCAL}}$, then this is a candidate event Gninenko:2320630 .

In 2017, the beam was run at an energy of $100\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. The length of the WCAL was varied from $290\text{\,}\mathrm{m}\mathrm{m}$, with $2.4\cdot 10^{10}$ EOT, to $220\text{\,}\mathrm{m}\mathrm{m}$, with $3\cdot 10^{10}$ EOT. The cutoff energy for the detector was set to $20\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. In 2018, the beam energy was increased to $150\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$, with $3\cdot 10^{10}$ EOT, with the cutoff energy being set to $30\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. No candidate events were observed NA64:2019auh ; Gninenko:2320630 .

Muon mode

The experimental setup for NA64$\mu$ is illustrated in ref. NA64:2024klw . As with the original NA64 run with an electron beam, the experiment is built to detect the invisible decay of BSM particles by measuring the missing energy of an event. A muon beam with an energy of $160\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ is directed at an ECAL, similar to the one in previous runs. The ECAL has a total length of $49\text{\,}\mathrm{c}\mathrm{m}$. Further downstream, a magnetic spectrometer is placed, which uses a bending magnet (MS2) to deflect the muon beam, as well as two HCALs to remove residuals from the upstream detectors NA64:2024klw . The muon momentum can then be reconstructed. Candidate events are defined as events in which the reconstructed muon momentum is smaller than $80\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$ and the energy deposition in the calorimeters is compatible with a minimum ionising particle, such that the missing energy is larger than $80\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$. This setup results in a shielding of approximately $11\text{\,}\mathrm{m}$, which the ALP must traverse before decaying. After the 2021 test run and the 2022-2023 runs, the experiment accumulated a total of $1.98\cdot 10^{10}$ muons on target (MOT). No candidate events were observed Crivelli:2907892 .

A.3 Summary of Experimental Parameters

For ease of reference, the values of all experimental parameters we use to calculate $N_{a}$ are listed in Tables 1 and 2.

Experiment	Target	Z	A	$X$	$\rho$	$M_{\textrm{target}}$
				$\quantity[$\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-2}$]$	$\quantity[$\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$]$	$[$\text{\,}\mathrm{g}\text{\,}{\mathrm{mol}}^{-1}$]$
E137	Al	$13$	$27$	$24.01$	$2.7$	$26.98$
NA64 (invisible,$\mu$)	Pb	$82$	$207$	$6.37$	$11.35$	$207.2$
NA64 (visible)	W	$74$	$184$	$6.76$	$19.3$	$183.84$

Table 1: Material parameters for E137 and NA64 ParticleDataGroup:2024cfk .

Experiment	EOT/MOT	$E_{0}$	$E_{\textrm{cut}}$	$L_{\textrm{sh}}$	$L_{\textrm{dec}}$	$T$	$\theta_{a,\text{max}}$
		[GeV]	[GeV]	[m]	[m]
E137	$1.87\cdot 10^{20}$	$20$	$2$	$179$	$204$	$2012.9$	$4.4\cdot 10^{-3}$
NA64 (invisible)	$9.37\cdot 10^{11}$	$100$	$50$	$5.5$	$\infty$	$979.8$	$0.11$
NA64 (visible):
Run 1 (2017)	$2.4\cdot 10^{10}$	$100$	$20$	$0.29$	$3.5$	$82.8$	$6.9\cdot 10^{-2}$
Run 2 (2017)	$3\cdot 10^{10}$	$100$	$20$	$0.22$	$3.5$	$62.8$	$6.9\cdot 10^{-2}$
Run 3 (2018)	$3\cdot 10^{10}$	$150$	$30$	$0.22$	$3.5$	$62.8$	$6.9\cdot 10^{-2}$
NA64$\mu$	$1.98\cdot 10^{10}$	$160$	$80$	11	$\infty$	$87.3$	$5.5\cdot 10^{-2}$

Table 2: Experimental parameters for E137 and NA64.

Since the experimental parameters were varied for the different runs of NA64 operated in visible mode, the sum over all runs is taken when determining the total number $N_{a}$

$$N_{a}=\sum_{\text{runs}}N_{a,i}\,.$$ (62)

The invisible mode of the NA64 experiment has three HCALs downstream of the ECAL, giving a total length of approximately $5.5\text{\,}\mathrm{m}$. In the first run in 2016, the experiment had a fourth HCAL, which was later moved to a zero-degree HCAL. This means the length of the shielding of the 2016 run was greater compared to the later runs. Nevertheless, we approximate $L_{\textrm{sh}}=$5.5\text{\,}\mathrm{m}$$ for all runs, as the error would be minimal on account of the relatively low contribution of the $4.5\cdot 10^{10}$ EOT from this run to the total $9.37\cdot 10^{11}$ EOT.

Appendix B Details on cross section calculations

B.1 Dark Bremsstrahlung and Primakoff Production

It is convenient to define modified Mandelstam variables $\tilde{s}$, $\tilde{u}$, and $t_{2}$ in the same way as Liu:2017htz , noting the use of the opposite sign of the metric. Assuming the incoming lepton has an energy that is much greater than the masses of the ALP and lepton, one can make the approximations Liu:2017htz

$$\tilde{s}=-\frac{\tilde{u}}{1-x}\,,\quad\tilde{u}=-xE^{2}\theta_{a}^{2}-m_{a}^{2}\frac{1-x}{x}-m_{\ell}^{2}x\,,\quad t_{2}=\frac{x}{1-x}\tilde{u}+m_{a}^{2}\,.$$ (63)

In the IWW approximation, the cross section is given by Liu:2017htz

$$\quantity(\derivative{\sigma}{x})_{\textrm{IWW}}=\frac{\alpha}{(4\pi)^{2}}\sqrt{x^{2}-\frac{m_{a}^{2}}{E^{2}}}\chi\frac{1-x}{x}\int_{\tilde{u}_{\textrm{min}}}^{\tilde{u}_{\textrm{max}}}\differential\tilde{u}\frac{\overline{\left|\mathcal{M}\right|}^{2}}{\tilde{u}^{2}}\,,$$ (64)

with the boundaries

$$\tilde{u}_{\textrm{min}}=-xE^{2}\theta_{a,\text{max}}^{2}-m_{a}^{2}\frac{1-x}{x}-m_{\ell}^{2}x,\quad\tilde{u}_{\textrm{max}}=-m_{a}^{2}\frac{1-x}{x}-m_{\ell}^{2}x\,,$$ (65)

where the maximum scattering angle $\theta_{a,\text{max}}$ depends on the geometry of the detector. The effective flux of photons $\chi$ is given by

$$\chi=\int_{t_{\textrm{min}}}^{t_{\textrm{max}}}\differential t\frac{t-t_{\textrm{min}}}{t^{2}}F^{2}(t)\,.$$ (66)

Here, $F$ is the form factor. The inelastic form factor can be neglected here as its contribution is negligible compared to the elastic form factor Liu:2017htz , given by

$$F(t)=\quantity(\frac{a^{2}t}{1+a^{2}t})\quantity(\frac{1}{1+t/d})Z\,.$$ (67)

The first part accounts for the atomic form factor, while the second part accounts for the nuclear form factor Liu:2023bby ; Fayet:2024ddk . Here, $a=111Z^{1/3}/m_{e}$ and $d=$0.164\text{\,}\mathrm{{}^{2}}{GeV}$A^{-2/3}$. The integers $Z$ and $A$ are the atomic number and mass number of the target atom, respectively.

Using FeynCalc Shtabovenko:2023idz , we calculate that the spin-averaged, squared amplitude for the simplified DB processes is

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{DB}}}}^{2}=\frac{e^{2}g_{a\ell}^{2}\Bigl{[}2m_{a}^{2}(\tilde{s}+\tilde{u})\left(m_{\ell}^{2}(\tilde{s}+\tilde{u})+\tilde{s}\tilde{u}\right)-\tilde{s}\tilde{u}(\tilde{s}+\tilde{u})^{2}-2m_{a}^{4}\tilde{s}\tilde{u}\Bigr{]}}{\tilde{s}^{2}\tilde{u}^{2}}\,.$$ (68)

Using eq. (64) and the substitutions in eq. (63), we then obtain eq. (37). Since the DB diagrams have tree-level couplings to the ALP, the choice between a pseudoscalar or derivative ALP-lepton coupling has no effect on the amplitude.

For the Primakoff process, the amplitude is

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{P}}}}^{2}=\frac{e^{2}\absolutevalue{g_{a\gamma^{*}\gamma}^{\textrm{eff}}(t_{2})}^{2}\Bigl{[}-t_{2}\left(\tilde{s}^{2}+\tilde{u}^{2}\right)-2m_{\ell}^{2}(\tilde{s}+\tilde{u})^{2}\Bigr{]}}{4t_{2}^{2}}\,,$$ (69)

where the effective off-shell coupling $g_{a\gamma^{*}\gamma}^{\textrm{eff}}(t_{2})$ is defined in eq. (25) and should be understood as the sum over the three effective couplings coming from the couplings to the leptons, i.e.

$$g_{a\gamma^{*}\gamma}^{\textrm{eff}}(t_{2})=\sum_{\ell=e,\mu,\tau}g_{a\gamma^{*}\gamma}^{\textrm{eff}}(t_{2},m_{\ell})\,.$$ (70)

Additionally, the interference term between the two processes gives the contribution

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{I}}}}^{2}=\frac{e^{2}\real\quantity(g_{a\gamma^{*}\gamma}^{\textrm{eff}}(t_{2}))g_{a\ell}m_{\ell}(\tilde{s}+\tilde{u})^{3}}{\tilde{s}t_{2}\tilde{u}}\,.$$ (71)

B.2 Non-Resonant Production

The Feynman diagrams for the non-resonant production of an ALP

$$e^{-}(p_{1})+e^{+}(p_{2})\rightarrow a(k_{1})+\gamma(k_{2})$$ (72)

are shown in figure 14.

Figure 14: Feynman diagrams contributing to the process $e^{-}e^{+}\rightarrow a\gamma$.

The positron can annihilate with the electron through either an s-channel, t-channel or u-channel diagram, of which the latter two are tree-level diagrams. The squared, spin-averaged matrix elements for the tree-level, loop-level, and interference terms are

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{NR,tree}}}}^{2}= \frac{e^{2}g_{ae}^{2}}{\left(t-m_{e}^{2}\right)^{2}\left(u-m_{e}^{2}\right)^{2}}$$

$$\times\Bigl{[}\left(t-m_{e}^{2}\right)\left(u-m_{e}^{2}\right)\left(2m_{a}^{4}+(s-m_{a}^{2})^{2}\right)+2m_{a}^{2}\left(s-m_{a}^{2}\right)\left(tu-m_{e}^{4}\right)\Bigr{]}\,,$$

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{NR,loop}}}}^{2}= \frac{e^{2}\absolutevalue{g_{a\gamma^{*}\gamma}^{\textrm{eff}}(s)}^{2}}{4s^{2}}$$ (73)

$$\times\Bigl{[}m_{a}^{4}\left(2m_{e}^{2}+s\right)-2m_{a}^{2}s\left(m_{e}^{2}+s+t\right)+2s\left(m_{e}^{4}-2m_{e}^{2}t+t(s+t)\right)+s^{3}\Bigr{]}\,,$$

$$\overline{\absolutevalue{\mathcal{M}_{\textrm{NR,inter}}}}^{2}= \frac{e^{2}\real\quantity(g_{a\gamma^{*}\gamma}^{\textrm{eff}}(s))g_{ae}m_{e}\left(s-m_{a}^{2}\right)^{3}}{s\left(t-m_{e}^{2}\right)\left(u-m_{e}^{2}\right)}\,,$$

respectively. The Mandelstam variables are defined in the standard way

$$s=(p_{1}+p_{2})^{2},\quad t=(p_{1}-k_{1})^{2},\quad u=(p_{1}-k_{2})^{2}\,.$$ (74)

The cross section for non-resonant production of ALPs is

$$\derivative{\sigma_{\textrm{NR}}}{x} =\frac{1}{16\pi s}\frac{1}{1-\frac{m_{e}}{E}}\overline{\absolutevalue{\mathcal{M}_{\textrm{NR}}}}^{2}\theta(\theta_{a,\text{max}}-\theta_{a})\,,$$ (75)

with the angle $\theta_{a}$ defined via

$$\cos(\theta_{a})\equiv\frac{2(E+m_{e})E_{a}-s-m_{a}^{2}}{2\absolutevalue{\mathbf{p}}\absolutevalue{\mathbf{k_{a}}}}=\frac{2(E+m_{e})(E_{a}-m_{e})-m_{a}^{2}}{2\absolutevalue{\mathbf{p}}\absolutevalue{\mathbf{k_{a}}}}\,.$$ (76)

The final Heaviside step function ensures that the trajectory of the ALP is within the angular coverage of the detector. For non-resonant annihilation, the integration boundaries in eqs. (30) and (31) are not applicable. Instead, using the Mandelstam variable $s=2m_{e}(E+m_{e})$, the integration boundaries are

$$x_{\textrm{min}} =\frac{1}{E}\cdot\max\quantity[E_{\textrm{cut}}\,,\,\frac{1}{4m_{e}}\quantity(s+m_{a}^{2}-\sqrt{1-\frac{4m_{e}^{2}}{s}}(s-m_{a}^{2}))]\,,$$ (77)

$$x_{\textrm{max}} =\frac{1}{E}\cdot\frac{1}{4m_{e}}\quantity(s+m_{a}^{2}+\sqrt{1-\frac{4m_{e}^{2}}{s}}(s-m_{a}^{2}))\,,$$ (78)

$$E_{\textrm{min}} =\frac{s_{\textrm{min}}}{2m_{e}}-m_{e}\,,\quad E_{\textrm{max}}=E_{0}\,,$$ (79)

$$s_{\textrm{min}} =\max\quantity[m_{a}^{2}\,,\,m_{e}\frac{(m_{a}^{2}-2m_{e}E_{\textrm{cut}})\quantity(\sqrt{E_{\textrm{cut}}^{2}-m_{a}^{2}}+E_{\textrm{cut}})+m_{a}^{2}m_{e}}{m_{a}^{2}-2m_{e}E_{\textrm{cut}}+m_{e}^{2}}]\,.$$ (80)

It is important to note that for $m_{a}>\sqrt{2m_{e}E_{\textrm{cut}}}$, the diagrams contributing to the non-resonant production of ALPs exhibit an infrared divergence at the lower integration bound $s_{\textrm{min}}=m_{a}^{2}$ due to the low-energy photon. This could be remedied by including one-loop corrections to the ALP-lepton vertex, which are, however, beyond the scope of the present work. We therefore neglect the sub-leading non-resonant production in this region of parameter space and focus solely on the dominant contribution from resonant production.

Appendix C Comparison of E137 with other implementations

Figure 15: Constraints on the tree-level electron coupling from E137. The result from this work is shown as the blue-shaded region, while the results from refs. Liu:2023bby and Altmannshofer:2022ckw are shown as solid red and dashed green lines, respectively.

In this appendix, we compare our implementation of E137 and the resulting constraints on the ALP-lepton couplings with previous studies in the literature. In ref. Liu:2023bby , the concurrent effects of an ALP-photon and ALP-electron coupling were analysed. The authors give the constraints on a derivative coupling $g_{a\overline{e}e}^{\textrm{eff}}$ for various values of $g_{a\gamma\gamma}^{\textrm{eff}}$, which are however taken to be constant. Our results can be compared by taking the constraints in ref. Liu:2023bby for $g_{a\gamma\gamma}^{\textrm{eff}}=0$ and comparing them to our tree-level model constraints from figure 9 (“Constant (P)” for Scenario E). This comparison is shown in figure 15. For low ALP masses, our constraints are in good agreement with those in ref. Liu:2023bby , but we obtain substantially stronger bounds for high masses.

It should be noted that ref. Liu:2023bby does not consider the contributions to ALP production arising from electron-positron annihilation. Therefore, we also compare our constraints with those given in ref. Altmannshofer:2022ckw , based on ref. Bjorken:1988as , where this contribution is included. The original analysis from ref. Bjorken:1988as did not discuss the bounds for large couplings, which were estimated in ref. Altmannshofer:2022ckw assuming the ALPs to be monoenergetic. Clearly, this approximation is not very accurate and leads to a very different upper bound. However, the shape of the bound coming from the annihilation of electrons and positrons matches that of ours relatively well.

The constant shift between the lower bound from ref. Bjorken:1988as and ours may be attributed to the different implementations of the track length. Indeed, the track length distribution that we use (see figure 3) is larger for low energies than the one given in figure 14 of ref. Bjorken:1988as , which is based on the (unpublished) EGS shower code. We note that our implementation agrees very well with the results from ref. Marsicano:2018krp , which has run GEANT4 simulations for E137. Another potential difference between our analysis and the one from ref. Bjorken:1988as is that the E137 analysis uses a much more sophisticated implementation of the detector acceptance, tracking the ALP decay products and their interactions with the detector, which is beyond the scope of our analysis.

As discussed in appendix B.2, our bounds have a discontinuity at $m_{a}=\sqrt{2m_{e}E_{\textrm{cut}}}\approx$45.2\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$$, because we conservatively neglect the non-resonant production channels. Including these channels is expected to further improve the agreement with ref. Bjorken:1988as , which finds a continuous bound.