# Probing Long-lived Axions at the KOTO Experiment

Yoav Afik<EMAIL_ADDRESS>Experimental Physics Department, CERN, 1211
Geneva, Switzerland  Babette Döbrich<EMAIL_ADDRESS>Max-Planck-
Institut für Physik (Werner-Heisenberg-Institut), Föhringer Ring 6, 80805
München, Germany  Jan Jerhot<EMAIL_ADDRESS>Centre for Cosmology,
Particle Physics and Phenomenology (CP3), Université catholique de Louvain,
Chemin du Cyclotron 2, B-1348 Louvain-la-Neuve, Belgium  Yotam Soreq
<EMAIL_ADDRESS>Physics Department, Technion–Institute of
Technology, Haifa 3200003, Israel  Kohsaku Tobioka<EMAIL_ADDRESS>Department of Physics, Florida State University, Tallahassee, FL 32306, USA
High Energy Accelerator Research Organization (KEK), Tsukuba 305-0801, Japan

(August 28, 2024)

###### Abstract

While the main goal of the J-PARC KOTO experiment is to measure the rare decay
$K_{L}\to\pi^{0}\nu\bar{\nu}$, the unique setup of KOTO raises the possibility
to search for physics beyond the Standard Model, in an attempt to probe parts
of the parameter space which are not covered by other experiments. In this
paper, we test the possibility of using KOTO to search for heavy QCD axions,
or axion-like particles, a well-motivated extension of the Standard Model
emerging in a variety of models. In particular, we estimate the sensitivity of
the current KOTO setup as well as the KOTO Step-2 for various benchmark
scenarios of axion coupling to the Standard Model. We find that KOTO Step-2
can probe new regions in the parameter space, while KOTO with its current form
can only reaffirm the existing bounds. The obtained axion datasets are
available as an update of the public code of the Alpinist framework, including
implementation of KOTO setups in the simulation, allowing for interpretation
of various analyses as searches for axions in custom models.

††preprint: IRMP-CP3-23-10, MPP-2023-40, KEK-TH-2499

##  I Introduction

Rare Kaon decays are a well-known test of the Standard Model (SM) and serve as
a very sensitive probe of New Physics (NP). Two golden channels are
$K^{+}\to\pi^{+}\nu\bar{\nu}$ and $K_{L}\to\pi^{0}\nu\bar{\nu}$, which are
very rare decays with branching ratio (BR) at $\sim 10^{-11}$ level. The NA62
[1] and KOTO [2] experiments aim at measuring these BRs for the first time,
with the latest results given in [3, 4]. Besides testing the SM, including the
Grossman-Nir bound [5], these measurements probe feebly interacting particles
(FIPs) which contribute to the $K\to\pi+\,{\rm invisible}$ decay, see e.g. the
recent review [6].

Both KOTO and NA62 are based on proton-fixed targets, with $30\,{\rm GeV}$ and
$400\,{\rm GeV}$ beams, respectively, and far detection systems. They can
effectively serve as beam-dump experiments probing NP without relying on Kaon
decays since NP particles can be produced already in the target. This was
pointed out in Ref. [7] in the context of NA62 by using a special running mode
of the experiment, see also [8]. This beam-dump potential was pointed out as a
possible explanation for the three candidate events in the KOTO 2019 data [9].

In this work, we study the potential of the KOTO experiment to serve as a
proton beam-dump for sub-GeV NP searches. Unlike NA62, KOTO can probe long-
lived new particles in a di-gamma final state during its Kaon physics running
without needing a dedicated trigger or run, in a completely parasitic scheme,
see Fig. 1. This limitation is caused by the NA62 trigger requiring the
presence of a Kaon in the standard data-taking. Due to the differences in the
beam energy and other geometrical factors, we expect that KOTO will explore a
different region of the parameter space than NA62 and past proton beam-dump
experiments. In particular, we show that the future run of KOTO can search for
NP particles, which are associated with solutions to the Strong CP problem.

Figure 1: A schematic illustration of the KOTO layout for Step-1 (Step-2). See
text for details.

Our primary benchmark model is the axion, $a$, which is a compelling addition
to the SM because it potentially solves the Strong CP problem [10, 11, 12, 13]
and it can be by itself a dark matter (DM) candidate [14, 15, 16]. A similarly
motivated case is the axion-like-particle (ALP) as a pseudo-Nambu-Goldstone
boson of spontaneously broken global symmetry at a high scale $f_{a}$. It can
also provide a portal to the dark sector [17, 18, 19, 20]. In both cases, we
focus on the SM gauge field interactions, which are given by

$\displaystyle\mathcal{L}\supset$
$\displaystyle\,c_{GG}\frac{\alpha_{s}a}{8\pi
f_{a}}G^{a}_{\mu\nu}\tilde{G}^{a\mu\nu}+c_{BB}\frac{\alpha_{Y}a}{8\pi
f_{a}}B_{\mu\nu}\tilde{B}^{\mu\nu}$
$\displaystyle+c_{WW}\frac{\alpha_{2}a}{8\pi
f_{a}}W_{\mu\nu}\tilde{W}^{\mu\nu}\,,$ (1)

where $c_{GG,BB,WW}$ are dimensionless parameters. The SM gauge field strength
is given by $G^{a}_{\mu\nu}$, $B_{\mu\nu}$ and $W_{\mu\nu}$ for the strong,
hypercharge, and weak interactions, respectively, and
$\alpha_{s}=g_{s}^{2}/(4\pi)$ is the strong gauge coupling and similarly for
$\alpha_{Y}$ and $\alpha_{2}$. We assume the axion/ALP mass is sub-GeV,
heavier than the QCD contribution in the light of heavy QCD axion models [21,
22, 23, 24, 25, 26, 27, 28], which revives the low decay constant from the
long-standing bounds [29, 30, 31, 32, 33, 34]. Furthermore, $f_{a}\lesssim
10\,{\rm TeV}$ is favored by the axion quality problem [23]. Since this
scenario is potentially discovered in laboratories, experimental data has been
reinterpreted, leading to additional constraints. Promising probes based on
the future experiments were also proposed [35, 36, 33, 34, 3, 37, 38, 39, 40,
41, 42, 43]. Hereafter, we collectively refer to the heavy QCD axion and ALPs
as axions, where the mass and the couplings are independent parameters.

In the following, we discuss the KOTO experimental setup and data-taking modes
in Sec. II. In Sec. III, we describe the axion production and decay. The
quantitative impact of this analysis is shown in Sec. IV, where we derive the
bounds from current data and estimate the projection for future data-taking.
We conclude in Sec. V.

##  II KOTO setup and data-taking

In this work, we exploit several past and future setups of the KOTO experiment
while accounting for the available information on the experimental conditions.
These setups fall into two independent categories: one regarding the
experiment layout and one regarding the data-taking mode.

We consider two experimental layouts as follows:

  * •

Step-1: the present 2022 layout, which was also used in the 2015 data-taking
[44].

  * •

Step-2: the proposed setup for the future as described in Ref. [45].

A schematic view of the setups is found in Fig. 1. For both steps we consider
two data-taking modes:

  * •

Kaon mode: the standard mode with a $K_{L}$ beam.

  * •

Beam-dump mode: a special run in a beam-dump mode, which includes a shield
that blocks the beam (beam plug) and different selection cuts.

###  II.1 The experimental setups

In both setups, the experiment uses a primary $30\,{\rm GeV}$ proton beam from
the J-PARC Main Ring. The proton beam impinges on a golden target T1 and
generates a secondary hadronic beam which, besides other particles, consists
of $K_{L}$. In the present setup, Step-1 [44], the experiment axis is under
$16^{\circ}$ angle with respect to the primary proton beam, and its front-end
is located $21\,{\rm m}$ from the T1 target with a set of collimators,
sweeping magnets in between for forming the neutral $K_{L}$ beam, and veto
detectors for upstream background suppression. The CsI calorimeter (ECAL),
located $27\,{\rm m}$ from the T1 target, has a $2\,{\rm m}$ diameter with
$15\,{\rm cm}\times 15\,{\rm cm}$ central hole for the beam. The decay volume
is $2.9\,{\rm m}$ long and precedes the ECAL.

The proposed KOTO Step-2 setup [45] with a higher intensity beam assumes a
$5^{\circ}$ angle between the detector and the primary beam. We assume the
beginning of the decay volume to be $45.75\,{\rm m}$ away from the target. The
calorimeter size is increased to a $3\,{\rm m}$ diameter with a $20\,{\rm
cm}\times 20\,{\rm cm}$ central hole for the beam and it is located $64\,{\rm
m}$ far from the target.

In a special beam-dump mode during the operation with the Step-1 setup [46], a
beam plug was placed to close the $K_{L}$ beamline, however, the sweeping
magnet was not functional. The dataset of this run is smaller than in the Kaon
mode. We do not know if the backgrounds stated in [46] could be further
reduced at the analysis level, but without a functioning sweeping magnet, a
0-background setting with the acquired data seems unlikely to us. With a
functional sweeping system or a dedicated run with optimized magnet sweeping,
one may be optimistic that a small background can be achieved.

The above-described experimental layouts and modes are implemented in the
Alpinist framework [47] together with specific selection conditions, such that
simplified simulations of axions production and decays can be performed for
the interpretation of KOTO sensitivity for axion detection. The updated code
and datasets are publicly available at https://github.com/jjerhot/ALPINIST.
Details of the datasets and selection conditions are given in the following
Sec. II.2.

###  II.2 KOTO data-taking modes and their interpretation

While operation in the beam-dump mode can potentially allow a direct search
for particles beyond the SM in a background-clean environment, the majority of
the data is collected in the Kaon mode, with the main aim to measure the
extremely rare $K_{L}\to\pi^{0}\nu\bar{\nu}$ decay.

####  II.2.1 Kaon mode

In the Kaon mode, searching for particles of different origin than in the
$K_{L}$ decay may prove to be challenging due to the backgrounds originating
from both the beam and from upstream. Nevertheless, since the
$K_{L}\to\pi^{0}\nu\bar{\nu}$ with a consequent $\pi^{0}\to\gamma\gamma$ decay
has the same signature as the axion $a\to\gamma\gamma$ decay, this provides an
opportunity to re-interpret the (non-)observation of the $\pi^{0}\nu\bar{\nu}$
signal to constrain also the axion parameter space, see e.g. [9]. Unlike SM
particles, axions can propagate through the beamline elements and upstream
detectors and they can decay off-axis from the neutral $K_{L}$ beam.
Obviously, the different kinematics of these two options would render a
dedicated analysis for the case of the axion possibly more sensitive than what
we state below.

For the re-interpretation of the $K_{L}\to\pi^{0}\nu\bar{\nu}$ search as an
axion search, we re-analyze the $a\to\gamma\gamma$ decays simulated using the
toy Monte Carlo implemented in the Alpinist framework, assuming the
$K_{L}\to\pi^{0}\nu\bar{\nu}$ selection conditions for Step-1 analysis of 2015
data [44], for the future Step-1 run [48] and for Step-2 [45]. These are
summarized in Tab. 1. In particular, we are implementing cuts on the following
kinematic variables. The photon cluster coordinates in the plane perpendicular
to the $K_{L}$ beam axis, $x_{\gamma_{1,2}}$ and $y_{\gamma_{1,2}}$, and $r$
which is the separation distance between the photon clusters. $R_{\text{COE}}$
is the center of energy deposited distance from the beam, based on the photon
position at the calorimeter ($x_{\gamma_{1,2}}$, $y_{\gamma_{1,2}}$) and the
final photon energies, $E_{\gamma_{1,2}}$. The photon separation angle
projection on the calorimeter plane and the angle between the beam axis and
the $\pi\nu\nu$-hypothesis-reconstructed photon momenta are denoted as
$\theta_{\gamma,{\rm calo}}$ and $\theta_{\gamma,{\rm beam}}$, respectively.
The $z_{\text{vtx}}$ position is calculated assuming that we are
reconstructing an on-axis $\pi^{0}\to\gamma\gamma$ decay, where
$z_{\text{vtx}}=0$ corresponds $z=21\,{\rm m}$ from the T1 target for Step-1
and $44\,{\rm m}$ for Step-2. Finally, we quote $\mathcal{A}_{\text{add}}$,
the selection efficiency of additional shape-related cuts (cluster shape,
pulse shape and shower depth). Since we do not have more detailed information
about $\mathcal{A}_{\text{add}}$, we assume a uniform distribution over the
whole signal region, using the numbers quoted in [44, 45]. As we can see from
Tab. 1, the main difference between the current Step-1 data and the future
planned run is the collected statistics in terms of the number of protons on
target $N_{\text{PoT}}$ and the selection efficiency of the shape-related cut
algorithms which has improved considerably while keeping a good rejection
power for hadronic backgrounds.

In the search for $a\to\gamma\gamma$ signal events, we are not limited to the
$K_{L}\to\pi^{0}\nu\bar{\nu}$ signal region, since the $a\to\gamma\gamma$
events have different kinematics, for details see Sec. III. Therefore, we need
to estimate the expected $N_{\rm exp}$ and observed $N_{\rm obs}$ number of SM
events in the whole $z_{\text{vtx}}$-$p_{T,\pi^{0}}$ plane for the various
KOTO datasets. For Step-1 2015 dataset, we make estimations in all regions in
the $z_{\text{vtx}}$-$p_{T,\pi^{0}}$ plane, using the $N_{\rm exp}$ and
$N_{\rm obs}$ shown in Fig. 3 of [44]. In addition, we estimate the
sensitivity of the region of $p_{T,\pi^{0}}>0.5\,{\rm GeV}$ that is out of the
range of the referential figure. Since the pion transverse momentum is
expected to be smaller than $0.5\,{\rm GeV}$ for most of the known physics
processes [49], we assume that there is no SM background in this region, i.e.
$N_{\rm exp}=0$.

In order to project the KOTO sensitivity for the future Step-1 dataset, we use
the $N_{\rm exp}$ backgrounds in the various $z_{\text{vtx}}$-$p_{T,\pi^{0}}$
regions which were presented in [4] for the $N_{\text{PoT}}=3.05\times
10^{19}$ statistics and rescale these numbers to the expected statistics
$N_{\text{PoT}}=14\times 10^{19}$ [48]. If we exclude the $\pi\nu\bar{\nu}$
and the surrounding region ($z_{\text{vtx},\pi\nu\bar{\nu}}<5.1\,{\rm m}$ and
$p_{T,\pi^{0}}<0.26\,{\rm GeV}$), where the number of background events is
large compared to the expected number of $a\to\gamma\gamma$ events,111The SM
$K_{L}\to\pi^{0}\nu\bar{\nu}$ is considered to be a background for the search
for $a\to\gamma\gamma$ decay. we get $N_{\rm exp}=3.35$ events with
$N_{\text{PoT}}=14\times 10^{19}$.

For KOTO Step-2, we use Ref. [49] to estimate the background rate in the
$1.75\;\mathrm{m}<z_{\text{vtx}}<15\,{\rm m}$ and $p_{T,\pi^{0}}>0.4\,{\rm
GeV}$ regions. We find $N_{\rm exp}\approx 1.38$ events with
$N_{\text{PoT}}=6\times 10^{20}$ statistics (assuming again that there is no
background in the $p_{T,\pi^{0}}>0.5\,{\rm GeV}$ region for
$z_{\text{vtx}}<15\,{\rm m}$).

| Step-1 | Step-2
---|---|---
$\sqrt{x_{\gamma_{1,2}}^{2}+y_{\gamma_{1,2}}^{2}}$ | $<0.85\,{\rm m}$ | $<1.35\,{\rm m}$
$\min(|x_{\gamma_{1,2}}|,|y_{\gamma_{1,2}}|)$ | $>0.15\,{\rm m}$ | $>0.175\,{\rm m}$
$R_{\text{COE}}$ | $>0.2\,{\rm m}$ | $-$
$r$ | $>0.3\,{\rm m}$ | $>0.3\,{\rm m}$
$\theta_{\gamma,{\rm calo}}$ | $<150^{\circ}$ | $<150^{\circ}$
$E_{\gamma_{1}}+E_{\gamma_{2}}$ | $>0.65\,{\rm GeV}$ | $>0.5\,{\rm GeV}$
$E_{\gamma_{1,2}}$ | $\in[0.1,2.0]\,{\rm GeV}$ | $>0.1\,{\rm GeV}$
$E_{\gamma_{1}}/E_{\gamma_{2}}$ | $>0.2$ | $-$
$z_{\text{vtx}}$ | $\in[2.9,6.0]\,{\rm m}$ | $\in[1.75,15]\,{\rm m}$
$\theta_{\gamma,{\rm beam}}$ | $>2.5^{\circ}\,{\rm GeV}$ | $-$
$\mathcal{A}_{\text{add}}$ | 0.52 (0.9) | 0.73
$N_{\text{PoT}}\times 10^{19}$ | $2.2\,(14)$ | $60$
Table 1: The selection conditions for $K_{L}\to\pi^{0}\nu\bar{\nu}$ for Step-1
[44] and Step-2 [45]. The $N_{\text{PoT}}$ and shape-related cut efficiency
(in parenthesis) are for the future Step-1 run [48].

Figure 2:  Expected axion production yield $N_{a}$ at the target as a function
of axion energy, $E_{a}$, and its angle to the incident proton beam for two
cases of axion masses and couplings. The respective angles at which the
detector is placed in Step-1 and Step-2 are shown by solid lines with dashed
lines indicating the part of the distribution of interest given the ECAL
edges.

####  II.2.2 Beam-dump mode

In the beam-dump mode, so far KOTO collected data corresponding to
$N_{\text{PoT}}=2.2\times 10^{17}$ [46], which is about two orders of
magnitude less than in the Kaon mode. For our projection of the beam-dump
mode, we assume that ten times more data will be collected in this mode, i.e.
$N_{\text{PoT}}=2.2\times 10^{18}$, while keeping the background under
control, i.e. we consider a background-free search. We explore this case for
both KOTO Step-1 and KOTO Step-2 layouts assuming
$\mathcal{A}_{\text{add}}=1$. The selection conditions are simply both photons
being in the calorimeter acceptance with cluster distance $>0.3\,{\rm m}$ (as
used in the $\pi^{0}\nu\bar{\nu}$ analysis). For KOTO Step-1 dump-mode we
require at least $50\,{\rm MeV}$ deposited on the calorimeter per photon and
for KOTO Step-2 at least $100\,{\rm MeV}$ per photon and at least $500\,{\rm
MeV}$ in total deposited on the calorimeter.

##  III Axion production and detection

The conventional axion production at the $K_{L}$ experiments is from
$K_{L}\to\pi^{0}a$ decay, which is CP violating and thereby suppressed. Here,
the relevant production to probe the long-lived axion occurs at the fixed
target T1 where $K_{L}$ is produced, i.e. in the proton-gold collisions. We
consider two production mechanisms: Primakoff production and axion-meson
mixing, both are implemented in the Alpinist framework, which uses Pythia 8
[50] to generate meson distributions. In the following, we show the axion
production yields for these mechanisms for the $30\,{\rm GeV}$ proton beam.
The validation of obtained yields with the detector under the $16^{\circ}$
angle is given in App. A. As shown in Ref. [8], the yield and the momentum
spectrum is well described by Pythia 8 also for angles smaller than
$16^{\circ}$ and higher beam energy.

###  III.1 Axion production in the target

The gluon coupling $c_{GG}$ induces the axion mixing to the neutral mesons of
the same quantum numbers, $P\in\left\\{\pi^{0},\eta,\eta^{\prime}\right\\}$.
The axion yield from mixing production is then approximately given by

$\displaystyle N_{a}^{\text{mix}}\approx
N_{\pi^{0}}\cdot|\theta_{a\pi}|^{2}+N_{\eta}\cdot|\theta_{a\eta}|^{2}+N_{\eta^{\prime}}\cdot|\theta_{a\eta^{\prime}}|^{2}$
(2)

where $N_{\pi^{0}}$ is the production yield of neutral pion, $\theta_{a\pi}$
is the pion-axion mixing angle, and similar notations are applied to $\eta$
and $\eta^{\prime}$. To leading order in $f_{\pi}/f_{a}$, the mixing angles
$\theta_{aP}$ are

$\displaystyle\theta_{aP}\approx\frac{f_{\pi}}{f_{a}}\frac{K_{aP}m_{a}^{2}+m_{aP}^{2}}{m_{a}^{2}-m_{P}^{2}},$
(3)

where $f_{\pi}\approx 93\,{\rm MeV}$ and we use the same notation for the
kinetic- and mass-mixing ($K_{aP}$ and $m_{aP}$) as in [38, 47]. The
$\eta$–$\eta^{\prime}$ mixing angle used is
$\sin\theta_{\eta\eta^{\prime}}=-1/3$. The kinematics for the axion mixing
production is treated according to App. D of [47]. The Primakoff process is
governed by the axion-photon interaction

$\displaystyle c_{\gamma\gamma}\frac{\alpha_{\rm EM}a}{8\pi
f_{a}}F_{\mu\nu}\tilde{F}^{\mu\nu},$ (4)

where $\alpha_{\rm EM}$ is the fine-structure constant and $c_{\gamma\gamma}$
is the effective axion-photon coupling. The axion-photon coupling is generated
by $c_{BB}$ and $c_{WW}$ couplings and for $m_{a}\lesssim m_{\rho}$ also by
the low-energy contribution of $c_{GG}$ [51, 37, 52]:

$\displaystyle c_{\gamma\gamma}=$ $\displaystyle
c_{BB}+c_{WW}-c_{GG}\left(1.92+2\sum_{P}\frac{f_{a}}{f_{P}}\theta_{aP}\right),$
(5)

where $f_{\eta^{\prime}}\approx 73\,{\rm MeV}$. There are two sources of
photons in the target that can produce axions in the interaction with the
target nuclei via the Primakoff process: off-shell photons from the proton of
the primary beam [7] and photons from decays of secondary neutral
pseudoscalars produced in the target. The latter has been shown to be
dominating [8]. The axion yield from these processes is given by their sum and
is denoted as $N_{a}^{\text{Prim}}$.

Additional production could be from flavor-changing Kaon decays near the fixed
target. With the gluon coupling $c_{GG}$, the
$K^{\pm}(K_{S})\to\pi^{\pm}(\pi^{0})a$ is not suppressed by a loop or CP
violation [53, 43, 6]. The $c_{WW}$ coupling induces the same processes at one
loop [54]. These production rates could be sizable because the total width of
Kaons is small, i.e., the BR is enhanced. However, most Kaons are removed by
the collimators or deflected by the magnetic field. Including these effects,
we estimate that $K^{+}\to\pi^{+}a$ is subdominant. Still, the production from
$K_{S}$ is potentially interesting due to the shorter lifetime, but it
requires a simulation of $K_{S}$ transport, which is beyond the scope of this
paper. Therefore, we neglect axion production from $K^{+}$ and $K_{S}$ decays.

To summarize, the total axion yield at the T1 target is approximated to be

$\displaystyle N_{a}\approx N_{a}^{\text{mix}}+N_{a}^{\text{Prim}}\,,$ (6)

where the exact relative contributions of the two mechanisms depend on the
specific model. We show the differential distributions with respect to the
axion energy ($E_{a}$) and its production angle to the proton beam
($\theta_{a}$, see Fig. 1) in Fig. 2 for two benchmarks with axion mass of 40
MeV or 400 MeV and the coupling being dominated by $c_{GG}$. In this coupling
benchmark for the $m_{a}=400$ MeV case, the overall production yield from the
mixing with $\pi^{0}$ and $\eta$ exceeds significantly the Primakoff
production while for smaller axion masses, the Primakoff production becomes
more relevant, similarly to what has been observed in [47] for experiments
operating with higher beam energies.

Figure 3:  A schematic picture of the axion event in the Kaon mode. The disk
represents the ECAL, and the axis is the $K_{L}$ beamline. The distance from
the ECAL to the upstream is $dz$. A typical axion trajectory is almost
parallel to the beam axis in the distance $\rho_{a}$ (see Fig. 1). If the
final state photons from the axion decay at $dz_{a}$ leave the energy of
$E_{1,2}$ with a separation of $D_{12}$ on the ECAL, the vertex position $dz$
is reconstructed on the beam axis.

###  III.2 Axion detection mechanism

Given the beam energy and the distance between the axion production location
and the detector, the KOTO experiment can search for long-lived axions. Below,
we see that for the relevant masses, the dominant decay channel is
$a\to\gamma\gamma$, which is given by

$\displaystyle\Gamma_{\gamma\gamma}=\frac{\alpha_{\rm
EM}^{2}m_{a}^{3}}{256\pi^{3}}\frac{c_{\gamma\gamma}^{2}}{f_{a}^{2}}\,,$ (7)

where the effective photon coupling $c_{\gamma\gamma}$ is defined in Eq. (5).
For an axion heavier than 1 GeV with non-zero $c_{GG}$ coupling, the width of
hadronic decay modes, such as $a\to\pi\pi\eta$, estimated in [38], dominate
the total width. However, the relevant final state for search at KOTO is di-
photon, which typically becomes a sub-dominant mode for $m_{a}\gtrsim
0.5\,{\rm GeV}$ (with large model dependence), resulting in reduced
sensitivity. Details on KOTO sensitivity for hadronic axion decays can be
found in App. B.

Long-lived axions produced at the fixed target can reach the distant decay
volume, and a di-photon decay leaves a characteristic signal, see a schematic
picture of the axion decay event in Fig. 3. When axions enter the decay
volume, they are almost parallel to the $K_{L}$ beam axis, because the
distance to the ECAL is larger than the ECAL size, but away from the axis with
the distance $\rho_{a}$. The two photons from a $a\to\gamma\gamma$ decay in
the decay volume can then hit the ECAL, mimicking the signal of
$K_{L}\to\pi^{0}\nu\bar{\nu}$. If the standard reconstruction algorithm for
$K_{L}\to\pi^{0}\nu\bar{\nu}$ is applied, the reconstructed position ($dz$)
will be different from the true distance between the axion decay point and the
ECAL ($dz_{a}$), but the event is not discarded. In this sense, this signal is
similar to the halo $K_{L}\to\gamma\gamma$ background that the KOTO
collaboration found in the earlier data, but the reconstructed distributions
are typically different.

Contrary to the Kaon mode, the detection mechanism in the beam-dump mode is
relatively straightforward as it is a dedicated run to search for long-lived
particles decaying to photons.

Figure 4:  Axion signal event distributions $N_{\text{sig}}$ as could be found
in the $\pi\nu\bar{\nu}$ analysis at KOTO Step-2 for specific axion models.
For the analytic understanding of the distributions, we show the lines
corresponding to the Eq. (11) with $\rho_{a}=0.175,1.35$ m and a specific
$E_{a}$ as dashed magenta lines.

###  III.3 Signal from axion decay in Kaon mode

The long-lived axion $a\to\gamma\gamma$ decay passes the event selections
described in Sec. II.2.1 when $\rho_{a}>0$ which introduces factitious
transverse momentum of the di-photon system. Some selection criteria are
universal, often limited by the experimental resolution, but the remaining
cuts assume the topology of $K_{L}\to\pi^{0}(\to\gamma\gamma)\nu\bar{\nu}$.
The KOTO experiment could implement a dedicated analysis for the long-lived
particles, but here we simply adopt the standard $K_{L}\to\pi^{0}\nu\bar{\nu}$
analysis because the backgrounds are well-investigated. Assuming
$K_{L}\to\pi^{0}\nu\bar{\nu}$ topology for the long-lived axion decays leads
to several non-trivial characteristics of the axion signal event distribution
$N_{\text{sig}}$ in the $z_{\text{vtx}}$-$p_{T,\pi^{0}}$ plane. Example
distributions for parameters for which we expect interesting sensitivities of
Step-2 are shown in Fig. 4. In the following, we give analytic understandings
of the characteristics based on several simplifications.

We assume that the axion enters the decay volume in parallel to the beamline
with distance $\rho_{a}$. The axion invariant mass is given by the photon
energies, $E_{1,2}$, and the opening angle, $\theta_{12}$,

$\displaystyle m_{a}^{2}=2E_{1}E_{2}(1-\cos\theta_{12})\simeq
E_{1}E_{2}\theta_{12}^{2}\,.$ (8)

Using this, the separation of the two photons on the ECAL, $D_{12}$, is
approximately

$\displaystyle D_{12}\simeq dz_{a}\theta_{12}\simeq
dz_{a}\frac{m_{a}}{\sqrt{E_{1}E_{2}}}\,.$ (9)

Then, the reconstructed vertex position, assuming the event topology of
$K_{L}\to\pi^{0}\nu\bar{\nu}$, is given by

$\displaystyle dz\equiv D_{12}\frac{\sqrt{E_{1}E_{2}}}{m_{\pi^{0}}}\simeq
dz_{a}\frac{m_{a}}{m_{\pi^{0}}}\,.$ (10)

The distance from the ECAL ($dz$) is translated to the standard coordinate
system by $z_{\text{vtx}}=L-dz$ where $L=$6 (20) m in Step-1 (Step-2). Because
only the axion decays in the decay volume are accepted, there exists a
limitation of $dz_{a}<dz_{a}^{\rm max}=2.9\,(18.25)\,{\rm m}$ in Step-1
(Step-2). This limitation together with Eq. (10) gives us a condition on the
maximum spread of the distribution over the reconstructed vertex position as
$z_{\text{vtx}}>L-dz_{a}^{\rm max}({m_{a}}/{m_{\pi^{0}}})$. The boundary
$z_{\text{vtx}}^{\text{min}}$ is shown as a blue vertical line of Fig. 4 left.

Another feature can be seen as a correlation with both $p_{T}^{\pi^{0}}$ and
$z_{\text{vtx}}$. In the case of $K_{L}\to\pi^{0}\nu\bar{\nu}$, the transverse
kick is from $K_{L}$ decay, and hence, $p_{T}^{\pi^{0}}<m_{K_{L}}$. However,
the transverse asymmetry of the ECAL hits is merely from the transverse
position of the incident axion. Suppose the distance from the axion to the
beamline is $\rho_{a}$, the reconstructed $p_{T}^{\pi^{0}}$ is roughly

$\displaystyle
p_{T}^{\pi^{0}}\simeq\frac{E_{a}\rho_{a}}{\sqrt{\rho_{a}^{2}+dz^{2}}}\simeq\frac{E_{a}}{\sqrt{1+(L-z_{\text{vtx}})^{2}/\rho_{a}^{2}}}.$
(11)

where $E_{a}=E_{1}+E_{2}$ is the axion energy. The correlation between
$p_{T}^{\pi^{0}}$ and $z_{\text{vtx}}$ is explained by the above formula. Note
that $p_{T}^{\pi^{0}}$ can easily exceed $m_{K_{L}}$ because
$E_{a}\sim\mathcal{O}(1)\,{\rm GeV}$.

Figure 5: 90% CL exclusion bounds (filled contours) and projected limits
(empty contours) for scenarios (i)-(iii) for all KOTO setups considered
compared to the exclusions from beam-dump experiments, exclusion from $B\to
Ka$ and $K\to\pi a$ decays and the bound from the supernova SN1987A (shown as
a dashed line as it is affected by significant uncertainties, see e.g. [51,
55]).

The last feature that can be observed is the vanishing of the distribution at
large $p_{T}^{\pi^{0}}$ and $z_{\text{vtx}}$, which is a consequence of the
two-photon separation cut, $D_{12}\leq 0.3$ m. This can be understood by
combining Eqs. (9), (10) and (11) with a simplification of $E_{1,2}\approx
E_{a}/2$ as

$\displaystyle
p_{T}^{\pi^{0}}\lesssim\frac{2m_{\pi^{0}}}{0.3\;\mathrm{m}\sqrt{(L-z_{\text{vtx}})^{-2}+\rho_{a}^{-2}}}\,.$
(12)

Therefore, the distribution between the two dashed lines from Eq. (11)
vanishes at large $p_{T}^{\pi^{0}}$ and $z_{\text{vtx}}$ at an approximate
bound corresponding to Eq. (12) with $\rho_{a}=1.35\,{\rm m}$.

Since $p_{T}^{\pi^{0}}$ of the long-lived axion events can be significantly
larger compared to the $K_{L}$ events or background, we use high
$p_{T}^{\pi^{0}}$ region. The detail of the signal region and the expected
background yield is discussed in Sec. II.2.1.

##  IV Bounds and projections for axions

Figure 6:  90% CL exclusion bounds and projected limits for fixed axion mass
and variable $c_{BB}$ and $c_{GG}$ couplings.

We derive the bounds for different axion models using the current KOTO Step-1
results and project the sensitivities of the future runs. For this purpose, we
simulate the signal, including all production processes, i.e. Primakoff and
meson mixing mechanisms. The various experimental setups are implemented in
the Alpinist framework where we have also generated the datasets222For
increasing the precision of the estimated number of observed events, each mass
and decay width bin is evaluated with 2 million axion decay events. with the
number of events expected $N_{\text{sig}}$ for each setup and production mode.
These tables can be further used for showing the $N_{\text{sig}}$
distributions for any setup of model-dependent parameters using the re-scaling
module of the framework. In Fig. 5, we show the current bounds and projections
for the expected sensitivity with future data corresponding to the following
three benchmark models:

  1. (i)

Heavy QCD axion, $c_{GG}\neq 0$ with $c_{BB}=c_{WW}=0$;

  2. (ii)

Hypercharge dominant, $c_{BB}\neq 0$ with $c_{WW}=c_{GG}=0$;

  3. (iii)

Codominance, $c_{BB}=c_{WW}=c_{GG}$;

and in Fig. 6, we show the current and projected results for several fixed
masses for variable $c_{GG}$ vs $c_{BB}$ couplings, which is expected to be
similar to the case of $c_{GG}$ vs $c_{WW}$, up to the FCNC production.

We compare our KOTO bounds and projections to the existing bounds from
different experiments. We consider electron beam-dumps E137 [35] and E141
[56], where we use the Alpinist framework for the interpretation of the data
provided in [57], implementation that has been already done for [58], and we
perform a dedicated analysis to derive the proton beam-dump bounds based on
CHARM [33] and NuCal [36], as done in [47] (bounds in gray shade).

The axion can be produced by flavor-changing meson decays, especially in the
presence of $c_{WW}$ and $c_{GG}$ (the corresponding bounds are shown in the
blue shade). We adopt the results of $K^{+}\to\pi^{+}a$ from Ref. [43] and
apply the bounds of E949 [59, 34] and NA62 [60, 3]. The scheme of the recast
is found in [61]. For $B\to Ka(\to\gamma\gamma)$, the BaBar bound can be used
[62, 42]. The two-loop production calculation with $c_{GG}$ is found in [41].
The contribution from $c_{WW}$ at one-loop is calculated in [54], but it is
numerically subdominant for the benchmark (iii), so for simplicity, we
approximate $B\to Ka$ by the $c_{GG}$ contribution. The total width $K^{+}$
would be modified significantly for the low $f_{a}/c_{GG}$, which leads to a
relevant bound at $m_{a}\sim m_{\pi^{0}}$. Requiring ${\rm
BR}(K^{+}\to\pi^{+}a)<3\times 10^{-3}$ based on Sec. 2.2.2 of [6] results in
the bound $f_{a}/c_{GG}\lesssim 5$ GeV. The $c_{BB}$ only scenario is not
significantly constrained by the meson decays since the production originates
from electroweak two-loop diagrams. Therefore, in Fig. 6 the meson decay
bounds are omitted because the corresponding bounds in the limit of the
benchmark (ii) are unknown. Finally, the shown SN1987A bounds are those
derived in [51] but are plotted with a dashed line as their robustness is
under debate [55].

We find that at Step-1, KOTO cannot probe new regions in the parameter space
and it is sensitive only in regions already covered by other proton beam-dump
experiments for all coupling scenarios. The non-observation of additional
signal on top of the expected backgrounds in the past KOTO
$K_{L}\to\pi^{0}\bar{\nu}\nu$ analyses only confirms the results of these past
experiments. While for scenarios with photon coupling domination in the future
also KOTO Step-2 cannot compete with the past electron beam-dump experiments
E137 and E141, it can probe new regions of parameter space for larger masses
($m_{a}\gtrsim m_{\pi^{0}}$) for scenarios with gluonic coupling thanks to
enhanced axion production through mixing with other neutral pseudoscalars.

##  V Conclusions

In this paper, we show that the KOTO experiment, beyond its conventional
purposes, can perform long-lived particle searches in its two different data-
taking modes, the Kaon and the beam-dump. In both modes, NP particles are
produced at the proton target interaction point and can decay in the whole
decay volume of the detector (see Figs. 1 and 3). We show that the future KOTO
runs will explore uncharted parameter space of sub-GeV axions which may
address the Strong CP problem.

Firstly, we show that the Kaon mode, where the majority of the KOTO data is
taken, is sensitive not only to $K_{L}\to\pi^{0}\nu\bar{\nu}$ and
$K_{L}\to\pi^{0}a$ but also to axions which are originating from the
interaction in the proton target and mimic the rare Kaon decay signal. The
main difference between the two signals is that the axion events extend the
distribution of $p_{T}^{\pi^{0}}$ greater than $m_{K_{L}}$, as shown in Fig.
4. This region is unphysical for the di-photon events from the Kaon decays,
thus, we assume no SM background there. Even though our derived constraints
based on KOTO 2015 dataset are currently not competitive, they reaffirm the
constraints obtained with experiments of very different topology and proton
impact energy. We have also shown in Fig. 5 that KOTO in Step-2 can indeed
explore new parameter space for axions with $c_{GG}$ coupling and
$m_{a}\gtrsim 100\,{\rm MeV}$, without changes of the main analysis steps.

Secondly, we have evaluated projections of KOTO running in the beam-dump mode
as recently presented by the collaboration [46]. Here we have shown that KOTO,
due to its low proton beam energy and large angle between the beam and the
detector can especially well explore parameter space at very low couplings,
complementary to such searches at higher energies, e.g. NA62 [7], FASER [63]
or DarkQuest [64]. Although the analysis in beam-dump mode is suitable for
searches for long-lived particles, the sensitivity to the axions is weaker
than in the Kaon mode because the expected statistics is significantly lower
considering the KOTO physics goals.

In this study, we have focused on the axions to demonstrate the proof of
concept, whereas similar analyses could be performed to update the bounds of
other long-lived particles from past proton beam experiments. Furthermore, a
dedicated analysis for long-lived particles rather than reinterpretation of
the $K_{L}\to\pi^{0}\nu\bar{\nu}$ analysis could further improve the
sensitivity although it requires additional background studies. This work also
updates the Alpinist framework [65] to include the KOTO geometry and the
kinematics of the various processes. Thereby, other scenarios, including
axions with different parameter combinations, can be easily studied.

###### Acknowledgements.

BD acknowledges funding through the European Research Council under grant
ERC-2018-StG-802836 (AxScale project) and the Lise Meitner program of the Max
Planck society. JJ acknowledges funding by the F.R.S.-FNRS, Belgium, through
grant FRIA/FC-36305. YS is supported by grants from NSF-BSF (No. 2021800), ISF
(No. 482/20), BSF (No. 2020300) and the Azrieli foundation. KT acknowledges
funding through the US Department of Energy grant DE-SC0010102 and Japan
Society for the Promotion of Science KAKENHI No. 21H01086.

##  Appendix A Validation of the Simulation

Figure 7:  Comparison of the $K_{L}$ total momentum distribution from the
simulation used in this work (Pythia 8) and the data derived by KOTO [66]
(KOTO Data). The red curve presents a fit of the data done by KOTO (KOTO Fit).

Figure 8:  Expected number of signal events $N_{\text{sig}}$ for given values
of axion mass $m_{a}$ and coupling $c_{GG}$ for scenario (i) with
$N_{\text{PoT}}=2.2\times 10^{17}$ in KOTO Step-2 in beam-dump mode.
Comparison between $N_{\text{sig}}$ from $a\to\pi^{0}\pi^{0}\eta$ and $a\to
3\pi^{0}$ decays (left) with $a\to\gamma\gamma$ decay (right).

For the estimation of the axion flux from both the Primakoff and the mixing
production, we first need to estimate the flux of the $\pi^{0}$, $\eta$ and
$\eta^{\prime}$ mesons. For this purpose, we use Pythia 8 [50] for the
generation of the $pp$ interactions. Since we cannot directly validate the
meson flux simulation, we use the measured $K_{L}$ flux at KOTO [66] to
normalize the meson yields. The measured number of $K_{L}$ collimated into
$8\times 8\,{\rm cm}^{2}$ [44] at the end of a $20\,{\rm m}$ long beamline is
$(4.2\pm 0.02_{\text{stat}}\pm 0.06_{\text{sys}})\times 10^{7}$ $K_{L}$ per
$2\times 10^{14}$ protons on target, which corresponds to measured number of
$K_{L}$ per proton on target $N_{K_{L}}^{\text{KOTO}}\sim 2.1\times 10^{-7}$.
Therefore, we normalize the number of mesons simulated by Pythia 8 to the
measured number as:

$\displaystyle N_{P}=\frac{N_{P}^{\text{sim}}}{N_{K_{L}}^{\text{sim}}}\times
N_{K_{L}}^{\text{KOTO}},$ (13)

where $N_{P}^{\text{sim}}$ is the number of ($\pi^{0}$, $\eta$,
$\eta^{\prime}$) mesons produced per $pp$ interaction in the simulation. For
$10^{8}$ simulated $pp$ interactions, after accounting for $K_{L}$ decays
assuming the peak $K_{L}$ momentum $1.4\,{\rm GeV}$ [44] and $60\%$ loss of
$K_{L}$ due to absorption in the beamline material [66], we obtain the
simulated number of $K_{L}$ per $pp$ interaction $N_{K_{L}}^{\text{sim}}\sim
2.0\times 10^{-7}$. Therefore
$N_{K_{L}}^{\text{KOTO}}/N_{K_{L}}^{\text{sim}}\sim 1.05$, showing a good
agreement between the total number of $K_{L}$ measured and simulated using
Pythia 8. For the multiplicities of ($\pi^{0}$, $\eta$, $\eta^{\prime}$)
mesons, based on the simulation, we obtain:

$\displaystyle\frac{N_{\pi^{0}}^{\text{sim}}}{N_{K_{L}}^{\text{sim}}}\sim
21\,,\quad\frac{N_{\eta}^{\text{sim}}}{N_{K_{L}}^{\text{sim}}}\sim
2.2\,,\quad\frac{N_{\eta^{\prime}}^{\text{sim}}}{N_{K_{L}}^{\text{sim}}}\sim
0.17\,.$ (14)

Furthermore, since axions interact rarely with ordinary matter, the surface of
the axion flux potentially entering the detector
$S_{\text{axion}}^{\text{sim}}$ is much larger, occupying the whole decay
volume plane. When compared to the $K_{L}$ flux, which is collimated to the
$S_{K_{L}}^{\text{KOTO}}=8\times 8\,{\rm cm}^{2}$ profile at the end of the
beamline, the ratio of the two surfaces is about
$S_{\text{axion}}^{\text{sim}}/S_{K_{L}}^{\text{KOTO}}\sim 490$.

Finally, in order to validate the kinematic distributions, we compare the
distributions of $K_{L}$ obtained with Pythia 8 with the distributions
measured by KOTO [66]. As shown in Fig. 7, a good agreement is observed
between the shapes of the distribution of the $K_{L}$ total momentum from the
simulation used in this paper and the data measured by KOTO. This validation
for $K_{L}$ gives some credibility also to the simulated ($\pi^{0}$, $\eta$,
$\eta^{\prime}$) distributions at $30\,{\rm GeV}$ and therefore to the
validity of the expected axion distributions which are used in this work. The
obtained axion distributions are publicly available at [65].

##  Appendix B Sensitivity for hadronic axion decays

KOTO could be in principle sensitive to decays $a\to\pi^{0}\pi^{0}\eta$ or
$a\to 3\pi^{0}$ with a subsequent $\pi^{0}(\eta)\to\gamma\gamma$ decay
resulting in a 6-cluster event. While studying hadronic decays in the Kaon
mode would require a dedicated analysis to address the various backgrounds, in
the case of beam-dump mode, it is stated in [46] that no 6-cluster events have
been found in the collected sample, indicating that the $K_{L}$ background is
kept under control in this case.

After running a simulation with Alpinist for $a\to 3\pi^{0}$ and
$a\to\pi^{0}\pi^{0}\eta$ decays with simple selection criteria on minimal
cluster energy and cluster separation as for the $a\to\gamma\gamma$ decay in
the beam-dump mode, we did not find KOTO sensitivity with hadronic decays to
surpass the sensitivity using $a\to\gamma\gamma$ even for larger axion masses.
As can be seen in Fig. 8 for KOTO Step-2, the difference in the number of
observable signal events $N_{\text{sig}}$ is by several orders of magnitude
smaller for combined search for $a\to 3\pi^{0}$ and $a\to\pi^{0}\pi^{0}\eta$
decays compared to $a\to\gamma\gamma$. Nevertheless, for convenience, we also
provide the resulting datasets for these hadronic decays.

## References

  * [1] Eduardo Cortina Gil et al.  The Beam and detector of the NA62 experiment at CERN.  JINST, 12(05):P05025, 2017, 1703.08501.
  * [2] Taku Yamanaka.  The J-PARC KOTO experiment.  PTEP, 2012:02B006, 2012.
  * [3] Eduardo Cortina Gil et al.  Measurement of the very rare K+→${\pi}^{+}\nu\overline{\nu}$ decay.  JHEP, 06:093, 2021, 2103.15389.
  * [4] J. K. Ahn et al.  Study of the $K_{L}\to\pi^{0}\nu\bar{\nu}$ Decay at the J-PARC KOTO Experiment.  Phys. Rev. Lett., 126(12):121801, 2021, 2012.07571.
  * [5] Yuval Grossman and Yosef Nir.  K(L) —$>$ pi0 neutrino anti-neutrino beyond the standard model.  Phys. Lett. B, 398:163–168, 1997, hep-ph/9701313.
  * [6] Evgueni Goudzovski et al.  New physics searches at kaon and hyperon factories.  Rept. Prog. Phys., 86(1):016201, 2023, 2201.07805.
  * [7] Babette Döbrich, Joerg Jaeckel, Felix Kahlhoefer, Andreas Ringwald, and Kai Schmidt-Hoberg.  ALPtraum: ALP production in proton beam dump experiments.  JHEP, 02:018, 2016, 1512.03069.
  * [8] Babette Döbrich, Joerg Jaeckel, and Tommaso Spadaro.  Light in the beam dump - ALP production from decay photons in proton beam-dumps.  JHEP, 05:213, 2019, 1904.02091.  [Erratum: JHEP 10, 046 (2020)].
  * [9] Teppei Kitahara, Takemichi Okui, Gilad Perez, Yotam Soreq, and Kohsaku Tobioka.  New physics implications of recent search for $K_{L}\to\pi^{0}\nu\bar{\nu}$ at KOTO.  Phys. Rev. Lett., 124(7):071801, 2020, 1909.11111.
  * [10] R. D. Peccei and Helen R. Quinn.  CP Conservation in the Presence of Instantons.  Phys. Rev. Lett., 38:1440–1443, 1977.
  * [11] R. D. Peccei and Helen R. Quinn.  Constraints Imposed by CP Conservation in the Presence of Instantons.  Phys. Rev. D, 16:1791–1797, 1977.
  * [12] Steven Weinberg.  A New Light Boson?  Phys. Rev. Lett., 40:223–226, 1978.
  * [13] Frank Wilczek.  Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons.  Phys. Rev. Lett., 40:279–282, 1978.
  * [14] John Preskill, Mark B. Wise, and Frank Wilczek.  Cosmology of the Invisible Axion.  Phys. Lett. B, 120:127–132, 1983.
  * [15] Michael Dine and Willy Fischler.  The Not So Harmless Axion.  Phys. Lett. B, 120:137–141, 1983.
  * [16] L. F. Abbott and P. Sikivie.  A Cosmological Bound on the Invisible Axion.  Phys. Lett. B, 120:133–136, 1983.
  * [17] Yasunori Nomura and Jesse Thaler.  Dark Matter through the Axion Portal.  Phys. Rev. D, 79:075008, 2009, 0810.5397.
  * [18] Marat Freytsis and Zoltan Ligeti.  On dark matter models with uniquely spin-dependent detection possibilities.  Phys. Rev. D, 83:115009, 2011, 1012.5317.
  * [19] Matthew J. Dolan, Felix Kahlhoefer, Christopher McCabe, and Kai Schmidt-Hoberg.  A taste of dark matter: Flavour constraints on pseudoscalar mediators.  JHEP, 03:171, 2015, 1412.5174.  [Erratum: JHEP 07, 103 (2015)].
  * [20] Yonit Hochberg, Eric Kuflik, Robert Mcgehee, Hitoshi Murayama, and Katelin Schutz.  Strongly interacting massive particles through the axion portal.  Phys. Rev. D, 98(11):115031, 2018, 1806.10139.
  * [21] Hajime Fukuda, Keisuke Harigaya, Masahiro Ibe, and Tsutomu T. Yanagida.  Model of visible QCD axion.  Phys. Rev. D, 92(1):015021, 2015, 1504.06084.
  * [22] Prateek Agrawal, Gustavo Marques-Tavares, and Wei Xue.  Opening up the QCD axion window.  JHEP, 03:049, 2018, 1708.05008.
  * [23] Prateek Agrawal and Kiel Howe.  Factoring the Strong CP Problem.  JHEP, 12:029, 2018, 1710.04213.
  * [24] M. K. Gaillard, M. B. Gavela, R. Houtz, P. Quilez, and R. Del Rey.  Color unified dynamical axion.  Eur. Phys. J. C, 78(11):972, 2018, 1805.06465.
  * [25] Tony Gherghetta, Valentin V. Khoze, Alex Pomarol, and Yuri Shirman.  The Axion Mass from 5D Small Instantons.  JHEP, 03:063, 2020, 2001.05610.
  * [26] R. S. Gupta, V. V. Khoze, and M. Spannowsky.  Small instantons and the strong CP problem in composite Higgs models.  Phys. Rev. D, 104(7):075011, 2021, 2012.00017.
  * [27] Tony Gherghetta and Minh D. Nguyen.  A Composite Higgs with a Heavy Composite Axion.  JHEP, 12:094, 2020, 2007.10875.
  * [28] Alessandro Valenti, Luca Vecchi, and Ling-Xiao Xu.  Grand Color axion.  JHEP, 10:025, 2022, 2206.04077.
  * [29] V. Anastassopoulos et al.  New CAST Limit on the Axion-Photon Interaction.  Nature Phys., 13:584–590, 2017, 1705.02290.
  * [30] Georg G. Raffelt.  Astrophysical axion bounds.  Lect. Notes Phys., 741:51–71, 2008, hep-ph/0611350.
  * [31] G. G. Raffelt.  Stars as laboratories for fundamental physics: The astrophysics of neutrinos, axions, and other weakly interacting particles.  University of Chicago Press, 5 1996.
  * [32] Alexander Friedland, Maurizio Giannotti, and Michael Wise.  Constraining the Axion-Photon Coupling with Massive Stars.  Phys. Rev. Lett., 110(6):061101, 2013, 1210.1271.
  * [33] F. Bergsma et al.  Search for Axion Like Particle Production in 400-GeV Proton - Copper Interactions.  Phys. Lett. B, 157:458–462, 1985.
  * [34] A. V. Artamonov et al.  Search for the decay K+ to pi+ gamma gamma in the pi+ momentum region P $>$ 213 MeV/c.  Phys. Lett. B, 623:192–199, 2005, hep-ex/0505069.
  * [35] J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu, L. W. Mo, T. A. Nunamaker, and P. Rassmann.  Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump.  Phys. Rev. D, 38:3375, 1988.
  * [36] J. Blumlein et al.  Limits on neutral light scalar and pseudoscalar particles in a proton beam dump experiment.  Z. Phys. C, 51:341–350, 1991.
  * [37] Martin Bauer, Matthias Neubert, and Andrea Thamm.  Collider Probes of Axion-Like Particles.  JHEP, 12:044, 2017, 1708.00443.
  * [38] Daniel Aloni, Yotam Soreq, and Mike Williams.  Coupling QCD-Scale Axionlike Particles to Gluons.  Phys. Rev. Lett., 123(3):031803, 2019, 1811.03474.
  * [39] Anson Hook, Soubhik Kumar, Zhen Liu, and Raman Sundrum.  High Quality QCD Axion and the LHC.  Phys. Rev. Lett., 124(22):221801, 2020, 1911.12364.
  * [40] Kevin J. Kelly, Soubhik Kumar, and Zhen Liu.  Heavy axion opportunities at the DUNE near detector.  Phys. Rev. D, 103(9):095002, 2021, 2011.05995.
  * [41] Sabyasachi Chakraborty, Manfred Kraus, Vazha Loladze, Takemichi Okui, and Kohsaku Tobioka.  Heavy QCD axion in b→s transition: Enhanced limits and projections.  Phys. Rev. D, 104(5):055036, 2021, 2102.04474.
  * [42] Emilie Bertholet, Sabyasachi Chakraborty, Vazha Loladze, Takemichi Okui, Abner Soffer, and Kohsaku Tobioka.  Heavy QCD axion at Belle II: Displaced and prompt signals.  Phys. Rev. D, 105(7):L071701, 2022, 2108.10331.
  * [43] Martin Bauer, Matthias Neubert, Sophie Renner, Marvin Schnubel, and Andrea Thamm.  Consistent Treatment of Axions in the Weak Chiral Lagrangian.  Phys. Rev. Lett., 127(8):081803, 2021, 2102.13112.
  * [44] J. K. Ahn et al.  Search for the $K_{L}\\!\to\\!\pi^{0}\nu\overline{\nu}$ and $K_{L}\\!\to\\!\pi^{0}X^{0}$ decays at the J-PARC KOTO experiment.  Phys. Rev. Lett., 122(2):021802, 2019, 1810.09655.
  * [45] Kazuya Aoki et al.  Extension of the J-PARC Hadron Experimental Facility: Third White Paper.  10 2021, 2110.04462.
  * [46] H. Nanjo.  Dark Sector search at KOTO experiment, 2021.  Contribution to “3rd meeting on Searches for Hidden Sectors at Kaon and Hyperon Factories”.
  * [47] Jan Jerhot, Babette Döbrich, Fatih Ertas, Felix Kahlhoefer, and Tommaso Spadaro.  ALPINIST: Axion-Like Particles In Numerous Interactions Simulated and Tabulated.  JHEP, 07:094, 2022, 2201.05170.
  * [48] Y.-C. Tung.  Status of $K_{L}\to\pi^{0}\nu\bar{\nu}$ Study at J-PARC KOTO Experiment, 2022.  Contribution to “41st International Conference on High Energy Physics, Bologna, Italy”.
  * [49] KOTO, LHCb, NA62/KLEVER, and US Kaon Interest Group.  Searches for new physics with high-intensity kaon beams.  In 2022 Snowmass Summer Study, 4 2022, 2204.13394.
  * [50] Torbjörn Sjöstrand, Stefan Ask, Jesper R. Christiansen, Richard Corke, Nishita Desai, Philip Ilten, Stephen Mrenna, Stefan Prestel, Christine O. Rasmussen, and Peter Z. Skands.  An introduction to PYTHIA 8.2.  Comput. Phys. Commun., 191:159–177, 2015, 1410.3012.
  * [51] Fatih Ertas and Felix Kahlhoefer.  On the interplay between astrophysical and laboratory probes of MeV-scale axion-like particles.  JHEP, 07:050, 2020, 2004.01193.
  * [52] Florian Domingo.  Decays of a NMSSM CP-odd Higgs in the low-mass region.  JHEP, 03:052, 2017, 1612.06538.
  * [53] Howard Georgi, David B. Kaplan, and Lisa Randall.  Manifesting the Invisible Axion at Low-energies.  Phys. Lett. B, 169:73–78, 1986.
  * [54] Eder Izaguirre, Tongyan Lin, and Brian Shuve.  Searching for Axionlike Particles in Flavor-Changing Neutral Current Processes.  Phys. Rev. Lett., 118(11):111802, 2017, 1611.09355.
  * [55] Nitsan Bar, Kfir Blum, and Guido D’Amico.  Is there a supernova bound on axions?  Phys. Rev. D, 101(12):123025, 2020, 1907.05020.
  * [56] E. M. Riordan et al.  A Search for Short Lived Axions in an Electron Beam Dump Experiment.  Phys. Rev. Lett., 59:755, 1987.
  * [57] Matthew J. Dolan, Torben Ferber, Christopher Hearty, Felix Kahlhoefer, and Kai Schmidt-Hoberg.  Revised constraints and Belle II sensitivity for visible and invisible axion-like particles.  JHEP, 12:094, 2017, 1709.00009.  [Erratum: JHEP 03, 190 (2021)].
  * [58] E. Cortina Gil et al.  HIKE, High Intensity Kaon Experiments at the CERN SPS.  11 2022, 2211.16586.
  * [59] 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:092004, 2009, 0903.0030.
  * [60] Eduardo Cortina Gil et al.  Search for $\pi^{0}$ decays to invisible particles.  JHEP, 02:201, 2021, 2010.07644.
  * [61] Stefania Gori, Gilad Perez, and Kohsaku Tobioka.  KOTO vs. NA62 Dark Scalar Searches.  JHEP, 08:110, 2020, 2005.05170.
  * [62] J. P. Lees et al.  Search for an Axionlike Particle in $B$ Meson Decays.  Phys. Rev. Lett., 128(13):131802, 2022, 2111.01800.
  * [63] Jonathan L. Feng, Iftah Galon, Felix Kling, and Sebastian Trojanowski.  Axionlike particles at FASER: The LHC as a photon beam dump.  Phys. Rev. D, 98(5):055021, 2018, 1806.02348.
  * [64] Nikita Blinov, Elizabeth Kowalczyk, and Margaret Wynne.  Axion-like particle searches at DarkQuest.  JHEP, 02:036, 2022, 2112.09814.
  * [65] Jan Jerhot.  jjerhot/alpinist: v1.2.0, March 2023.
  * [66] T. Masuda et al.  Long-lived neutral-kaon flux measurement for the KOTO experiment.  PTEP, 2016(1):013C03, 2016, 1509.03386.