aainstitutetext: Instituto de Física Teórica UAM-CSIC, Calle de Nicolás
Cabrera 13–15, Universidad Autónoma de Madrid, Cantoblanco, E-28049 Madrid,
Spainbbinstitutetext: Instituto de Física Corpuscular, Universidad de Valencia
& CSIC, Edificio Institutos de Investigación, Calle Catedrático José Beltrán
2, 46980 Paterna, Spain

# New bounds on axion-like particles from MicroBooNE

Pilar Coloma<EMAIL_ADDRESS>b Pilar Hernández
<EMAIL_ADDRESS>b Salvador Urrea<EMAIL_ADDRESS>
###### Abstract

Neutrino experiments lie at the edge of the intensity frontier and therefore
can be exploited to search for new light particles weakly coupled to the
visible sector. In this work we derive new constraints on axion-like particles
(ALPs) using data from the MicroBooNE experiment, from a search for
$e^{+}e^{-}$ pairs pointing in the direction of the NuMI absorber. In
particular, we consider the addition of higher-dimensional effective operators
coupling the ALP to the electroweak gauge bosons. These would induce $K\to\pi
a$ from kaon decay at rest in the NuMI absorber, as well as ALP decays into
pairs of leptons or photons. We discuss in detail and compare various results
obtained for the decay width $K\rightarrow\pi a$ in previous literature. For
the operator involving the Higgs, MicroBooNE already sets competitive bounds
(comparable to those of NA62) for ALP masses between 100 and 200 MeV. We also
compute the expected sensitivities from the full NuMI dataset recorded at
MicroBooNE. Our results show that a search for a $a\to\gamma\gamma$ signal may
be able to improve over current constraints from beam-dump experiments on the
operator involving the ALP coupling to the $W$.

###### Keywords:

axion-like particles, neutrino experiments

††preprint: FTUV-22-0121.2990, IFT-UAM/CSIC-22-4, IFIC/22-03

##  1 Introduction

Light pseudoscalars may arise as pseudo-Nambu-Goldstone bosons of a
spontaneous global symmetry breaking, and are therefore ubiquitous in
extensions of the Standard Model (SM). These particles as often referred to as
axion-like particles (ALPs) since the best motivated example is the axion
Weinberg:1977ma ; Wilczek:1977pj , in models that solve the strong CP problem
and could explain dark matter. If the scale of symmetry breaking, $f_{a}$, is
much higher than the electroweak (EW) scale, these particles might be the only
accessible BSM relics of the high-energy theory.

In recent years, model-independent constraints on these particles have been
systematically studied using an Effective Field Theory (EFT) framework
Georgi:1986df ; Izaguirre:2016dfi ; Gavela:2019wzg ; Bauer:2020jbp ;
Chala:2020wvs ; Bauer:2021wjo ; Bauer:2021mvw ; Bonilla:2021ufe ;
Guerrera:2021yss , that includes all the possible couplings of the ALPs
allowed by the SM symmetries and the shift symmetry. At the renormalizable
level, only a soft-breaking ALP mass term is included, while the leading non-
trivial interactions are of dimension five, and therefore suppressed as
$f^{-1}_{a}$. As a result, ALPs fall in the category of _feebly interacting
particles_. Although the Peccei-Quinn solutions to the strong CP problem
Peccei:1977hh  imply a very restrictive relation between the axion mass111The
physical mass depends on the tree-level mass arising in some UV extensions,
and generically on the non-perturbative QCD contribution to the axion
potential.  and $f_{a}$, fixed by QCD dynamics, alternative scenarios have
been considered where this relation may be significantly modified, increasing
the target parameter space accordingly (see for example Ref. DiLuzio:2020wdo
for a recent review).

The experimental search for axions and ALPs has intensified enormously in
recent years. This has been driven in part by the operation of new axion
search experiments, but also by the identification of new and complementary
opportunities offered by experiments in the intensity frontier at accelerator
facilities, both in beam-dumps and colliders. Stringent limits have been
recently obtained from $B$-meson decays Belle-II:2020jti ; LHCb:2015nkv ;
LHCb:2016awg , from kaon decays BNL-E949:2009dza ; NA62:2021zjw  and from
neutrino experiments Essig:2010gu . Very important constraints from past beam-
dump experiments CHARM:1983ayi ; Bjorken:1988as  have also been recasted in
the context of more general ALP searches Dolan:2017osp . Furthermore many
future experiments are being proposed that can improve the limits further (for
a recent review see Ref. Agrawal:2021dbo ).

In this paper we derive a new constraint on ALPs from recent data by the
MicroBooNE experiment MicroBooNE:2021usw . In neutrino experiments using a
conventional beam (such as MicroBooNE), most of the mesons produced in the
proton collisions in the target are focused by the magnetic horns and decay
within the decay pipe. However, a fraction of the protons (which can be as
large as 10-15%) does not get stopped by the target and ends up hitting the
absorber at the end of the decay pipe. The kaons produced in such collisions,
after losing energy as they interact with the medium, eventually decay at rest
within the absorber. Although the fraction of protons which reach the absorber
is much smaller than the ones hitting the target, the detector acceptance will
be much larger simply because of its proximity. This fact was precisely used
for the analysis performed in Ref. MicroBooNE:2021usw , where the MicroBooNE
collaboration searched for monoenergetic scalars coming from the direction of
the NuMI hadron absorber, located at a distance of only 100 m from the
detector.

In this work, we concentrate on the ALP couplings to EW gauge bosons and on
the most relevant ALP production channel, $K\rightarrow\pi a$. The theoretical
prediction of this process has been recently revisited in Refs. Bauer:2021wjo
; Bauer:2021mvw , including running effects and treating the long-distance
dynamics in the context of the chiral effective theory, as originally proposed
in Georgi:1986df . Alternative estimates for this channel have been used
before in Refs. Izaguirre:2016dfi ; Gavela:2019wzg ; Guerrera:2021yss . We
compare the two approaches in detail, and use the results of Ref.
Bauer:2021wjo  to derive our bounds. We then provide a recasting of the
results of the MicroBooNE analysis in Ref. MicroBooNE:2021usw  and re-derive a
competitive bound on the effective EW couplings of ALPs, using the same decay
topology as in their analysis.

The paper is organized as follows. In Sec. 2 we review the ALP EFT as well as
the computation of $K\rightarrow\pi a$ when EW couplings are present. Section
3 discusses the most relevant bounds on this scenario. Our results are then
presented in Sec. 4: first, we derive new bounds using published MicroBooNE
data, comparing them to current constraints; next, we compute the projected
sensitivity using the full statistics expected, and study the interference
effects between different effective operators. Our summary and conclusions are
then presented in Sec. 5. Appendices A and B contain technical details
relevant for the discussion in Sec. 2.

##  2 Electroweak Lagrangian with effective ALP interactions

We consider the following set of effective operators describing ALP
interactions with EW bosons at some high scale, $\Lambda\propto f_{a}$

$\delta\mathcal{L}_{\rm
EW}=c_{\phi}\mathcal{O}_{\phi}+c_{B}\mathcal{O}_{B}+c_{W}\mathcal{O}_{W}\,,$
(1)

where $c_{i}$ stand for the Wilson coefficients of each operator:

$\displaystyle\mathcal{O}_{\phi}$
$\displaystyle=i\frac{\partial^{\mu}a}{f_{a}}\phi^{\dagger}\overleftrightarrow{D}_{\mu}\phi\,,$
$\displaystyle\mathcal{O}_{B}$
$\displaystyle=-\frac{a}{f_{a}}B_{\mu\nu}\widetilde{B}_{\mu\nu}\,,$ (2)
$\displaystyle\mathcal{O}_{W}$
$\displaystyle=-\frac{a}{f_{a}}W^{a}_{\mu\nu}\widetilde{W}^{a}_{\mu\nu}\,.$

Here, $\phi$ is the Higgs doublet while $B$ and $W^{a}$ stand for the EW
vector bosons, and $a$ is the ALP field. The dual field strengths are defined
as
$\widetilde{X}^{\mu\nu}\equiv\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}X_{\rho\sigma}$,
with $\epsilon^{0123}=1$, and
$\phi^{\dagger}\overleftrightarrow{D}_{\mu}\phi\equiv\phi^{\dagger}\big{(}D_{\mu}\phi\big{)}-\big{(}D_{\mu}\phi\big{)}^{\dagger}\phi{}$.

It can be shown that the operator $\mathcal{O}_{\phi}$ can be eliminated
through a hypercharge rotation Georgi:1986df

$\begin{array}[]{rcl}\phi&\to&e^{ic_{\phi}\frac{a}{f_{a}}}\phi\,,\\\
\Psi_{F}&\to&e^{2iY_{F}c_{\phi}\frac{a}{f_{a}}}\Psi_{F}\,,\end{array}$ (3)

where $F$ runs over left- and right-handed quarks and charged leptons (with
hypercharge $Y_{F}$). This trades the $\mathcal{O}_{\phi}$ operator by a set
of fermionic operators of the form:

${\partial_{\mu}a(x)\over
f_{a}}\sum_{F}~{}\bar{\Psi}_{F}\gamma^{\mu}\Psi_{F}\,.$ (4)

###  2.1 ALP production in kaon decays

While the $\mathcal{O}_{\phi}$ and $\mathcal{O}_{W}$ operators are flavour-
blind, they can lead to flavor-changing neutral-current (FCNC) processes via
the exchange of gauge bosons at one loop, as discussed in detail in Refs.
Izaguirre:2016dfi ; Gavela:2019wzg . This will induce new kaon decay modes
(_e.g.,_ $K\to\pi a$).

We are interested in these hadronic processes at energies below the EW scale.
The relevant couplings are the induced ALP couplings to the light quark
currents in Eq. (4)

${\partial_{\mu}a(x)\over
f_{a}}\left(\sum_{q}~{}\bar{q}_{R}~{}k_{q}\gamma^{\mu}q_{R}+\sum_{Q}\bar{Q}_{L}~{}k_{Q}\gamma^{\mu}Q_{L}\right)\,,$
(5)

where following the notation of Ref. Bauer:2021wjo , lower and upper case for
the quark field refers to the right-handed/left-handed quarks respectively,
and $k_{q}$ and $k_{Q}$ are $3\times 3$ matrices with indices $(u,d,s)$. The
flavour-diagonal axial-current couplings are222Here we neglect the differences
between the couplings in the flavor and mass bases as they will only introduce
subleading corrections in our case (we refer the interested reader to Ref.
Bauer:2020jbp  for a detailed discussion).

$\displaystyle{\partial_{\mu}a(x)\over
2f_{a}}\sum_{q=u,d,s}c_{qq}~{}\bar{q}~{}\gamma^{\mu}\gamma^{5}q\,,\quad\mathrm{with}\quad
c_{qq}\equiv[k_{q}-k_{Q}]_{qq}\,.$ (6)

The standard approach to incorporate these non-standard interactions in
hadronic physics is to match the theory to Chiral Perturbation Theory
($\chi$PT) as first proposed in Ref. Georgi:1986df . Recently, this approach
has been revisited and some previous inconsistencies related to weak hadron
decays have been corrected Bauer:2021wjo , in particular for the important
decay $K\to\pi a$. The result for the amplitude found in Ref. Bauer:2021wjo
is

$\displaystyle i\mathcal{A}_{K\to\pi a}=-$
$\displaystyle\frac{m_{K}^{2}-m_{\pi}^{2}}{2f_{a}}[k_{q}+k_{Q}]_{ds}$
$\displaystyle+\frac{N_{8}}{4f_{a}}$
$\displaystyle\left\\{6\left[c_{uu}+c_{dd}-2c_{ss}\right]\frac{m_{a}^{2}(m_{K}^{2}-m_{a}^{2})}{4m_{K}^{2}-m_{\pi}^{2}-3m_{a}^{2}}\right.$
$\displaystyle+\left[2c_{uu}+c_{dd}+c_{ss}\right](m_{K}^{2}-m_{\pi}^{2}-m_{a}^{2})+4c_{ss}m_{a}^{2}$
$\displaystyle\left.+\left([k_{q}+k_{Q}]_{dd}-[k_{q}+k_{Q}]_{ss}\right)(m_{K}^{2}+m_{\pi}^{2}-m_{a}^{2})\right\\}\,,$
(7)

where $m_{a},m_{\pi},m_{K}$ stand for the ALP, pion and kaon masses,
respectively, and all the couplings in this amplitude are assumed at an energy
scale $\mu=2$ GeV. The first term results from the induced flavor-changing
vector currents in the EFT at low energies, while the rest of the terms are
proportional to $N_{8}\simeq-
g_{8}G_{F}/\sqrt{2}V_{us}V_{ud}^{*}f_{\pi}^{2}\sim\mathcal{O}(1.5\times
10^{-7})$, which is the dominant coupling of the weak Hamiltonian in $\chi$PT
that mediates the standard ($\Delta I=1/2$) $K\rightarrow\pi\pi$ decays (see
Refs. Bauer:2021wjo ; Cirigliano:2011ny  for details).

Next we need to connect the couplings $k_{q}(\mu),k_{Q}(\mu),c_{qq}(\mu)$ to
the high-energy parameters, $c_{\phi}(\Lambda),c_{W}(\Lambda)$ and
$c_{B}(\Lambda)$ in Eq. (1). At tree level the only contribution comes from
the operator ${\mathcal{O}}_{\phi}$, leading to flavour-diagonal and universal
couplings of the form:

$k_{u}(\Lambda)=-{4\over 3}c_{\phi}(\Lambda),\quad k_{d}(\Lambda)={2\over
3}c_{\phi}(\Lambda),\quad k_{Q}(\Lambda)=-{1\over 3}c_{\phi}(\Lambda),$ (8)

for quarks and

$k_{e}(\Lambda)=2c_{\phi}(\Lambda),\quad k_{L}(\Lambda)=c_{\phi}(\Lambda),$
(9)

for charged leptons. This implies

$\displaystyle
c_{ff}(\Lambda)=\left\\{\begin{array}[]{ll}-c_{\phi}(\Lambda)&\mathrm{for}~{}f=u,c,t\,,\\\
c_{\phi}(\Lambda)&\mathrm{for}~{}f=d,s,b\,,\\\
c_{\phi}(\Lambda)&\mathrm{for}~{}f=e,\mu,\tau\,.\end{array}\right.$ (13)

The axion operators mix under renormalization and, at the weak scale
$\mu_{w}$, flavour non-diagonal contributions result from the matching to the
low-energy EFT (_i.e.,_ without $W,Z,t$). Let us stress that that the
normalization used for the operators in Ref. Bauer:2021wjo  differs from the
one in Eq. (1). Denoting the former as $\tilde{c}_{i}$, they are related as

$\displaystyle{\alpha_{2}(\Lambda)\over 4\pi}\tilde{c}_{W}(\Lambda)\equiv-
c_{W}(\Lambda),~{}~{}{\alpha_{1}(\Lambda)\over
4\pi}\tilde{c}_{B}(\Lambda)\equiv-c_{B}(\Lambda)\,,$ (14)

with $\alpha_{1}\equiv\alpha/c_{w}^{2}$, $\alpha_{2}\equiv\alpha/s_{w}^{2}$,
where $s_{w}^{2}\equiv\sin^{2}\theta_{w}\simeq 0.231$ is the Weinberg angle
and $c_{w}^{2}=1-s_{w}^{2}$. Setting $\mu_{w}$ to the top mass ($m_{t}$), the
effective coupling $[k_{Q}]_{ds}$ reads Bauer:2021wjo

$[k_{Q}(m_{t})]_{ds}=V_{td}^{*}V_{ts}\left\\{-\frac{1}{6}I_{t}(m_{t},\Lambda)+\frac{y_{t}^{2}(m_{t})}{16\pi^{2}}c_{tt}(m_{t})f(x_{t})-\frac{y_{t}^{2}(m_{t})}{16\pi^{2}}\frac{3\alpha(m_{t})}{2\pi
s_{w}^{2}}\tilde{c}_{W}(\Lambda)h(x_{t})\right\\}\,,$ (15)

where

$h(x)\equiv\frac{1-x+x\log
x}{(1-x)^{2}},~{}f(x)\equiv-\frac{1}{2}\left[\frac{1}{2}+3h(x)\right]\,,$ (16)

while $V_{qq^{\prime}}$ stands for the CKM matrix elements, $y_{t}$ is the
top-quark Yukawa, $x_{t}\equiv m_{t}^{2}/m_{W}^{2}$ ($m_{W}$ being the mass of
the $W$), and $\alpha$ is the fine structure constant. The function $I_{t}$
and the top axial current $c_{tt}$ depend on the couplings at $\Lambda$ as
follows:

$\displaystyle I_{t}(\mu_{w},\Lambda)$ $\displaystyle=$
$\displaystyle-\tilde{c}_{W}(\Lambda)I_{1}(\mu_{w},\Lambda)-\tilde{c}_{B}(\Lambda)I_{3}(\mu_{w},\Lambda)$
(17) $\displaystyle-$ $\displaystyle
c_{\phi}(\Lambda)\left[I_{5}(\mu_{w},\Lambda)-\frac{2}{3}A(\mu_{w},\Lambda)\right],$

and

$\displaystyle c_{tt}(\mu_{w})$ $\displaystyle=$
$\displaystyle\tilde{c}_{W}(\Lambda)I_{2}(\mu_{w},\Lambda)+\tilde{c}_{B}(\Lambda)I_{4}(\mu_{w},\Lambda)$
(18) $\displaystyle-$ $\displaystyle
c_{\phi}(\Lambda)\left[1-A(\mu_{w},\Lambda)-I_{6}(\mu_{w},\Lambda)\right]\,,$

where the integrals $I_{1-6}$ and the function $A$ are given in App. A.

Substituting Eqs. (17) and (18) into Eq. (15) we obtain the following
expression in terms of the high-energy parameters:

$\displaystyle\frac{[k_{Q}(\mu_{w})]_{ds}}{V_{td}^{*}V_{ts}}=\tilde{c}_{W}(\Lambda)$
$\displaystyle\left[\frac{1}{6}I_{1}+\frac{\alpha_{t}(\mu_{w})}{4\pi}f(x_{t})I_{2}-\frac{\alpha_{t}(\mu_{w})}{4\pi}\frac{3\alpha(\mu_{w})}{2\pi
s_{w}^{2}}h(x_{t})\right]+$ $\displaystyle+\tilde{c}_{B}(\Lambda)$
$\displaystyle\left[\frac{1}{6}I_{3}+\frac{\alpha_{t}(\mu_{w})}{4\pi}f(x_{t})I_{4}\right]$
$\displaystyle-c_{\phi}(\Lambda)$
$\displaystyle\left[-\frac{1}{6}I_{5}+\frac{1}{9}A(\mu_{w},\Lambda)+\frac{\alpha_{t}(\mu_{w})}{4\pi}f(x_{t})\left[1-A(\mu_{w},\Lambda)-I_{6}\right]\,,\right]$
(19)

where $\alpha_{t}\equiv y_{t}^{2}/(4\pi)$. We are now in a position to
evaluate the effective coupling numerically. Setting $\Lambda=1~{}\rm{TeV}$,
and substituting Eqs. (14) into Eq. (2.1) we find

$\displaystyle\frac{[k_{Q}(\mu_{w})]_{ds}}{V_{td}^{*}V_{ts}}\bigg{|}_{\Lambda=\rm{1TeV}}\simeq-9.7\times
10^{-3}c_{W}(\Lambda)+8.2\times 10^{-3}c_{\phi}(\Lambda)-3.5\times
10^{-5}c_{B}(\Lambda)\,.$ (20)

From this expression we see that the contribution from the $\mathcal{O}_{B}$
operator to the effective coupling is strongly suppressed, so the production
mechanism will be dominated by the $\mathcal{O}_{\phi}$ and $\mathcal{O}_{W}$
operators. Furthermore Eq. (20) is the leading contribution to the amplitude,
Eq. (7), being much larger than $N_{8}\sim 10^{-7}$. No other flavour-changing
coupling is generated, in particular $[k_{q}(\mu_{w})]_{ds}=0$. Moreover, once
the EW bosons have been integrated out of the theory there is no additional
running of the coupling constants, and we can therefore take
$[k_{Q}(\mu_{w})]_{ds}$ at energies below the EW scale, as discussed in Ref.
Bauer:2020jbp . Hence to a good approximation we have

$i\mathcal{A}_{K\to\pi
a}\simeq-\frac{m_{K}^{2}-m_{\pi}^{2}}{2f}[k_{Q}(\mu_{w})]_{ds}\,.$ (21)

The decay width obtained from this amplitude reads

$\displaystyle\Gamma(K^{+}\to\pi^{+}a)$
$\displaystyle=\dfrac{m_{K}^{3}\big{|}[k_{Q}(\mu_{w})]_{sd}\big{|}^{2}}{64\pi}\lambda_{\pi
a}^{1/2}\left(1-\dfrac{m_{\pi}^{2}}{m_{K}^{2}}\right)^{2}\,,$ (22)

with the definitions

$\displaystyle\lambda_{\pi a}$
$\displaystyle\equiv\lambda(1,m_{a}^{2}/m_{K}^{2},m_{\pi}^{2}/m_{K}^{2})\,,$
(23) $\displaystyle\lambda(a,b,c)$
$\displaystyle=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc\,.$ (24)

Other recent estimates of this width (including tree-level and loop
corrections) have been presented in Refs. Izaguirre:2016dfi ; Gavela:2019wzg ;
Guerrera:2021yss . In Refs. Izaguirre:2016dfi ; Gavela:2019wzg , the one-loop
correction to $[k_{Q}]_{ds}$ from the operators ${\mathcal{O}}_{W}$ and
${\mathcal{O}}_{\phi}$ is assumed to be the dominant contribution. The
physical amplitude mediated by this coupling is written in terms of the scalar
form factor between a kaon and pion, and extracted from lattice QCD (see Refs.
Izaguirre:2016dfi ; Gavela:2019wzg  for further details). This result is
equivalent to that obtained in $\chi$PT as long as the effective coupling
$[k_{Q}]_{ds}$ is the same, and if the form factor is set to one333In the case
of $B\to Ka$ decays, the decay width can be obtained from Eq. (22) replacing
the CKM matrix elements and meson masses by the corresponding ones in this
case, and noting that the relevant form factor is $f_{0}^{B\to
K}(q^{2}=0)\simeq 0.37$ (see _e.g.,_ Fig. 33 in Ref. Aoki:2021kgd )., which is
a good approximation. The result quoted in Ref. Gavela:2019wzg  for the
effective coupling is

$\displaystyle\frac{[k_{Q}(\mu_{w})]_{ds}}{V_{td}^{*}V_{ts}}\bigg{|}_{\Lambda=\rm{1TeV}}\simeq-{g_{2}(\mu_{w})^{2}\over
16\pi^{2}}\left(3c_{W}(\Lambda)g(x_{t})-{c_{\phi}(\Lambda)\over
4}x_{t}\log{\Lambda^{2}\over\mu_{w}^{2}}\right),$ (25)

where $g(x)=xh(x)$, and $g_{2}(\mu_{w})$ is the weak gauge coupling evaluated
at the weak scale. Using Eq. (14) and the relations
$m_{t}=y_{t}v/\sqrt{2},\alpha_{2}=\alpha/s_{w}^{2},v=2m_{W}/g_{2}$, it is
straightforward to show that Eq. (26) is equivalent to the dominant
contributions in Eq. (2.1) up to a global sign:

$\displaystyle\frac{[k_{Q}]_{ds}}{V_{td}^{*}V_{ts}}\bigg{|}_{\Lambda=\rm{1TeV}}\simeq-\tilde{c}_{W}(\Lambda)$
$\displaystyle\left[\frac{\alpha_{t}(\mu_{w})}{4\pi}\frac{3\alpha}{2\pi
s_{w}^{2}}g(x_{t})\right]-c_{\phi}(\Lambda)\left[\frac{1}{9}A(\mu_{w},\Lambda)\right],$
(26)

since (see App. A)

$A(\mu,\Lambda)\simeq{9g_{2}^{2}\over
64\pi^{2}}x_{t}\log{\Lambda^{2}\over\mu^{2}}\,.$ (27)

Finally, in Ref. Guerrera:2021yss  the authors considered additional
contributions to $K\to\pi a$ at tree level, stemming from a non-universal ALP
coupling to fermions. Even in that case the one-loop contributions are found
to dominate over the tree-level diagrams, as long as the ALP coupling to the
top quark in the loop is non-vanishing. It should also be noted that both the
tree-level and one-loop diagrams (plus additional contributions not considered
in from Ref. Guerrera:2021yss ) are accounted for in the calculation using
$\chi$PT, Eq. (7).

###  2.2 ALP decay channels

For ALP masses below 400 MeV, the decay channels that are kinematically open
are $a\to\gamma\gamma$, $a\to e^{+}e^{-}$ and $a\to\mu^{+}\mu^{-}$. The decay
width into leptons can be written as:

$\Gamma(a\to\ell^{+}\ell^{-})=|c_{\ell\ell}|^{2}\dfrac{m_{a}m_{\ell}^{2}}{8\pi
f_{a}^{2}}\sqrt{1-\dfrac{4m_{\ell}^{2}}{m_{a}^{2}}}\,,$ (28)

where at low energies ($\mu\sim 2~{}\mathrm{GeV}$) $c_{\ell\ell}$ is given at
one loop by Gavela:2019wzg ; Bauer:2017ris

$\displaystyle\begin{split}c_{\ell\ell}=c_{\phi}&+\frac{3\,\alpha}{4\pi}\left(\frac{3\,c_{W}}{s_{w}^{2}}+\frac{5\,c_{B}}{c_{w}^{2}}\right)\log\dfrac{f_{a}}{m_{W}}+\dfrac{6\,\alpha}{\pi}\left(c_{B}\,c_{w}^{2}+c_{W}\,s_{w}^{2}\right)\log\dfrac{m_{W}}{m_{\ell}}\,,\end{split}$
(29)

and to simplify the notation we have written $c_{i}\equiv c_{i}(\Lambda)$.
Similarly, the decay width into two photons reads

$\displaystyle\Gamma(a\to\gamma\gamma)$
$\displaystyle=|c_{\gamma\gamma}|^{2}\dfrac{m_{a}^{3}}{4\pi f_{a}^{2}}\,,$
(30)

where the effective coupling at low energies is given at one loop by
Gavela:2019wzg ; Bauer:2017ris

$\displaystyle\begin{split}c_{\gamma\gamma}=&c_{W}\,\Big{[}s_{w}^{2}\,+\frac{2\,\alpha}{\pi}B_{2}(\tau_{W})\Big{]}+c_{B}\,c_{w}^{2}-c_{\phi}\,\frac{\alpha}{4\pi}\,\bigg{(}B_{0}-\frac{m_{a}^{2}}{m_{\pi}^{2}-m_{a}^{2}}\bigg{)}\,.\end{split}$
(31)

Here, $B_{0}$ and $B_{2}$ are loop functions (which can be found in App. B),
and $\tau_{W}=4m_{W}^{2}/m_{a}^{2}$.

A comparison of the different terms entering the effective couplings
$c_{\gamma\gamma}$ and $c_{\ell\ell}$ allows to see that:

  * •

the partial width $\Gamma(a\to\gamma\gamma)$ depends exclusively on the mass
of the ALP, while $\Gamma(a\to\ell^{+}\ell^{-})$ depends also on the lepton
mass. This means that, for similar values of the effective couplings
$c_{\gamma\gamma}$ and $c_{\ell\ell}$, then
$\Gamma(a\to\gamma\gamma)\gg\Gamma(a\to\ell^{+}\ell^{-})$;

  * •

once the decay channel $a\to\mu^{+}\mu^{-}$ is open, it will completely
dominate over the decay channel $a\to e^{+}e^{-}$, due to the much larger muon
mass;

  * •

if $c_{\phi}\gg c_{W},c_{B}$ the ALP will predominantly decay into lepton
pairs since the terms in $c_{\ell\ell}$ which are proportional to
$c_{B},c_{W}$ are suppressed by $\alpha$. Conversely, if $c_{W}\gg c_{\phi}$
then the ALP will decay mostly into photons;

  * •

if loop-corrections are neglected, the effective coupling for
$\Gamma(a\to\gamma\gamma)$ reads $c_{\gamma\gamma}\sim
s_{w}^{2}c_{W}+c_{w}^{2}c_{B}$. Thus, the operator $\mathcal{O}_{B}$ can have
a significant impact on (and even lead) the decay width for this channel,
since $s_{w}^{2}<c_{w}^{2}$.

##  3 Previous constraints

In this section we provide a brief summary of the most relevant constraints
for ALPs in the mass range $2m_{e}<m_{a}<m_{K}-m_{\pi}$, and explain how these
have been recasted to the specific scenario studied in this work. For more
detailed compilations of bounds, see _e.g.,_ Refs. Dolan:2017osp ;
Dobrich:2019dxc ; Essig:2010gu ; Bauer:2017ris ; Bauer:2021mvw ;
Izaguirre:2016dfi ; Jaeckel:2015jla ; Goudzovski:2022vbt .

### Visible ALP decays

At fixed-target experiments, the ALP can be produced through Primakoff
scattering (that is, the conversion of a photon into an ALP in the vicinity of
a nucleus), or through its mixing with pseudoscalar mesons produced in the
target ($\pi^{0},\eta,\eta^{\prime}$). If it couples directly to quarks or
electrons, it can also be produced by proton or electron Bremsstrahlung. The
analyses performed in the literature typically assume that the ALP is coupled
predominantly to either photons or electrons and therefore the leading
production and decay mechanisms will be different in the two scenarios. If the
ALP predominantly couples to photons, the strongest constraints below 500 MeV
come from the E137 experiment Bjorken:1988as , see _e.g.,_ Refs. Dolan:2017osp
; Dobrich:2019dxc ; Essig:2010gu . We take the limit from Ref. Dolan:2017osp ,
where the E137 bound was recasted on the ALP coupling to photons after EW
symmetry breaking, $g_{a\gamma\gamma}$. Since they were obtained assuming that
the ALP couples predominantly to photons, they can be directly applied444The
bound is rescaled as $g_{a\gamma\gamma}/4=c_{W}s_{w}^{2}$, considering the
different normalization used in Ref. Dolan:2017osp  for the $\mathcal{O}_{W}$
operator. to the $c_{W}$ coupling. Conversely, if $\mathcal{O}_{\phi}$
dominates the ALP couples mostly to electrons. In this case we use the
analysis of Ref. Essig:2010gu , which computed the excluded regions for E137
and CHARM CHARM:1985anb  assuming the ALP is predominantly produced through
its mixing with pseudoscalar mesons, and that it decays into $e^{+}e^{-}$
pairs. Here we use their CHARM bound, which is stronger than the one obtained
from E137.

At higher masses, significant constraints are also obtained from CHARM
CHARM:1985anb  and LHCb LHCb:2015nkv ; LHCb:2016awg , assuming that the ALP is
produced from $B$-meson decays (via $B\to Ka$) and that it decays visibly
within the detector as $a\to\mu^{+}\mu^{-}$. Note that, although the obtained
constraints in this case are less stringent they are applicable in a wider
mass range, $2m_{\mu}<m_{a}<m_{B}-m_{K}$. We take these from Refs.
Dobrich:2018jyi ; Gavela:2019wzg .

Finally, at electron-positron colliders ALPs could be produced through an off-
shell photon (for example, via $e^{+}e^{-}\to\gamma^{*}\to\gamma a$) or in
photon fusion ($e^{+}e^{-}\to e^{+}e^{-}a$) and decay to two photons. The
strongest bound in this case is obtained from a reanalysis of LEP data
Jaeckel:2015jla .

### Invisible ALP decays.

If the ALP is sufficiently long-lived (for very light masses, or small enough
couplings) it may exit the detector without decaying. This would lead to an
excess of decays of kaons into pions plus missing energy
($K\to\pi+\mathrm{inv}$), or $B$ mesons to kaons plus missing energy ($B\to
K+\mathrm{inv}$), which can be constrained using data from
$K\to\pi\nu\bar{\nu}$ or $B\to K\bar{\nu}\nu$ searches Essig:2010gu ;
Izaguirre:2016dfi ; Gavela:2019wzg . In this work, unlike in Refs.
Izaguirre:2016dfi ; Gavela:2019wzg  we do not assume significant ALP couplings
to the dark sector and, therefore, we need to take into account that if the
ALP is sufficiently short-lived it will not contribute to the mentioned
observables. Thus, an upper bound is imposed on $\text{BR}(M\to
M^{\prime}a)\times P_{\text{exit}}$, where $M,M^{\prime}$ stand for the parent
and daughter mesons, respectively, and $P_{\text{exit}}\sim
e^{-\Delta\ell_{\mathrm{det}}/L_{a}}$ is the probability of the ALP to exit
the detector without decaying. Here $\Delta\ell_{\mathrm{det}}$ stands for the
approximate detector size and $L_{a}\equiv\gamma_{a}\beta_{a}c\tau_{a}$ is the
ALP decay length in the lab frame ($\tau_{a}$ denotes the proper lifetime of
the ALP and $\gamma_{a},\beta_{a}$ correspond to the boost variables, while
$c$ is the speed of light).

The strongest limits on $K\to\pi a$ are obtained from the NA62 experiment
NA62:2020xlg ; NA62:2020pwi , where we assume $\Delta\ell_{\text{det}}\simeq
100\;\text{m}$ as detector size and $p_{K}\simeq 75~{}\mathrm{GeV}$ as the
kaon momentum NA62:2020pwi . However, the very competitive bounds from E787 &
E949 BNL-E949:2009dza  are comparable (and even dominate) in the region close
to the pion mass. In this case we assume the kaon decays at rest, and we take
$\Delta\ell_{\text{det}}\sim 1.5~{}\mathrm{m}$. In the case of $B\to Ka$
decays, we saturate the upper bound from the Belle experiment on $B^{+}\to
K^{+}\nu\bar{\nu}$ Belle:2017oht , taking $\Delta\ell_{\mathrm{det}}\simeq 5\
\mathrm{m}$ Belle:2012iwr  and $p_{a}\simeq 2.5\ \mathrm{GeV}$.

In principle, searches for mono-photon signals at colliders (_e.g.,_ at BaBar
BaBar:2008aby ; BaBar:2017tiz  or LEP L3:1997exg ) may be reinterpreted in the
context of an ALP that decays outside the detector Izaguirre:2013uxa ;
Essig:2013vha ; Izaguirre:2016dfi . However, these are milder than the rest of
the limits considered in this work.

### ALP contributions to kaon three-body decays.

A third set of constraints can be derived from the decays $K\to\pi ee$ or
$K\to\pi\gamma\gamma$. An ALP would contribute to these through the same
penguin diagram leading to $K\to\pi$, with an extra vertex coupling the ALP
propagator (which can be either on-shell Izaguirre:2016dfi  or off-shell
Gavela:2019wzg ) to either two photons or an electron-positron pair. However,
such constraints are superseded by the ones from invisible ALP decays outlined
above, and will not be considered here.

### Astrophysical bounds.

Three main bounds are obtained from supernovae data: (1) from the requirement
that the energy loss induced by ALP emission does not exceed the energy loss
from neutrino emission Raffelt:1990yz ; Raffelt:1987yt ; (2) from the visible
signal resulting from the ALP burst, in the case where the ALP decays into
pairs of photons Jaeckel:2017tud ; and (3) from the observation of low-
luminosity core-collapse supernovae, which constrains the total energy
deposition in the progenitor star from radiative ALP decays, such as
$a\to\gamma\gamma$ or $a\to e^{+}e^{-}$ Caputo:2022mah . In this work we focus
on laboratory experiments, and refer the interested reader to Refs.
Chang:2018rso ; Croon:2020lrf ; Caputo:2022mah  for recent updates on
astrophysical constraints.

##  4 Axion-like particles at MicroBooNE

The results from Sec. 2 can be used to compute the branching ratio into the
$K\to\pi a$ channel as well as the differential angular distribution of the
ALP flux.

In order to compute the expected number of events at MicroBooNE from the NuMI
absorber, we use the same NuMI kaon distributions as in Ref. Coloma:2015pih ,
which were derived from a Monte Carlo simulation of the NuMI target when
exposed to a 120 GeV proton beam. As outlined in the introduction, the
expected ALP flux from the NuMI absorber would correspond to that of kaons
decaying at rest. Since the production takes place via a two-body decay an
isotropic flux is expected, with energy
$E_{a}=(-m_{\pi}^{2}+m_{a}^{2}+m_{K}^{2})/(2m_{K})$.

For a total number of $N_{K}$ kaon decays, the event rate expected from ALP
decays into $e^{+}e^{-}$ pairs inside the MicroBooNE detector can be computed
as

$N_{events}=\frac{N_{K}\times\text{BR}(K\to\pi a)}{4\pi}\;\text{BR}(a\to
e^{+}e^{-})\;\epsilon_{eff}\int_{\Delta\Omega_{\text{det}}}d\Omega~{}P_{\text{decay}}(\Omega)\,,$
(32)

where the integral runs over all trajectories with solid angle $\Omega$
intersecting the detector, and $\Delta\Omega_{\text{det}}$ is the solid angle
of the detector as seen from the absorber. Here $\epsilon_{eff}$ stands for
the detection efficiency (which depends on $m_{a}$) and $\text{BR}(a\to
e^{+}e^{-})$ is the branching ratio for the ALP to decay into an electron-
positron pair. Moreover, we have assumed an isotropic flux of kaons as well as
a point-like absorber. Finally, $P_{\text{decay}}$ represents the probability
of an ALP to decay inside the detector:

$P_{\text{decay}}=e^{-\frac{\ell_{\text{det}}}{L_{a}}}\left[1-e^{-\frac{\Delta\ell_{\text{det}}}{L_{a}}}\right]\,,$
(33)

where $\ell_{\text{det}}$ is the distance traveled before it reaches the
detector, and $\Delta\ell_{\text{det}}$ is the length of the ALP trajectory
intersecting the detector. In practice, both $\ell_{\text{det}}$ and
$\Delta\ell_{\text{det}}$ depend on the angle of the ALP trajectory. In this
work we perform a numerical integration over all ALP trajectories, using the
same detector dimensions and orientation with respect to the absorber as in
Ref. MicroBooNE:2021usw . However, for estimation purposes it is useful to
neglect the dependence of $P_{\text{decay}}$ with the solid angle and write
the expected number of events as

$N_{events}\sim N_{K}\times\text{BR}(K\to\pi
a)\times\varepsilon_{\text{det}}\times\epsilon_{eff}\times\langle
P_{\text{decay}}\rangle\times\text{BR}(a\to e^{+}e^{-})\,,$ (34)

where $\langle P_{\text{decay}}\rangle$ represents the average probability of
an ALP to decay inside the detector (obtained assuming as typical values
$\ell_{\text{det}}\sim 100$ m and $\Delta\ell_{\text{det}}\sim 4$ m), while
$\varepsilon_{\text{det}}=\Delta\Omega_{\text{det}}/(4\pi)\sim 1.5\times
10^{-3}$ stands for the detector geometric acceptance. Assuming that the
$\mathcal{O}_{\phi}$ coupling dominates, the branching ratio into electrons
dominates, $BR(a\to e^{+}e^{-})\simeq 1$, see Sec. 2.2. Regarding the total
number of kaons produced, numerically we find $N_{K}\sim 0.2N_{\text{PoT}}$
kaons produced in the NuMI absorber, where $N_{\text{PoT}}$ stands for the
number of Protons on Target (PoT) considered. According to this estimate, we
find that a few signal events are a priori expected for $m_{a}=100$ MeV and
$c_{\phi}\sim 2\times 10^{-3}$, for an exposure of $N_{\text{PoT}}\sim 2\times
10^{20}$. As we will see below, this agrees very well with the results from
the exact numerical calculation.

The rest of this section is structured as follows. First, we present the
current constraints derived from MicroBooNE data in Sec. 4.1. In Sec. 4.2 we
proceed and determine the future sensitivity to this scenario using the whole
collected (unprocessed) data, and consider additional search channels in
addition to $a\to e^{+}e^{-}$. Finally, in Sec. 4.3 we study possible
interference effects when two operators are simultaneously included in the
analysis.

###  4.1 Bounds from current MicroBooNE data

In order to determine the limit on the effective Wilson coefficients, we
follow the same approach as in Ref. MicroBooNE:2021usw  and perform an
unbinned $\chi^{2}$ analysis, taking $n_{bg}=1.9$ as the expected total number
of background events and $n_{\text{obs}}=1$ as the observed number of
events555The collaboration observed two candidate events in the dataset
passing the cuts; however, one of them was rejected.. We also include two
nuisance parameters in order to account for systematic uncertainties affecting
the overall normalization of the events for the signal ($\xi_{s}$) and
background ($\xi_{b}$). For each uncertainty a penalty term is added to the
$\chi^{2}$, which is then minimized over the nuisance parameters. Our
$\chi^{2}$ reads:

$\chi^{2}=\text{min}_{\xi_{s},\xi_{b}}\left\\{2(n_{\text{pred}}-n_{\text{obs}})+2n_{\text{obs}}\log\frac{n_{\text{obs}}}{n_{\text{pred}}}\right\\}+\frac{\xi_{s}^{2}}{\sigma_{s}^{2}}+\frac{\xi_{b}^{2}}{\sigma_{b}^{2}}\,,$
(35)

where we take $\sigma_{s}=30\%$ for the signal and $\sigma_{b}=89\%$ for the
backgrounds (according to Tab. I in Ref. MicroBooNE:2021usw ). Here
$n_{\text{pred}}$ is the predicted number of events

$n_{\text{pred}}(\\{c\\},m_{a},\xi_{s},\xi_{b})=\left[1+\xi_{s}\right]\,n_{sig}(\\{c\\},m_{a})+\left[1+\xi_{b}\right]\,n_{bg}\,,$
(36)

where $n_{sig}$ stands for the signal predicted by the model for an ALP with
mass $m_{a}$ and a given set of operator coefficients $\\{c\\}$. When
computing $n_{sig}$ we take the selection efficiencies provided in the
supplemental material of Ref. MicroBooNE:2021usw  as a function of the scalar
mass. We believe these efficiencies are directly applicable to our model given
that the final topology for the decay is exactly the same.

Our derived constraints using the $\chi^{2}$ definition in Eq. (35) are
presented in Fig. 1, under the assumption that the production branching ratio
and the decay width are completely independent.

Figure 1:  MicroBooNE bounds (at 95% C.L., for 2 d.o.f.) in the plane of
$c\tau$ vs the product of the branching ratios for the ALP production and
decay, for two different ALP masses as indicated by the labels (blue regions).
Our bounds are obtained from the recasting of the analysis in Ref.
MicroBooNE:2021usw , which used $1.93\times 10^{20}$ Protons on Target (PoT)
and looked for $e^{+}e^{-}$ pairs. For comparison we also show the constraints
(at 90% CL) from invisible decays ($B\to K+\mathrm{inv}$ Belle:2017oht  and
$K\to\pi+\mathrm{inv}$ NA62:2020xlg ; NA62:2020pwi , as well as the
constraints (at 95% CL) from visible ALP decays into $\mu^{+}\mu^{-}$ at LHCb
LHCb:2015nkv ; LHCb:2016awg  and CHARM CHARM:1985anb ; CHARM:1983ayi , taken
from Ref. Dobrich:2018jyi .

Although in a given model the two quantities are typically correlated, this
allows us to derive a model-independent constraint so that our results can be
easily recasted for other scenarios. Our results are shown for two
illustrative values of the mass of the ALP, below (left panel) and above the
kaon mass (right panel). The MicroBooNE region (shown in blue) shows the
portion of parameter space where we obtain a $\Delta\chi^{2}>5.99$,
corresponding to 95% confidence level (CL) for two degrees of freedom
(d.o.f.).

Figure 2:  MicroBooNE bounds (at 95% C.L., for 2 d.o.f) on $c_{\phi}$, for
$f_{a}=1\;\text{TeV}$ and setting $c_{W}=c_{B}=0$. These are obtained from the
recasting of the analysis in Ref. MicroBooNE:2021usw , which used $1.93\times
10^{20}$ Protons on Target (PoT). For comparison we also show the constraints
from the NA62 experiment NA62:2021zjw ; NA62:2020xlg , E787 & E949
BNL-E949:2009dza , CHARM Essig:2010gu ; CHARM:1985anb  and LHCb LHCb:2015nkv ;
LHCb:2016awg ; Gavela:2019wzg . The shaded vertical band is excluded due to
the large $a-\pi^{0}$ mixing, and is taken from Ref. Gavela:2019wzg .

In a given model, however, the values of the branching ratios and the lifetime
of the ALP will be correlated since they both will depend on the same set of
couplings as well as on the mass of the ALP. This typically leads to decrease
in sensitivity for experiments searching for visible ALP decays, as is the
case of MicroBooNE in this work. For example, from Fig. 1 we see that
MicroBooNE would be in principle sensitive to values of the production
branching ratio of the ALP at the level of $\mathcal{O}(2\times 10^{-11})$.
Taking into account the expression for the branching ratio in Eq. (22), for a
mass of 50 MeV this can be obtained setting $c_{W}\sim\mathcal{O}(5\times
10^{-4})$; however, this leads to very large values of $c\tau\sim 1500$ m (and
even larger values if we assume that the ALP is coupled via the
$\mathcal{O}_{\phi}$ operator instead of $\mathcal{O}_{W}$). In this regime
the expected number of decays in the detector scales as
$\Delta\ell_{\text{det}}/L_{a}$, see Eq. (33), suppressing the sensitivity for
very large values of $L_{a}$. Therefore, once we impose this restriction the
relative comparison between the sensitivity of visible and invisible decay
searches changes dramatically, as shown in Fig. 2. However, in spite of this
we still find that MicroBooNE is generally competitive with current
constraints, and is able to improve over these for ALP masses close to the
pion mass.

###  4.2 Future sensitivity

The limits shown in Sec. 4.1 were obtained using only a reduced subset of the
available NuMI data, and only searching for decays into $e^{+}e^{-}$ pairs.
However, according to Ref. MicroBooNE:2021usw , MicroBooNE has collected about
10 times more NuMI data, which has not been processed yet. Moreover, the
MicroBooNE detector consists on a Liquid Argon Time Projection Chamber (LAr
TPC), which allows to study not only final state topologies involving
$e^{+}e^{-}$ pairs, but also $\gamma\gamma$ and $\mu^{+}\mu^{-}$ pairs. Here
we compute the expected sensitivity using the full NuMI dataset recorded at
MicroBooNE (corresponding to $\approx 2.2\times 10^{21}$ PoT), and for three
different searches using the three possible final states for the decay of the
ALP.

Figure 3:  MicroBooNE sensitivity projections for $c_{\phi}$ (left panel) and
$c_{W}$ (right panel), assuming $f_{a}=1\;\text{TeV}$, considering all
(processed and unprocessed) NuMI data recorded at MicroBooNE
($N_{\text{PoT}}=2.2\times 10^{21}$). Results are shown separately for the
decay channels with the largest branching ratios in each case:
$a\to\gamma\gamma$, $a\to e^{+}e^{-}$ and $a\to\mu^{+}\mu^{-}$, as indicated
by the labels. Solid lines assume 89% and 30% signal and background systematic
uncertainties, while dotted lines assume 20% uncertainties for both. For
comparison, the shaded areas show the current bounds from E137 Bjorken:1988as
; Dolan:2017osp , NA62 NA62:2021zjw ; NA62:2020xlg , E787 & E949
BNL-E949:2009dza , LHCb LHCb:2015nkv ; LHCb:2016awg ; Gavela:2019wzg , CHARM
CHARM:1985anb ; Essig:2010gu  and LEP Jaeckel:2015jla , as well as our bound
derived from MicroBooNE data in Sec. 4.1. The shaded vertical band is
disfavored due to the large $a-\pi^{0}$ mixing, and is taken from Ref.
Gavela:2019wzg . In each panel only one operator is switched on at a time,
setting the remaining operator coefficients to zero.

The same $\chi^{2}$ definition as in the previous section, Eq. (35), will be
used here. In this case we assume that the experiment will observe a number of
events fully consistent with the background expectation
($n_{\text{obs}}=n_{bg}$). The number of background events is estimated
rescaling the predicted background in Ref. MicroBooNE:2021usw  by a factor of
$\approx 10$ in order to account for the much larger exposure, and we assume
this number to be the same for $e^{+}e^{-}$, $\mu^{+}\mu^{-}$ and
$\gamma\gamma$ events. Additionally, we use the same same detection
efficiencies as in Ref.MicroBooNE:2021usw  up to $m_{a}=210$ MeV, and a flat
efficiency at 9.6% for higher masses. We believe that this is a good estimate
for the $\gamma\gamma$ and $\mu^{+}\mu^{-}$ decay channels, given the good
resolution of the MicroBooNE detector and that the signal follows a two-body
decay topology in all cases.

Figure 3 shows the expected sensitivity from a search using all recorded NuMI
data at MicroBooNE. Our projections are shown for the ALP decay channels
yielding the largest branching ratios in each case: $a\to e^{+}e^{-}$ and
$a\to\mu^{+}\mu^{-}$ in the case of the $\mathcal{O}_{\phi}$ operator (left
panel), and $a\to\gamma\gamma$ in the case of the $\mathcal{O}_{W}$ operator
(right panel), see Sec. 2.2. For comparison, the regions of parameter space
disfavored by current experiments (including the MicroBooNE bound from Fig. 2)
are indicated by the shaded areas, while the colored lines indicate the future
sensitivity contours for MicroBooNE at 95% CL (2 d.o.f.). Note that, while the
$\mathcal{O}_{\phi}$ operator would lead to $\ell^{+}\ell^{-}$ signals,
sensitivity to the $\mathcal{O}_{W}$ is possible only from a search in the
$\gamma\gamma$ channel since the contribution of this operator to
$c_{\ell\ell}$ is suppressed by $\alpha$, see Eq. (29). This illustrates the
importance of conducting searches with different final state topologies,
making full use of the LAr TPC capabilities. Also, note the $\sim 10$ times
better sensitivity to $c_{W}$ compared to $c_{\phi}$. Although the two
couplings enter on equal terms at the production level, the decay width for
$c_{W}$ is larger by roughly ${\mathcal{O}}(m_{a}^{2}/m_{l}^{2})$ compared to
that of $c_{\phi}$ leading to a shorter ALP lifetime. This eventually enhances
the event rate in the detector, since in the limit of small couplings the
decay probability can be approximated as
$P_{\text{decay}}\propto\Delta\ell_{\text{det}}/L_{a}$, see Eq. (33).

Due to the very large exposure considered, we find that the bound is no longer
statistics-limited but limited by systematics instead. While the results in
Ref. MicroBooNE:2021usw  were obtained for very conservative systematic
uncertainties, it may be possible to reduce these considerably. For example,
their background uncertainties were dominated by the detector modeling (and,
in particular, by the low statistics of simulated neutrino events passing
selection cuts), while cross section uncertainties may be further reduced with
dedicated studies and improvement of theoretical models. Therefore, Fig. 3
shows two sets of lines for each decay channel: solid lines, corresponding to
the same set of uncertainties as in Fig. 2
($\sigma_{s}=30\%,\sigma_{b}=89\%$); and dotted lines, corresponding to
$\sigma_{s}=\sigma_{b}=20\%$. While in the former case the improvement over
the present MicroBooNE constraint for $c_{\phi}$ would only be marginal (left
panel), from the comparison of the solid and dotted lines it can be seen that
there is room for improvement provided that systematic uncertainties can be
significantly reduced below current values. Additionally, our results show the
potential to constrain an unexplored region of parameter space for $c_{\phi}$
by looking for ALP decays into $\mu^{+}\mu^{-}$ final states. Finally, for
$c_{W}$ (right panel) we find that the bounds from MicroBooNE would be
competitive with those obtained from beam-dump experiments and may even be
able to improve over these, depending on the level of systematic uncertainties
assumed.

###  4.3 Interference effects between different operators

So far our results have been obtained under the assumption that only one of
the operators which control the production of the ALP ($\mathcal{O}_{W}$, or
$\mathcal{O}_{\phi}$) is present in the effective Lagrangian and, for
simplicity, the coefficient that goes with $\mathcal{O}_{B}$ has also been set
to zero. However, it has been pointed out previously in the literature
Gavela:2019wzg  that destructive interference between different EW operators
could significantly affect the effective couplings controlling both the ALP
production and decay. Thus, in Fig. 4 we compute the expected sensitivity for
the full data set when two operators are switched on simultaneously. In each
row the left (right) panel shows to the result when the relative sign between
the two operators introduced is positive (negative), since this either allows
or forbids a destructive interference for the production and decay mechanisms,
depending on the operators considered. For concreteness, these results are
obtained assuming 20% systematic uncertainties on both signal and background,
and for a fixed mass of the ALP ($m_{a}=200$ MeV).

Figure 4:  MicroBooNE sensitivity projections, assuming two couplings are non-
zero at a time. The shaded regions show the sensitivity of the experiment at
95% CL (2 d.o.f.), from a search for a $\gamma\gamma$ (red) or $e^{+}e^{-}$
(blue) signal. Results are shown for $f_{a}=1\;\text{TeV}$ and
$m_{a}=200~{}\mathrm{MeV}$, for $N_{\text{PoT}}=2.2\times 10^{21}$, and
assuming $\sigma_{sig}=\sigma_{bg}=20\%$. Left (right) panels are obtained
assuming a positive (negative) relative sign between the operators introduced.
The diagonal lines indicate destructive interference between different
operators entering the production rate (dotted lines) or decay width (dashed
lines), see text for details.

The first row in Fig. 4 shows the sensitivity when the two operators that
control the production of the ALP, $c_{\phi}$ and $c_{W}$, are simultaneously
present. As can be seen, in the limit $c_{W}\to 0$ ($c_{\phi}\to 0)$ we
recover the sensitivity to $c_{\phi}$ ($c_{W}$) given in Fig. 3, as expected.
In the right panel we see that, as $c_{\phi}$ increases, the sensitivity to
$c_{W}$ is enhanced because of the larger ALP production rate, see Eq. (15).
However, in the left panel two destructive interference patterns take place.
First, when the magnitude of the two coefficients is the same (along the
diagonal dotted line) the sensitivity is lost because they interfere
destructively in the ALP production rate. Secondly, for $c_{W}\sim
m_{a}^{2}/(m_{a}^{2}-m_{\pi}^{2})c_{\phi}\alpha/(4s_{w}^{2}\pi)\sim
2c_{\phi}\alpha/\pi$ (along the diagonal dashed line) the contributions from
the two operators to $c_{\gamma\gamma}$ cancel each other (see Eq. (31)),
leading to a suppressed decay width into the $a\to\gamma\gamma$ decay channel.
As a result, the sensitivity to the ALP using a di-photon signal is lost, but
the experiment would still be sensitive to $a\to e^{+}e^{-}$.

In the middle row of Fig. 4 we see the effect due to the interplay between
$\mathcal{O}_{B}$ and $\mathcal{O}_{\phi}$. As already mentioned, since the
$\mathcal{O}_{B}$ operator does not significantly affect the ALP production
rate, it will only have an impact on its lifetime. In the right panel the two
coefficients have opposite signs and, thus, no destructive interference arises
in the decay width. As we can see, the inclusion of the $\mathcal{O}_{B}$
operator allows for a better sensitivity to $c_{\phi}$ up to a factor of $\sim
10$, for searches using the $\gamma\gamma$ channel. This is due to the larger
decay width obtained in this case (which makes the ALP shorter-lived,
increasing the decay rate within the detector) combined with the larger
branching ratio into the $\gamma\gamma$ channel, see Eqs. (29) and (31). In
the left panel a similar effect can be seen, leading to an increased
sensitivity to $c_{\phi}$ for values of $c_{B}\sim\mathcal{O}(10^{-4})$.
However since in this case the two coefficients enter with the same sign in
Eq. (31) they can interfere destructively and suppress the ALP decay width
into $\gamma\gamma$, for values of $c_{B}\sim
m_{a}^{2}/(m_{a}^{2}-m_{\pi}^{2})c_{\phi}\alpha/(4c_{w}^{2}\pi)\sim
2c_{\phi}\alpha/(3\pi)$. Moreover, once $c_{B}$ gets too large the sensitivity
to $c_{\phi}$ is lost as the ALP decays too fast and the decay probability is
exponentially suppressed.

Finally, the lower row in Fig. 4 shows the interplay between $c_{B}$ and
$c_{W}$. In this case, the absence of the $\mathcal{O}_{\phi}$ operator
implies that the signal to $\ell^{+}\ell^{-}$ is heavily suppressed and the
sensitivity is obtained only in the $\gamma\gamma$ channel. In the left panel
no interference is expected since the two coefficients enter with the same
sign in $c_{\gamma\gamma}$; however, for values of $c_{B}\gtrsim 2\times
10^{-5}$ the ALP becomes too short-lived and the sensitivity is lost. In the
right panel, on the other hand, since the two couplings enter with a different
sign a destructive interference can arise in $c_{\gamma\gamma}$, leading to
longer lifetimes. Thus, sensitivity to larger values of $c_{W}$ and $c_{B}$ is
still possible along the line where $c_{W}s_{w}^{2}\sim c_{B}c_{w}^{2}$, see
Eq. (31).

##  5 Summary and conclusions

In spite of the strong experimental evidence pointing towards the existence of
BSM physics, our efforts to discover it at colliders and direct detection
experiments have been unfruitful so far. While it is possible that the new
physics is too heavy and lies outside of our reach at the LHC, an interesting
alternative is that the new physics is light but weakly coupled to the visible
sector, making it very elusive.

Neutrino experiments, counting on very massive detectors and powerful sources,
lie at the edge of the intensity frontier and are therefore well-suited to
search for weakly coupled light degrees of freedom. In this work we focused on
the ALP scenario, which is well-motivated from the theoretical point of view:
light pseudoscalars generally arise as pseudo-Nambu-Goldstone bosons of BSM
theories with spontaneous breaking of a global symmetry, and might address
some of the most relevant open questions in particle physics, such as the
strong CP problem or the origin of dark matter.

For concreteness, we have considered a set of higher-dimensional effective
operators coupling the ALP to the electroweak gauge bosons, which would lead
to ALP production from kaon decays, via $K\to\pi a$. Working in chiral
perturbation theory ($\chi$PT) allows to derive the effective coupling of ALP
to mesons at low energies in all generality. While this had been done in the
literature before, we have obtained the relevant expressions for the
particular set of operators considered in this work. In doing so, we have also
shown explicitly that previous calculations taking only into account the
penguin diagrams agree with the full result obtained in $\chi$PT for the case
under consideration here, as they capture the leading order contributions to
this process.

Next, we have recasted a recent MicroBooNE analysis MicroBooNE:2021usw  which
used data taken for the NuMI beam, searching for electron-positron pairs
pointing towards the NuMI absorber. These results can be directly applied to
our model, if the ALP is coupled to the SM through the $\mathcal{O}_{\phi}$
operator. Our results show that MicroBooNE data already sets competitive
bounds on this operator (comparable to those of NA62) for ALP masses between
100 and 200 MeV (Fig. 2). For completeness, we also present our results in the
plane of production branching ratio vs the ALP lifetime, without taking into
account that these are typically correlated within a given model (Fig. 1).
Being model-independent, this allows to easily recast our constraints to other
models (including a different set of ALP couplings), or even for a different
long-lived particle as long as its production and decay mechanisms are the
same as considered here. For optimal values of the lifetime of the long-lived
particle, current MicroBooNE data sets tight constraints on the production
branching ratio, $\mathrm{BR}(K\to\pi a)<\mathcal{O}(\mathrm{few}\times
10^{-11})$.

Finally, we have also computed the sensitivity using the full NuMI dataset
recorded at MicroBooNE (Fig. 3). Due to the excellent particle identification
and resolution capabilities of the LArTPC technology, we have presented our
sensitivities for three different searches as indicated, $a\to\mu\mu$, $a\to
ee$ and $a\to\gamma\gamma$. Our results show that, depending on the level of
systematic uncertainties assumed, MicroBooNE might be able to improve over
current constraints for masses in the range between 100 and 250 MeV. We point
out the complementarity among searches using different final state topologies,
which takes full advantage of the unique LAr TPC capabilities: while searches
for an excess in the $\ell^{+}\ell^{-}$ channels are mostly sensitive to the
$\mathcal{O}_{\phi}$ operator, searches for an excess in the di-photon channel
would be sensitive to the $\mathcal{O}_{W}$ operator instead. Finally, we also
explored the possible interference effects arising when two operators are
switched on simultaneously (Fig. 4). Here it should be noted that, while the
$\mathcal{O}_{B}$ operator does not induce ALP production from kaon decays,
MicroBooNE can be sensitive to its impact on the ALP decay rate. Additional
interference effects can take place in the production vertex between the
$\mathcal{O}_{\phi}$ and $\mathcal{O}_{W}$ operators.

In summary, this work stands out as a clear example of the multiple
capabilities of neutrino experiments to search for new physics, not only in
the neutrino sector but in other sectors as well. Needless to say, that the
type of analysis performed here may be applicable to other neutrino beam
experiments using near detectors. An obvious example is the case of the DUNE
experiment, which will also make use of the LAr TPC technology. In this case
the near detector will be placed on axis with respect to the direction of the
beam, resulting in a larger background level. However, the use of a gas TPC
(instead of a LAr TPC) would reduce it significantly; additionally, the
possibility to move the near detectors off-axis may allow to enhance the
signal sensitivity. A study of the DUNE sensitivity to this scenario cannot be
done without a careful assessment of the background levels and is left for
future work.

###### Acknowledgements.

We warmly thank Belen Gavela, Luca Merlo and Olcyr Sumensari for useful
discussions. This project has received funding/support from the European
Union’s Horizon 2020 research and innovation program under the Marie
Skłodowska-Curie grant agreement No 860881-HIDDeN, as well as from Grants
PID2019-108892RB-I00, PID2020-113644GB-I00 and and CEX2020-001007-S, funded by
MCIN/AEI/10.13039/501100011033. The authors acknowledge support from
Generalitat Valenciana through the plan GenT program (CIDEGENT/2018/019) and
PROMETEO/2019/083. The work of PC is supported by Grant RYC2018-024240-I
funded by MCIN/AEI/10.13039/501100011033 and by “ESF Investing in your
future”.

##  Appendix A Integrals

We define

$\displaystyle A(\mu,\mu^{\prime})\equiv 1-e^{-18U(\mu,\mu^{\prime})}\,,$ (37)
$\displaystyle
U(\mu,\Lambda)\equiv-\int_{\Lambda}^{\mu}\frac{d\mu^{\prime}}{\mu^{\prime}}\frac{y_{t}^{2}(\mu^{\prime})}{32\pi^{2}}\,.$
(38)

and the following integrals, which take into account the running of the SM
coupling constants with the energy scale:

$\displaystyle I_{1}(\mu_{w},\Lambda)$ $\displaystyle\equiv$
$\displaystyle\int_{\Lambda}^{\mu_{w}}\frac{d\mu^{\prime}}{\mu^{\prime}}A(\mu_{w},\mu^{\prime})\frac{3\alpha_{2}^{2}(\mu^{\prime})}{8\pi^{2}}\,,$
(39) $\displaystyle I_{2}(\mu_{w},\Lambda)$ $\displaystyle\equiv$
$\displaystyle\int_{\Lambda}^{\mu_{w}}\frac{d\mu^{\prime}}{\mu^{\prime}}\left[1-A(\mu_{w},\mu^{\prime})\right]\frac{9\alpha_{2}^{2}(\mu^{\prime})}{16\pi^{2}}\,,$
(40) $\displaystyle I_{3}(\mu_{w},\Lambda)$ $\displaystyle\equiv$
$\displaystyle\int_{\Lambda}^{\mu_{w}}\frac{d\mu^{\prime}}{\mu^{\prime}}A(\mu_{w},\mu^{\prime})\frac{17\alpha_{1}^{2}(\mu^{\prime})}{72\pi^{2}}\,,$
(41) $\displaystyle I_{4}(\mu_{w},\Lambda)$ $\displaystyle\equiv$
$\displaystyle\int_{\Lambda}^{\mu_{w}}\frac{d\mu^{\prime}}{\mu^{\prime}}\left[1-A(\mu_{w},\mu^{\prime})\right]\frac{17\alpha_{1}^{2}(\mu^{\prime})}{48\pi^{2}}\,,$
(42) $\displaystyle I_{5}(\mu_{w},\Lambda)$ $\displaystyle\equiv$
$\displaystyle\int_{\Lambda}^{\mu_{w}}\frac{d\mu^{\prime}}{\mu^{\prime}}A(\mu_{w},\mu^{\prime})A(\mu^{\prime},\Lambda)\left[\frac{8\alpha_{s}^{2}(\mu^{\prime})}{27\pi^{2}}+\frac{\alpha_{2}^{2}(\mu^{\prime})}{16\pi^{2}}+\frac{17^{2}}{54^{2}}\frac{3\alpha_{1}^{2}(\mu^{\prime})}{4\pi^{2}}\right],$
(43) $\displaystyle I_{6}(\mu_{w},\Lambda)$ $\displaystyle\equiv$
$\displaystyle\int_{\Lambda}^{\mu_{w}}\frac{d\mu^{\prime}}{\mu^{\prime}}\left[1-A(\mu_{w},\mu^{\prime})\right]A(\mu^{\prime},\Lambda)\left[\frac{4\alpha_{s}^{2}(\mu^{\prime})}{9\pi^{2}}+\frac{3\alpha_{2}^{2}(\mu^{\prime})}{32\pi^{2}}+\frac{17^{2}}{54^{2}}\frac{9\alpha_{1}^{2}(\mu^{\prime})}{8\pi^{2}}\right].$

In Eqs. (39)-(LABEL:eq:Is6), $\alpha_{1}\equiv\alpha/c_{w}^{2}$,
$\alpha_{2}\equiv\alpha/s_{w}^{2}$, $\alpha_{s}\equiv g_{s}^{2}/(4\pi)$
correspond to the different gauge coupling strength functions.

The running of the coupling strength functions is computed solving the RGE
equations $d\alpha_{i}(\mu)/d\ln\mu=-\beta^{(i)}_{0}\alpha_{i}^{2}/(2\pi)$,
where $\beta^{(i)}_{0}$ are the coefficients at one-loop order and
$\alpha_{i}\in\left\\{\alpha_{1},\alpha_{2},\alpha_{s},\alpha_{t}\right\\}$,
with $\alpha_{t}\equiv y_{t}^{2}/(4\pi)$. Specifically, we take the running of
$\alpha_{t}$ to follow that of $\alpha_{s}$ and therefore use
$\beta^{(1)}_{0}=41/6;\beta^{(2)}_{0}=-19/6;\beta^{(3)}_{0}=-7;\beta^{(t)}_{0}=-7$.
In doing this, the cross terms between $\alpha_{t}$ and $\alpha_{s}$ are
neglected. However, we believe this is a good approximation since do find a
good numerical agreement with the results in Ref. Bauer:2020jbp  where this
effect has been considered.

Finally, let us point out that using Eq. (37) in Ref. Bauer:2020jbp  it is
straightforward to show that $U(\mu,\Lambda)$ can be approximated as

$U(\mu,\Lambda)\simeq\frac{1}{64\pi^{2}}\frac{g_{2}(\mu)^{2}x_{t}}{2}\ln\frac{\Lambda^{2}}{\mu^{2}}\,,$
(45)

where we have used the relations $m_{t}=y_{t}v/\sqrt{2},v=2m_{W}/g_{2}$. Thus,
Eq. (27) follows directly from the substitution of Eq. (45) into Eq. (37).

##  Appendix B Loop functions

The loop functions in Eq. (31) read

$\displaystyle
B_{0}=\bigg{(}\sum_{{f\,=\,c,t}}N_{c}Q_{f}^{2}\,B_{1}(\tau_{f})-\sum_{f\,=\,b,\ell^{-}_{\alpha}}N_{c}Q_{f}^{2}\,B_{1}(\tau_{f})\bigg{)}$
(46)

where

$\begin{array}[]{l}B_{1}(\tau)=1-\tau\,f^{2}(\tau)\,,\\\
B_{2}(\tau)=1-(\tau-1)\,f^{2}(\tau)\,,\end{array}$ (47)

with

$\displaystyle
f(\tau)=\left\\{\begin{array}[]{ll}\arcsin\frac{1}{\sqrt{\tau}}\,;&~{}\tau\geq
1\,,\\\
\frac{\pi}{2}+\frac{i}{2}\ln\frac{1+\sqrt{1-\tau}}{1-\sqrt{1-\tau}}\,;&~{}\tau<1\,.\end{array}\right.$
(50)

Here, $\tau_{f}\equiv 4m_{f}^{2}/m_{a}^{2}$, $Q_{f}$ denotes the electric
charge of the fermion $f$ and $N_{c}^{f}$ is the color multiplicity ($3$ for
quarks, and $1$ for leptons).

## References

  * (1) S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * (2) F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * (3) H. Georgi, D. B. Kaplan and L. Randall, _Manifesting the Invisible Axion at Low-energies_ , _Phys. Lett. B_ 169 (1986) 73.
  * (4) E. Izaguirre, T. Lin and B. Shuve, _Searching for Axionlike Particles in Flavor-Changing Neutral Current Processes_ , _Phys. Rev. Lett._ 118 (2017) 111802 [1611.09355].
  * (5) M. B. Gavela, R. Houtz, P. Quilez, R. Del Rey and O. Sumensari, _Flavor constraints on electroweak ALP couplings_ , _Eur. Phys. J. C_ 79 (2019) 369 [1901.02031].
  * (6) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _The Low-Energy Effective Theory of Axions and ALPs_ , _JHEP_ 04 (2021) 063 [2012.12272].
  * (7) M. Chala, G. Guedes, M. Ramos and J. Santiago, _Running in the ALPs_ , _Eur. Phys. J. C_ 81 (2021) 181 [2012.09017].
  * (8) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Consistent Treatment of Axions in the Weak Chiral Lagrangian_ , _Phys. Rev. Lett._ 127 (2021) 081803 [2102.13112].
  * (9) M. Bauer, M. Neubert, S. Renner, M. Schnubel and A. Thamm, _Flavor probes of axion-like particles_ , 2110.10698.
  * (10) J. Bonilla, I. Brivio, M. B. Gavela and V. Sanz, _One-loop corrections to ALP couplings_ , _JHEP_ 11 (2021) 168 [2107.11392].
  * (11) A. W. M. Guerrera and S. Rigolin, _Revisiting $K\rightarrow\pi a$ decays_, 2106.05910.
  * (12) R. D. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * (13) L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, _The landscape of QCD axion models_ , _Phys. Rept._ 870 (2020) 1 [2003.01100].
  * (14) Belle-II collaboration, F. Abudinén et al., _Search for Axion-Like Particles produced in $e^{+}e^{-}$ collisions at Belle II_, _Phys. Rev. Lett._ 125 (2020) 161806 [2007.13071].
  * (15) LHCb collaboration, R. Aaij et al., _Search for hidden-sector bosons in $B^{0}\\!\to K^{*0}\mu^{+}\mu^{-}$ decays_, _Phys. Rev. Lett._ 115 (2015) 161802 [1508.04094].
  * (16) LHCb collaboration, R. Aaij et al., _Search for long-lived scalar particles in $B^{+}\to K^{+}\chi(\mu^{+}\mu^{-})$ decays_, _Phys. Rev. D_ 95 (2017) 071101 [1612.07818].
  * (17) BNL-E949 collaboration, A. V. Artamonov et al., _Study of the decay $K^{+}\to\pi^{+}\nu\bar{\nu}$ in the momentum region $140<P_{\pi}<199$ MeV/c_, _Phys. Rev. D_ 79 (2009) 092004 [0903.0030].
  * (18) NA62 collaboration, E. Cortina Gil et al., _Measurement of the very rare K ${}^{+}\to{\pi}^{+}\nu\overline{\nu}$ decay_, _JHEP_ 06 (2021) 093 [2103.15389].
  * (19) R. Essig, R. Harnik, J. Kaplan and N. Toro, _Discovering New Light States at Neutrino Experiments_ , _Phys. Rev. D_ 82 (2010) 113008 [1008.0636].
  * (20) CHARM collaboration, F. Bergsma et al., _A Search for Decays of Heavy Neutrinos_ , _Phys. Lett. B_ 128 (1983) 361.
  * (21) J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu et al., _Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump_ , _Phys. Rev. D_ 38 (1988) 3375.
  * (22) M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer and K. Schmidt-Hoberg, _Revised constraints and Belle II sensitivity for visible and invisible axion-like particles_ , _JHEP_ 12 (2017) 094 [1709.00009].
  * (23) P. Agrawal et al., _Feebly-interacting particles: FIPs 2020 workshop report_ , _Eur. Phys. J. C_ 81 (2021) 1015 [2102.12143].
  * (24) MicroBooNE collaboration, P. Abratenko et al., _Search for a Higgs Portal Scalar Decaying to Electron-Positron Pairs in the MicroBooNE Detector_ , _Phys. Rev. Lett._ 127 (2021) 151803 [2106.00568].
  * (25) V. Cirigliano, G. Ecker, H. Neufeld, A. Pich and J. Portoles, _Kaon Decays in the Standard Model_ , _Rev. Mod. Phys._ 84 (2012) 399 [1107.6001].
  * (26) Y. Aoki et al., _FLAG Review 2021_ , 2111.09849.
  * (27) M. Bauer, M. Neubert and A. Thamm, _Collider Probes of Axion-Like Particles_ , _JHEP_ 12 (2017) 044 [1708.00443].
  * (28) B. Döbrich, J. Jaeckel and T. Spadaro, _Light in the beam dump - ALP production from decay photons in proton beam-dumps_ , _JHEP_ 05 (2019) 213 [1904.02091].
  * (29) J. Jaeckel and M. Spannowsky, _Probing MeV to 90 GeV axion-like particles with LEP and LHC_ , _Phys. Lett. B_ 753 (2016) 482 [1509.00476].
  * (30) E. Goudzovski et al., _New Physics Searches at Kaon and Hyperon Factories_ , 2201.07805.
  * (31) CHARM collaboration, F. Bergsma et al., _Search for Axion Like Particle Production in 400-GeV Proton - Copper Interactions_ , _Phys. Lett. B_ 157 (1985) 458.
  * (32) B. Döbrich, F. Ertas, F. Kahlhoefer and T. Spadaro, _Model-independent bounds on light pseudoscalars from rare B-meson decays_ , _Phys. Lett. B_ 790 (2019) 537 [1810.11336].
  * (33) NA62 collaboration, E. Cortina Gil et al., _Search for a feebly interacting particle $X$ in the decay $K^{+}\rightarrow\pi^{+}X$_, _JHEP_ 03 (2021) 058 [2011.11329].
  * (34) NA62 collaboration, E. Cortina Gil et al., _Search for $\pi^{0}$ decays to invisible particles_, _JHEP_ 02 (2021) 201 [2010.07644].
  * (35) Belle collaboration, J. Grygier et al., _Search for $\boldsymbol{B\to h\nu\bar{\nu}}$ decays with semileptonic tagging at Belle_, _Phys. Rev. D_ 96 (2017) 091101 [1702.03224].
  * (36) Belle collaboration, J. Brodzicka et al., _Physics Achievements from the Belle Experiment_ , _PTEP_ 2012 (2012) 04D001 [1212.5342].
  * (37) BaBar collaboration, B. Aubert et al., _Search for Invisible Decays of a Light Scalar in Radiative Transitions $\upsilon_{3S}\to\gamma$ A0_, in _34th International Conference on High Energy Physics_ , 7, 2008, 0808.0017.
  * (38) BaBar collaboration, J. P. Lees et al., _Search for Invisible Decays of a Dark Photon Produced in ${e}^{+}{e}^{-}$ Collisions at BaBar_, _Phys. Rev. Lett._ 119 (2017) 131804 [1702.03327].
  * (39) L3 collaboration, M. Acciarri et al., _Search for new physics in energetic single photon production in $e^{+}e^{-}$ annihilation at the $Z$ resonance_, _Phys. Lett. B_ 412 (1997) 201.
  * (40) E. Izaguirre, G. Krnjaic, P. Schuster and N. Toro, _New Electron Beam-Dump Experiments to Search for MeV to few-GeV Dark Matter_ , _Phys. Rev. D_ 88 (2013) 114015 [1307.6554].
  * (41) R. Essig, J. Mardon, M. Papucci, T. Volansky and Y.-M. Zhong, _Constraining Light Dark Matter with Low-Energy $e^{+}e^{-}$ Colliders_, _JHEP_ 11 (2013) 167 [1309.5084].
  * (42) G. G. Raffelt, _Astrophysical methods to constrain axions and other novel particle phenomena_ , _Phys. Rept._ 198 (1990) 1.
  * (43) G. Raffelt and D. Seckel, _Bounds on Exotic Particle Interactions from SN 1987a_ , _Phys. Rev. Lett._ 60 (1988) 1793.
  * (44) J. Jaeckel, P. C. Malta and J. Redondo, _Decay photons from the axionlike particles burst of type II supernovae_ , _Phys. Rev. D_ 98 (2018) 055032 [1702.02964].
  * (45) A. Caputo, H.-T. Janka, G. Raffelt and E. Vitagliano, _Low-Energy Supernovae Severely Constrain Radiative Particle Decays_ , 2201.09890.
  * (46) J. H. Chang, R. Essig and S. D. McDermott, _Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle_ , _JHEP_ 09 (2018) 051 [1803.00993].
  * (47) D. Croon, G. Elor, R. K. Leane and S. D. McDermott, _Supernova Muons: New Constraints on $Z$’ Bosons, Axions and ALPs_, _JHEP_ 01 (2021) 107 [2006.13942].
  * (48) P. Coloma, B. A. Dobrescu, C. Frugiuele and R. Harnik, _Dark matter beams at LBNF_ , _JHEP_ 04 (2016) 047 [1512.03852].