# Captured molecules could make a Bose star visible

V. V. Flambaum<EMAIL_ADDRESS>I. B. Samsonov
<EMAIL_ADDRESS>School of Physics, University of New South Wales,
Sydney 2052, Australia

###### Abstract

A Bose star passing through cold molecular clouds may capture atoms, molecules
and dust particles. The observational signature of such an event would be a
relatively small amount of matter that is gravitationally bound. This binding
may actually be provided by invisible dark matter forming the Bose star. We
may expect a relative excess of heavier atoms, molecules, and solid dust
compared to the content of giant cold molecular clouds since the velocity of
heavy particles at a given temperature is lower and it may be small compared
to the escape velocity,
$v_{\mathrm{rms}}=\sqrt{3k_{\mathrm{B}}T/m_{\mathrm{gas}}}\ll
v_{\mathrm{esc}}=\sqrt{2GM/R}$. Finally, the velocity of this captured matter
cloud may correlate with the expected velocity of free dark matter particles
(e.g. expected axion wind velocity relative to Earth).

##  I Introduction

Bose stars are assumed to be formed from cold dark matter fields (scalars,
axions) which may form Bose condensate and solitons bound mainly by the
gravitational attraction [1, 2]. Bose condensation temperature $T_{B}\propto
n^{2/3}/m$, so in thermal equilibrium axions (or scalars) form Bose condensate
for sufficiently small mass $m$ and large density $n$.

A simple estimate, similar to that for hydrogen atom, makes relation between
mass $M\approx Nm$ and size $R$ of Bose star, where $N$ is the number of
bosons and $m$ is their mass. Estimating kinetic energy of the ground state
using uncertainty relation for minimal momentum $p\sim\hbar/r$ and estimating
gravitational energy we obtain equation for the non-relativistic energy of
Bose star:

$E=a_{1}N\frac{\hbar^{2}}{mr^{2}}-a_{2}\frac{N(N-1)}{2}\frac{Gm^{2}}{r},$ (1)

where $a_{1}$, $a_{2}$ are constants $\sim 1$. Minimizing energy by varying
size $r$ or using virial theorem ($U=-2E_{k}$, $E=-E_{k}$) we obtain size of
gravitationally bound Bose star

$R=\frac{4a_{1}\hbar^{2}}{a_{2}(N-1)Gm^{3}},$ (2)

binding energy

$E=N(N-1)^{2}\frac{a_{2}^{2}G^{2}m^{5}}{16a_{1}\hbar^{2}},$ (3)

and number density

$n\sim\frac{N}{V}\sim 0.004N(N-1)^{3}\frac{G^{3}m^{9}}{\hbar^{6}},$ (4)

where the star volume is $V\sim 4\pi R^{3}/3$. We see very strong dependence
on the boson mass $m$. We also see that the density $n$ in the ground state is
proportional to $N^{4}$ so any, even small non-gravitational self interaction
between bosons may become important for large $N$. We may add the self-
interaction contribution of the self-interaction term $\lambda\phi^{4}$ to the
Bose star energy

$E=a_{1}N\frac{\hbar^{2}}{mr^{2}}-a_{2}\frac{N(N-1)}{2}\frac{Gm^{2}}{r}+a_{3}\lambda\frac{N(N-1)}{2m^{2}r^{3}},$
(5)

where $a_{1}$, $a_{2}$ $a_{3}$ are constants $\sim 1$. Using $dE/dr=0$, this
expression gives quadratic equation for $r$ with two solutions. For the QCD
axion, the selfinteraction constant $\lambda=-m^{2}/(24f_{a}^{2})$ is
expressed via the axion decay constant $f_{a}$ (see e.g. [3, 4]). In the case
of $\lambda<0$ the energy tends to minus infinity for $r\to 0$, so the Bose
star is metastable. Solution with larger $r$ corresponds to a minimum whereas
the smaller $r$ solution is maximum. More accurate analytical calculations
using variational approach will be given in the next section.

Distribution of ordinary matter is given by Boltzmann distribution in the
gravitational field of Bose star. Velocity distribution is given by Maxwell
distribution. We can add molecular matter terms $E_{m}$ to the energy $E$:

$\displaystyle E_{m}=$ $\displaystyle
a_{4}N_{m}\frac{5}{2}kT(r)-a_{5}N_{M}N\frac{Gm_{m}m}{r}$ (6) $\displaystyle-
a_{6}\frac{N_{m}(N_{m}-1)}{2}\frac{Gm_{m}^{2}}{r}\,.$

If the total mass of the molecular cloud bound to the Bose star,
$M_{m}=N_{m}m_{m}$, is relatively small as compared with the mass of the Bose
star, $N_{m}m_{m}\ll Nm$, we may assume that in the leading approximation
molecular distribution does not affect Bose star parameters. Then we may
include gravitational effect of molecules using iterations.

Time of formation of Bose stars very strongly depends on the Boson mass $m$.
For $m>10^{-22}$ eV the formation time is smaller than the lifetime of the
Universe 13.7 billion years (see e.g. [3]).

Total mass of the Bose star $M$ is strongly model dependent and depends on
time. There are estimates of a typical scale for expected mass $M$. The QCD
axion dark matter field with axion mass $m\sim 10^{-5}$ eV makes Bose star
with mass $M\sim 10^{-15}M_{\odot}$ [3, 5, 6]. Fuzzy dark matter with particle
mass $m\sim 10^{-22}$ eV may compose a Bose star with mass $M\sim
10^{8}M_{\odot}$ [7, 8]. Bose star accumulating certain critical mass becomes
unstable, collapses and decays to relativistic jets [9].

Axions and scalars are neutral, so in zero approximation Bose star is
transparent and is hardly observable by telescopes. In principle, Bose stars
may be observed using gravitational lensing of passing light and axion decays
to photons. We discuss a different possibility: Bose star may capture atoms,
molecules and dust. Radiation and absorption of light by this matter could
make Bose stars visible.

##  II Variational approach

The variational approach allows one to get an approximate analytic solution to
the particle density which minimizes total energy of the system of interacting
particles. We consider $N$ bosonic DM particles with position vectors
$\vec{r}_{i}$ and momenta $\vec{p}_{i}$. Given that each particle has mass
$m$, they can interact gravitationally. The corresponding non-relativistic
Hamiltonian is

$H=-\frac{\hbar^{2}}{2m}\sum_{i=1}^{N}\nabla_{k}^{2}-\sum_{i<j}^{N}\frac{Gm^{2}}{r_{ij}}\,,$
(7)

where $G$ is the gravitational constant and
$r_{ij}\equiv|\vec{r}_{i}-\vec{r}_{j}|$.

In the ground state, all particles are in the same state with the wave
function $\psi(\vec{r}_{i})$, and the full system is described by a product of
these wave functions,

$\Psi=\prod_{i=1}^{N}\psi(\vec{r}_{i})\,.$ (8)

Below we consider two approximate analytic solutions to the function
$\psi(\vec{r})$.

###  II.1 Exponential ansatz

The ground state wave function may be searched for within the ansatz

$\psi(r)=\frac{1}{\sqrt{\pi}}R^{-3/2}e^{-r/R}\,,$ (9)

where $R$ is a free parameter. The function $\psi(r)$ is the 1s solution of
the stationary Schrödinger equation with the potential
$V=-\frac{\hbar^{2}}{mRr}$ and energy $E_{0}=-\frac{\hbar^{2}}{2mR^{2}}$.

Making use of the function (8) we find the expectation value of the
Hamiltonian (7):

$E=\langle\Psi|H|\Psi\rangle=\frac{N}{2}\left[\frac{\hbar^{2}}{mR^{2}}-\frac{5}{8}\frac{Gm^{2}}{R}(N-1)\right].$
(10)

This energy reaches its minimum value

$E_{\text{min}}=-\frac{25}{512}\frac{G^{2}m^{5}}{\hbar^{2}}N(N-1)^{2}$ (11)

at

$R=\frac{16}{5}\frac{\hbar^{2}}{Gm^{3}(N-1)}\,.$ (12)

The energy of the solution (11) agrees at 90% level with the corresponding
numerical result in Ref. [10]. The root-mean-squared (RMS) radius of this
solution is

$r_{\text{rms}}=\left(\langle\Psi|r_{i}^{2}|\Psi\rangle\right)^{1/2}=\sqrt{3}R\,.$
(13)

As an example, consider a Bose star with mass of Bennu asteroid, $M=7\times
10^{10}$ kg composed of QCD axions with mass $m=10^{-5}$ eV. It consists of
$N=4\times 10^{51}$ particles. Its size is

$r_{\text{rms}}=4\times 10^{10}\,\mbox{km}.$ (14)

Note that the diameter of the Bennu asteroid is 0.5 km.

The RMS speed:

$v_{\text{rms}}=\frac{1}{m}\left(\langle\Psi|\vec{p}_{i}^{2}|\Psi\rangle\right)^{1/2}=\frac{\hbar}{mR}\,.$
(15)

For the axion star considered above, it is $v\approx 8\times 10^{-16}c$. Thus,
the non-relativistic approximation is justified. However, for solutions with
boson mass $m\sim 10^{-5}$ eV and total mass of order the Solar mass, a
relativistic treatment is needed.

###  II.2 Gaussian ansatz

The simple exponential ansatz (9) for the wave function provides a rather low
accuracy for the energy. A better accuracy may be achieved within the Gaussian
ansatz [11]

$\psi(r)=\left(\frac{2}{\pi R^{2}}\right)^{3/4}e^{-r^{2}/R^{2}}\,,$ (16)

where $R$ is the parameter to be found. Making use of this wave function, we
find the expectation value of the Hamiltonian (7):

$E=\langle
H\rangle=\frac{3N\hbar^{2}}{2mR^{2}}-\frac{N(N-1)Gm^{2}}{\sqrt{\pi}R}\,.$ (17)

This energy reaches its minimum value

$E_{\text{min}}=-\frac{N(N-1)^{2}}{6\pi}\frac{G^{2}m^{5}}{\hbar^{2}}$ (18)

at

$R=\frac{3\sqrt{\pi}\hbar^{2}}{(N-1)Gm^{3}}\,.$ (19)

The energy (18) agrees at 98% level with the numerical solution presented in
Ref. [10]. Therefore, we will consider only the Gaussian ansatz in what
follows.

The RMS radius calculated with the wave function (16) is

$r_{\text{rms}}=\frac{\sqrt{3}R}{2}=\frac{3\sqrt{3\pi}\hbar^{2}}{2(N-1)Gm^{3}}\,,$
(20)

and the RMS speed is

$v_{\text{rms}}=\frac{\sqrt{3}\hbar}{mR}\,.$ (21)

The non-relativistic approximation is applicable for $v_{\mathrm{rms}}\lesssim
0.1c$. Making use of equations (19) and (21), this relation translates into
the constraint on the particle number and mass,

$Nm^{2}<\frac{0.1\hbar c\sqrt{3\pi}}{G}\approx 1.5\times
10^{-10}\,\mathrm{g}^{2}\,.$ (22)

###  II.3 Gravitational potential

The Bose star solution is given by the product of $N$ identical wave functions
(16) each describing a particle with mass $m$. They create a gravitational
potential $V(r)$ obeying the Poisson equation

$\nabla^{2}V(r)=-4\pi G\rho(r)\,,$ (23)

where the mass density is

$\rho=Nm\psi^{\dagger}\psi=Nm\left(\frac{2}{\pi
R^{2}}\right)^{3/2}e^{-2r^{2}/R^{2}}\,.$ (24)

The solution to Eq. (23) reads

$\displaystyle V(r)$ $\displaystyle=-GNm\frac{\mathrm{erf}(\sqrt{2}r/R)}{r}$
$\displaystyle\approx-
GNm\left\\{\begin{array}[]{ll}\frac{2\sqrt{2}}{\sqrt{\pi}R}-\frac{4\sqrt{2}r^{2}}{3\sqrt{\pi}R^{3}}&(r\ll
R)\,,\\\ \frac{1}{r}&(r\gg R)\,.\end{array}\right.$ (27)

On the boundary $r=r_{\mathrm{rms}}$, the gravitational potential is

$V(r_{\mathrm{rms}})=-\frac{2GNm}{\sqrt{3}R}\mathrm{erf}\left(\sqrt{3/2}\right)\approx-0.2N^{2}\frac{G^{2}m^{4}}{\hbar^{2}}\,.$
(28)

Thus, the escape velocity on the boundary is

$v_{\mathrm{esc}}=\sqrt{2|V(r_{\mathrm{rms}})|}\approx 0.6NGm^{2}/\hbar\,.$
(29)

Consider, for example, the QCD axion with mass $m=10^{-5}$ eV and the total
mass comparable with the mass of the Earth, $Nm=M_{\oplus}$. The size of this
object is about $r_{\mathrm{rms}}\approx 0.4$ m. The corresponding escape
velocity at the boundary is $v_{\mathrm{esc}}\approx 0.14c$.

If the DM particle mass is $m=10^{-9}$ eV, the corresponding RMS radius and
escape velocity are $r_{\mathrm{rms}}\approx 6.3R_{\oplus}$,
$v_{\mathrm{esc}}\approx 4.3$ km/s.

###  II.4 Adding $\lambda\phi^{4}$ selfinteraction

The scalar field $\phi$ is related to the wave function (8) as

$\phi=\frac{1}{\sqrt{m}}\Psi\,.$ (30)

Making use of the Gaussian ansatz (16), we find

$\lambda\int\phi^{4}d^{3}r=\frac{\lambda}{\pi^{3/2}m^{2}R^{3}}\,.$ (31)

In the contribution to the energy, this expression should be multiplied by
$N(N-1)/2$ which accounts for the pairwise particle interaction. Thus, the
total energy of this field configuration reads

$\displaystyle E$
$\displaystyle=E_{\text{kin}}+E_{\text{grav}}+E_{\lambda}\,,$ (32a)
$\displaystyle E_{\text{kin}}$ $\displaystyle=\frac{3N\hbar^{2}}{2mR^{2}},$
(32b) $\displaystyle E_{\text{grav}}$
$\displaystyle=-\frac{N(N-1)Gm^{2}}{\sqrt{\pi}R}\,,$ (32c) $\displaystyle
E_{\lambda}$
$\displaystyle=\frac{N(N-1)}{2}\frac{\lambda\hbar^{3}}{\pi^{3/2}cm^{2}R^{3}}\,.$
(32d)

This expression reaches its extrema at

$R_{\pm}=\frac{3\sqrt{\pi}\hbar^{2}}{2Gm^{3}(N-1)}\left(1\pm\sqrt{1+\frac{2\lambda
Gm^{2}(N-1)^{2}}{3\pi^{2}\hbar c}}\right).$ (33)

The case $\lambda<0$ corresponds to unstable solution because the energy (32)
is unbounded from below in this case. However, there is a potential barrier
separating energy minimum at $R_{+}$ in Eq. (33) from the collapse area close
to $r=0$, with the barrier maximum at $R_{-}$. In the case of small $\lambda$,
we obtain minimum of energy

$E_{\text{min}}=-\frac{N(N-1)^{2}}{6\pi}\frac{G^{2}m^{5}}{\hbar^{2}}+O(\lambda)$
(34)

at

$R=\frac{3\sqrt{\pi}\hbar^{2}}{Gm^{3}(N-1)}+\frac{\lambda\hbar(N-1)}{2\pi^{3/2}mc}+O(\lambda^{2})\,,$
(35)

and maximum at small $R\propto\lambda$,

$R=-\frac{\lambda\hbar(N-1)}{2\pi^{3/2}mc}+O(\lambda^{2})\,.$ (36)

In the case of small $\lambda$ this barrier will be sufficiently high to make
the decay time of the metastable Bose star exponentially long, making this
kind of Bose star potentially observable.

For positive selfinteraction constant, $\lambda>0$, only “+” solution in Eq.
(33) is physical. If $\lambda$ is sufficiently small, the leading contribution
has the same form as in Eq. (36). For strong selfinteraction, $\frac{2\lambda
Gm^{2}(N-1)^{2}}{3\pi^{2}\hbar c}\gg 1$, the parameter $R$ is independent of
$N$,

$R=\frac{\hbar^{3/2}}{m^{2}}\sqrt{\frac{3\lambda}{2\pi cG}}\,.$ (37)

##  III Ideal gas distribution around the Bose star

Let us consider an ideal gas with particles of mass $m_{\mathrm{gas}}$ in the
gravitational potential $V(r)<0$ at temperature $T$. Assuming that the total
mass of the gas is sufficiently small so that the gravitational potential is
not changed significantly, we can write the density distribution in the gas
cloud,

$\rho(r)=\rho_{0}\exp\left(\frac{m_{\mathrm{gas}}V(r)}{k_{\mathrm{B}}T}\right)\,.$
(38)

The total mass of the gas cloud is

$M_{\mathrm{gas}}=4\pi\rho_{0}\int_{0}^{\infty}\exp\left(\frac{m_{\mathrm{gas}}V(r)}{k_{\mathrm{B}}T}\right)r^{2}dr\,.$
(39)

The gas at every point should follow the Maxwell-Boltzmann distribution

$f(v)=4\pi v^{2}\left(\frac{m_{\mathrm{gas}}}{2\pi
k_{\mathrm{B}}T}\right)^{3/2}e^{-\frac{m_{\mathrm{gas}}v^{2}}{2k_{\mathrm{B}}T}}\,.$
(40)

This distribution, however, should be truncated at $v=v_{\mathrm{esc}}$. Note
that $v_{\mathrm{rms}}=\sqrt{3k_{\mathrm{B}}T/m_{\mathrm{gas}}}$.

At the Earth’s surface, the temperature is about $T=300$ K. Empirically, we
find the condition for the gas molecules to be captured by the gravitational
field of the Earth:

$v_{\mathrm{rms}}\lesssim 0.1v_{\mathrm{esc}}\,.$ (41)

This corresponds to the escape probability suppressed by the exponential
factor in the distribution (40)
$\exp{(-{\frac{3}{2}\frac{v_{\mathrm{esc}}^{2}}{v_{\mathrm{rms}}^{2}}})}\lesssim\exp{(-150)}$.
For the Earth, this condition implies that the gas molecules with mass number
$A<6$ escape from the gravitational attraction. This is confirmed by the fact
that the fractions of hydrogen and helium gases in the Earth’s atmosphere are
very low.

In the case of Bose stars, making use of Eq. (29) the above condition may be
cast in the form

$m_{\mathrm{gas}}>\frac{830k_{\mathrm{B}}T\hbar^{2}}{N^{2}G^{2}m^{4}}\,,$ (42)

Representing the gas particle mass as $m_{\mathrm{gas}}=Am_{p}$, with $m_{p}$
the proton mass, we find the condition on the atomic number $A$ of gas
particles captured near the surface of the Bose star,

$\displaystyle A>$
$\displaystyle\frac{830k_{\mathrm{B}}T\hbar^{2}}{N^{2}G^{2}m^{4}m_{p}}$
$\displaystyle=1.35\cdot
10^{-29}\left(\frac{M_{\odot}}{M}\right)^{2}\left(\frac{\mathrm{eV}}{m}\right)^{2}\frac{T}{10K}\,.$
(43)

Figure 1: Relation between the mass and radius of the Bose star $M$ (in units
of solar mass and radius) and the scalar particle mass $m$, assuming that this
object can gravitationally bound gas particles with atomic mass numbers $A$
exceeding 1, 10 and 100. The gas is assumed to have temperature $T=10$ K.

The temperature in cold molecular clouds in our Galaxy is typically in the
range from 10 to 20 K. Assuming $T=10$ K, we estimate the possibility of
capture of gases in the fuzzy dark matter with $m=10^{-22}$ eV and total mass
of the object $Nm=10^{7}M_{\odot}$:

$A>13\,.$ (44)

Thus, this object can accumulate only heavy atoms and molecules, as well as
dust particles.

In Fig. 1, we present plots of mass (and radius) of the Bose star as a
function of scalar particle mass $m$ which can capture molecules with atomic
mass numbers $A=1$, 10, and 100.

If radius $R$ of Bose star is small, there will be too little matter inside
the star and impossible to observe it. However, matter may be captured to
distant orbits, but this will require bigger mass $M$ of the star, since the
escape velocity will be smaller; for matter at distance $r$, the required star
mass is $M(r)=M_{0}R/r$, where $M_{0}$ is the star mass required to capture
matter inside.

Cold molecular gas and dust may also be captured by large axion clusters.

If matter capture conditions are fulfilled inside Bose star for all A, there
will be violation of such conditions outside Bose star at some distance $r$
from the centre. Firstly, this condition will be violated for A=1, going
further, for higher A. Therefore, there will be separation of heavy and light
elements outside the Bose star. This is important for Bose star of large mass
and small radius.

##  IV Summary

We applied the analytic Gaussian solution for Bose star density to study the
distribution of interstellar gas particles within the gravitational field of a
Bose star, assuming the gas does not significantly perturb the star’s
gravitational field. We examined the conditions under which interstellar gas
particles with a given atomic number $A$ could escape the gravitational
attraction of a Bose star with a total mass $M$ composed of scalar dark matter
particles of mass $m$. Our findings indicate that, for a range of values of
$m$ and $M$, the Bose star can capture atoms, molecules, and dust particles.

A small cloud of interstellar gas cannot form a gravitationally bound state
without an additional gravitational source. Thus, the observation of such a
gravitationally bound gas cloud would suggest it is held by the gravitational
field of a Bose star or a dark matter cluster. Due to differences in escape
velocities between light and heavy molecules, molecular clouds captured by a
Bose star may exhibit higher metallicity compared to typical giant molecular
clouds. The discovery of small, gravitationally bound molecular clouds with
high metallicity in the interstellar medium would provide compelling evidence
supporting this dark matter model.

Acknowledgements.— We are grateful to John Webb for interest in this work and
useful discussions. The work was supported by the Australian Research Council
Grants No. DP230101058 and No. DP200100150.

## References

  * Ruffini and Bonazzola [1969] R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).
  * Tkachev [1986] I. Tkachev, Sov. Astron. Lett. 12, 305–308 (1986).
  * Levkov  _et al._ [2017] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Phys. Rev. Lett. 118, 011301 (2017).
  * Marsh [2016] D. J. E. Marsh, Physics Reports 643, 1 (2016).
  * Eggemeier and Niemeyer [2019] B. Eggemeier and J. C. Niemeyer, Phys. Rev. D 100, 063528 (2019).
  * Chen  _et al._ [2021] J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh, and J. C. Niemeyer, Phys. Rev. D 104, 083022 (2021), arXiv:2011.01333 [astro-ph.CO] .
  * Schive  _et al._ [2014] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Veltmaat and Niemeyer [2016] J. Veltmaat and J. C. Niemeyer, Phys. Rev. D 94, 123523 (2016), arXiv:1608.00802 [astro-ph.CO] .
  * Visinelli  _et al._ [2018] L. Visinelli, S. Baum, J. Redondo, K. Freese, and F. Wilczek, Physics Letters B 777, 64 (2018).
  * Membrado  _et al._ [1989] M. Membrado, A. F. Pacheco, and J. Sañudo, Phys. Rev. A 39, 4207 (1989).
  * Chavanis [2011] P.-H. Chavanis, Phys. Rev. D 84, 043531 (2011).


††thanks: 0000-0002-3143-1976

# Searching for axions with kaon decay at rest

Yohei Ema<EMAIL_ADDRESS>William I. Fine Theoretical Physics Institute,
University of Minnesota, Minneapolis, MN 55455, USA  School of Physics and
Astronomy, University of Minnesota, Minneapolis, MN 55455, USA  Zhen Liu
<EMAIL_ADDRESS>School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA  Ryan Plestid<EMAIL_ADDRESS>Walter Burke
Institute for Theoretical Physics, California Institute of Technology,
Pasadena, CA 91125, USA

###### Abstract

We describe a novel search strategy for axions (or hadronically coupled axion-
like particles) in the mass range of $m_{a}\lesssim 350\,{\rm MeV}$. The
search relies on kaon decay at rest, which produces a mono-energetic signal in
a large volume detector (e.g. a tank of liquid scintillator) from axion decays
$a\rightarrow\gamma\gamma$ or $a\rightarrow e^{+}e^{-}$. The decay modes
$K^{+}\to\pi^{+}a$ and $a\to\gamma\gamma$ are induced by the axion’s coupling
to gluons, which is generic to any model which addresses the strong CP
problem. We recast a recent search from MicroBooNE for $e^{+}e^{-}$ pairs, and
study prospects at JSNS2 and other near-term facilities. We find that JSNS2
will have world-leading sensitivity to hadronically coupled axions in the mass
range of $40\,\mathrm{MeV}\lesssim m_{a}\lesssim 350\,\mathrm{MeV}$.

††preprint: UMN-TH-4221/23††preprint: FTPI-MINN-23-13††preprint: CALT-
TH/2023-028

Introduction: The neutron’s electron dipole moment (EDM), $d_{n}\lesssim
2\times 10^{-26}~{}e~{}{\rm cm}$ Graner  _et al._ (2016); Abel  _et al._
(2020), is ten orders of magnitude smaller than its naive estimate of $5\times
10^{-16}e~{}{\rm cm}$. This is unexpected since charge-parity (CP) is violated
within the Standard Model (SM), and therefore no fundamental symmetry forbids
a neutron EDM. This is the strong CP problem, and its microscopic origin can
be traced to the minuscule value of the QCD $\bar{\theta}$ parameter,
$\bar{\theta}\lesssim 10^{-10}$, which controls the unique CP-violating QCD
coupling in the SM. Axions are a popular solution that provide a dynamical
mechanism for the relaxation of $\bar{\theta}$ Peccei and Quinn (1977a, b);
Weinberg (1978); Wilczek (1978); if the axion’s potential is generated
exclusively by QCD, it naturally aligns the axion’s ground state with
$\bar{\theta}=0$.

Axions generically suffer from the so-called quality problem Kamionkowski and
March-Russell (1992); Barr and Seckel (1992); Ghigna  _et al._ (1992); Holman
_et al._ (1992) wherein high energy contributions to the axion potential can
(and often do) displace the minimum away from $\bar{\theta}=0$. These problem
does not arise if the QCD potential is “strenghtened”, for instance by
introducing a mirror QCD sector. These non-minimal axion models produce a
potential that is robust against the high energy corrections discussed above
and often predict heavier axions as compared to minimal axion models
Dimopoulos (1979); Tye (1981); Holdom and Peskin (1982); Flynn and Randall
(1987); Rubakov (1997); Berezhiani  _et al._ (2001); Hook (2015); Fukuda  _et
al._ (2015); Gherghetta  _et al._ (2016); Dimopoulos  _et al._ (2016); Agrawal
and Howe (2018a, b); Gaillard  _et al._ (2018); Lillard and Tait (2018);
Fuentes-Martín  _et al._ (2019); Csáki  _et al._ (2020); Hook  _et al._
(2020); Gherghetta  _et al._ (2020); Gherghetta and Nguyen (2020); Valenti
_et al._ (2022); Kivel  _et al._ (2022).

Recent investigations have reignited interest in these so-called heavy axion
models with $m_{a}\gtrsim 1\,{\rm MeV}$. More generally, there is a broad
interest in axion like particles (ALPs) Bauer  _et al._ (2017); Irastorza and
Redondo (2018) which may serve as IR messengers of UV completions such as a
string landscape Arvanitaki  _et al._ (2010); Jaeckel and Ringwald (2010). A
wide range of search strategies have been proposed ranging from beam dumps, to
flavor facilities, and collider experiments Essig  _et al._ (2010); Döbrich
_et al._ (2016); Dolan  _et al._ (2017); Döbrich  _et al._ (2019); Harland-
Lang  _et al._ (2019); Dent  _et al._ (2020); Brdar  _et al._ (2021); Jaeckel
_et al._ (2013); Mimasu and Sanz (2015); Jaeckel and Spannowsky (2016);
Izaguirre  _et al._ (2017); Knapen  _et al._ (2017); Bauer  _et al._ (2017);
Brivio  _et al._ (2017); Mariotti  _et al._ (2018); Cid Vidal  _et al._
(2019); Beacham  _et al._ (2020); Alonso-Álvarez  _et al._ (2019); Ebadi  _et
al._ (2019); Gavela  _et al._ (2020, 2020); Altmannshofer  _et al._ (2020);
Gershtein  _et al._ (2021); Di Luzio  _et al._ (2020); Knapen  _et al._
(2022); Co  _et al._ (2023); Acciarri  _et al._ (2023).

In this work, we point out that kaon decay at rest (KDAR) offers a powerful
probe of axions (or hadronically coupled ALPs) in the mass range of
$m_{a}=40\,\mathchar 45\relax\mathchar 45\relax\,350~{}{\rm MeV}$. Unlike
previous studies of KDAR Coloma  _et al._ (2022); Batell  _et al._ (2019) the
axions we consider are naturally coupled to quarks and/or gluons such that the
hadronic decays of kaons serve as a powerful axion factory. The signal is a
visible decay of either $a\rightarrow\gamma\gamma$ or $a\rightarrow
e^{+}e^{-}$. The dominant SM kaon decay modes are
$K^{+}\rightarrow\mu^{+}\nu_{\mu}$, as well as
$K^{+}\rightarrow\pi^{0}e^{+}\nu_{e}$,111Negatively charged $K^{-}$ decays
are subdominant to $K^{-}$ capture on nuclei. This cannot occur for $K^{+}$
because it does not form bound states with nuclei. and both channels produce a
neutrino flux that can be measured Spitz (2014); Nikolakopoulos  _et al._
(2021). In addition to predictions of hadronic cascade simulations, this
provides an experiment with an in situ measurement of their KDAR population.
The axion rate is then calculable, and a counting experiment can be performed.

For concreteness, we consider an axion coupled to Standard Model (SM) gauge
bosons via

$\begin{split}\mathcal{L}_{a}&=\frac{1}{2}(\partial
a)^{2}-\frac{m_{a}^{2}}{2}a^{2}+c_{GG}\frac{\alpha_{s}}{4\pi}\frac{a}{f}G_{\mu\nu}\tilde{G}^{\mu\nu}\\\
&+c_{WW}\frac{\alpha_{2}}{4\pi}\frac{a}{f}W_{\mu\nu}\tilde{W}^{\mu\nu}+c_{BB}\frac{\alpha_{Y}}{4\pi}\frac{a}{f}B_{\mu\nu}\tilde{B}^{\mu\nu},\end{split}$
(1)

at a high scale, where $a$ is the axion, and $\alpha_{s}=g_{s}^{2}/4\pi$,
$\alpha_{2}=g^{2}/4\pi$ and $\alpha_{Y}={g^{\prime}}^{2}/4\pi$. The coupling
$c_{GG}$ can be absorbed into $f$ and we define the axion decay constant as
$f_{a}=f/2c_{GG}$. As two useful benchmarks, we focus on the so-called “gluon
dominance” ($c_{GG}\neq 0$, $c_{WW}=c_{BB}=0$), and “co-dominance”
($c_{GG}=c_{WW}=c_{BB}\neq 0$) scenarios following Kelly  _et al._ (2021);
Jerhot  _et al._ (2022); Afik  _et al._ (2023) (using the relative
normalization of $c_{BB}$ in Jerhot  _et al._ (2022); Afik  _et al._ (2023)).
We do not consider, e.g. axion couplings to quarks, though their inclusion is
straightforward. In what follows, we discuss the theory of
$K^{+}\rightarrow\pi^{+}a$, and project sensitivities for JSNS2 and a
MicroBooNE $\gamma\gamma$ search.

Kaon decay at rest: In models with a hadronically coupled axion,
$K^{+}\rightarrow\pi^{+}a$ serves as a powerful axion production channel.
Since axions are long-lived, they decay in flight to visible final states. The
resulting signal is a mono-energetic peak at an energy of

$E_{a}=\frac{m_{K}^{2}+m_{a}^{2}-m_{\pi}^{2}}{2m_{K}}~{}.$ (2)

The branching ratio for $K^{+}\rightarrow\pi^{+}a$ can be reliably predicted
in chiral perturbation theory Bauer  _et al._ (2021a), with the result
[neglecting terms of $\mathcal{O}(m_{a}^{2}/m_{K}^{2})$ and
$\mathcal{O}(m_{\pi}^{2}/m_{K}^{2})$]

$\displaystyle\mathrm{BR}(K^{+}\to\pi^{+}a)$
$\displaystyle=\frac{\tau_{K^{-}}}{\tau_{K_{S}}}\frac{f_{\pi}^{2}}{8f_{a}^{2}}\times\mathrm{BR}(K_{S}\to\pi^{+}\pi^{-})~{}.$
(3)

Equation 3 receives corrections from finite axion and pion mass effects in
both the matrix element and phase space, both of which are taken into account
in our numerical estimates (see Bauer  _et al._ (2021a) for the complete
expression). In the co-dominance case, the axion coupling to the $W$-boson
induces an extra contribution Izaguirre  _et al._ (2017). However, this
contribution is negligible for $c_{WW}=c_{GG}$ and we do not consider it
further in what follows.

Figure 1:  (Gluon Dominance) Sensitivities of MicroBooNE and JSNS2 compared
with existing limits and other projected sensitivities when all couplings are
induced by a gluon coupling $c_{GG}$ at a high scale. The MicroBooNE
sensitivity is cut at 210 MeV because that is the range that appears in
Abratenko  _et al._ (2021). Existing limits include constraints from SN1987A
Chang  _et al._ (2018); Ertas and Kahlhoefer (2020) and cosmology Depta  _et
al._ (2020) adapted from Kelly  _et al._ (2021), Kaon decays (E949 Artamonov
_et al._ (2009) and NA62 Cortina Gil  _et al._ (2021)), and beam dump searches
(CHARM Bergsma  _et al._ (1985) and NuCal Blumlein  _et al._ (1991)), with
data adapted from Jerhot  _et al._ (2022). We also show projected sensitivity
of DUNE Kelly  _et al._ (2021). Other projected sensitivities not shown here
include DarkQuest Aidala  _et al._ (2019); Berlin  _et al._ (2018); Blinov
_et al._ (2022), FASER Ariga  _et al._ (2019), KOTO Yamanaka (2012); Afik  _et
al._ (2023), and SHiP Anelli  _et al._ (2015). We have re-derived constraints
from E949 and NA62 and our results agree with Afik  _et al._ (2023) (but
disagree with Ertas and Kahlhoefer (2020) and therefore the curves in Kelly
_et al._ (2021)).

Axions produced from KDAR can decay inside nearby detectors, either to
$\gamma\gamma$ or $e^{+}e^{-}$. The number of axions produced at a KDAR source
is given by

$N_{a}=\frac{{\rm BR}(K^{+}\rightarrow\pi^{+}a)~{}}{{\rm
BR}(K^{+}\rightarrow\mu^{+}\nu_{\mu})}N_{\nu_{\mu}}~{},$ (4)

where $N_{\nu_{\mu}}$ is the number of muon neutrinos from KDAR that are
produced in the beam stop. In the limit of a long decay length,
$\lambda_{a}\gg L$ where $L$ is the distance from the KDAR source to the
detector, the number of axions that decay in the detector is given by222 In
our numerical computation, we do not rely on this approximation. Instead, we
include the finite axion decay length properly to obtain the upper limit of
the sensitivity correctly.

$N_{a\rightarrow{\rm vis}}=N_{a}\times\frac{1}{4\pi
L^{2}}\frac{V}{\lambda_{a}(m_{a})}~{},$ (5)

where $V$ is the volume of the detector,
$\lambda_{a}=\beta_{a}\gamma_{a}\tau_{a}$ is the decay length of the axion in
the lab frame, and we have assumed all final states of the axion decay are
visible.

Assuming $a\rightarrow\gamma\gamma$ dominates over $a\rightarrow ee$, which is
generically true for theories without direct axion-electron couplings, the
lifetime of the axion is given by

$\tau_{a}^{-1}=\frac{\alpha^{2}m_{a}^{3}}{256\pi^{3}f^{2}_{a}}\absolutevalue{c_{\gamma\gamma}^{\rm
eff}}^{2}~{},$ (6)

where $\alpha$ is the electromagnetic fine structure constant. The effective
coupling to the photon is given by Bauer  _et al._ (2017); Blinov  _et al._
(2022)

$c_{\gamma\gamma}^{\rm eff}\approx
c_{\gamma\gamma}-\quantity(\frac{5}{3}+\frac{m_{\pi}^{2}}{m_{\pi}^{2}-m_{a}^{2}}\frac{m_{d}-m_{u}}{m_{u}+m_{d}})~{},$
(7)

where $c_{\gamma\gamma}=0$ in the gluon dominance case and
$c_{\gamma\gamma}=2$ in the co-dominance case, and our constraints are
expressible in terms of $1/f_{a}^{4}$. Equation 7 assumes a two-flavor
approximation, which is reasonably accurate for the masses we consider,
$m_{a}\lesssim 350~{}{\rm MeV}$. For these masses we never approach regions of
resonant mixing with $\eta$ and $\eta^{\prime}$; more complete three flavor
expressions can be found in Appendix B of Blinov  _et al._ (2022). The
constraints can be easily re-scaled for other model-dependent choices (e.g.
with the axion couplings to the up and down quarks). It is also
straightforward to project sensitivities for models where $a\rightarrow
e^{+}e^{-}$ is the dominant decay channel.

For $m_{a}$ close to $m_{\pi}$ cancellations can occur such that
$c_{\gamma\gamma}^{\rm eff}$ vanishes. In this case the decay length of the
axion is set by $a\rightarrow e^{+}e^{-}$. Such a coupling is always generated
by radiative effects. For our numerical estimates, we include the
$a\rightarrow e^{+}e^{-}$ decay channel, taking $g_{aee}=0.37\times
10^{-3}~{}c_{GG}$ (see Eq. (62) of Bauer  _et al._ (2021b)), although the
effect of this coupling is almost invisible in our plots. Substantially
stronger coupling to electrons is possible if it arises at tree-level in the
UV or is induced by top-quarks Bauer  _et al._ (2021b); in these cases our
constraints would strengthen. In principle for $m_{a}\geq 2m_{\mu}$ the muon
decay channel can also be included, but this only affects the ceiling of our
constraint by an $\mathcal{O}(1)$ factor since
$\Gamma(a\rightarrow\gamma\gamma)$ is comparable to
$\Gamma(a\rightarrow\mu^{+}\mu^{-})$. If a detector can measure muon tracks,
then $a\rightarrow\mu^{+}\mu^{-}$ would present an additional signal channel
with much lower backgrounds. However, since the impact of muon couplings on
our results is relatively modest for the co-dominance and gluon dominance
scenarios, and this region can be probed by other experiments Co  _et al._
(2023); Acciarri  _et al._ (2023), we do not discuss it further.

MicroBooNE search for heavy scalars: Using the above formulae, we can recast a
recent search by MicroBooNE Abratenko  _et al._ (2021). Their search channel
was $K\rightarrow\pi h_{D}$ followed by $h_{D}\rightarrow e^{+}e^{-}$ with
$h_{D}$ a dark Higgs Batell  _et al._ (2019). If we consider the $a\rightarrow
e^{+}e^{-}$ search, then mapping their result to a heavy axion is an immediate
constraint on $g_{aee}$. Using the induced $g_{aee}$ from gluon couplings
mentioned above, we then obtain a limit on $f_{a}$. We find that the
constraints on $f_{a}$ obtained in this gluon-dominance scenario are very
weak, excluding $f_{a}\lesssim 1~{}{\rm TeV}$, which is already ruled out by
other experiments, and so we do not include this in our summary plot (see
Coloma  _et al._ (2022) for ALPs coupled exclusively to electroweak bosons).

For the $a\rightarrow\gamma\gamma$ channel, it is not possible to interpret
the MicroBooNE result as a constraint. The search in Abratenko  _et al._
(2021) made use of a boosted decision tree (BDT) in classifying their events
and should reject $\gamma\gamma$ topologies.333In principle there should be
some probability of a $\gamma\gamma$ pair contaminating the BDT-tagged
$e^{+}e^{-}$ sample, but this would require collaboration input. As a crude
estimate of the sensitivity, we may assume that a dedicated analysis for
$\gamma\gamma$ final states is performed with comparable BDT performance. Then
the sensitivity to $f_{a}$ (from the same dataset) would be given by equating
$\Gamma(K\to\pi h_{D})\times\Gamma(h_{D}\to ee)$, evaluated with the upper
bound on the mixing angle in Abratenko  _et al._ (2021), to our
$\Gamma(K\to\pi a)\times\Gamma(a\to\gamma\gamma)\times e^{-L/\lambda_{a}}$
with $L\simeq 100\,\mathrm{m}$ the distance between MicroBooNE and the NuMI
absorber (i.e. the KDAR source). This treatment is valid since the decay
length of the dark Higgs is much longer than $L$ for the mixing angles in
Abratenko  _et al._ (2021), and because the MicroBooNE detector is much
smaller than $L$.

In Figure 1 (gluon dominance) and Figure 2 (co-dominance), we plot the
sensitivity of MicroBooNE estimated in this way by the orange lines.
MicroBooNE may be able to explore certain small regions of parameter space not
covered by existing experiments if a dedicated search for $\gamma\gamma$ final
states is performed; this is qualitatively similar to the situation with a
dark Higgs Batell  _et al._ (2019); Abratenko  _et al._ (2021). As alluded to
above, searches for dimuon final states may also be of interest for $m_{a}\geq
2m_{\mu}$ since MicroBooNE can easily reconstruct $\mu^{+}\mu^{-}$ pairs.

Figure 2: (Co-Dominance) The same figure as Figure 1 but with
$c_{WW}=c_{BB}=c_{GG}$. The sensitivity for $m_{a}\ll m_{\pi}$ is worse than
the gluon dominance case since $|c_{\gamma\gamma}^{\mathrm{eff}}|\ll 1$ for
this specific choice of the parameters.

Searches at JSNS2: Next, we consider JSNS2 Ajimura  _et al._ (2017, 2021);
Maruyama (2022), which is designed to test the excess of events seen at LSND
Aguilar  _et al._ (2001). A crucial difference between JSNS2 and LSND is the
proton beam energy, $3~{}{\rm GeV}$ vs $0.8~{}{\rm GeV}$, such that JSNS2
serves as both a $\pi$DAR and KDAR facility, whereas LSND had much fewer KDAR
events (if any) Gudima  _et al._ (2010). We, therefore, find that JSNS2 offers
much more compelling sensitivity to heavy axions and do not consider LSND
further in what follows.

To estimate the sensitivity at JSNS2 we take the number of KDAR neutrinos per
proton on target ($N_{\nu_{\mu}}=0.0034/{\rm POT}$) from Ajimura  _et al._
(2017). Although the kaon production rate is uncertain ($N_{\nu_{\mu}}$ can be
twice as large Axani  _et al._ (2015)), as we mentioned in the Introduction,
in the actual experiment the KDAR neutrino flux can be measured _in situ_.
Using estimates from Fig. 3 of Jordan  _et al._ (2018), we find that the
search will be nearly background free $N_{\rm bkg}\leq 1$ in three years of
operation, except for $E_{a}\leq 238~{}{\rm MeV}$ (corresponding to
$m_{a}<104\,\mathrm{MeV}$) where $N_{\rm bkg}\approx 2.5$ over three years.
Ref. Jordan  _et al._ (2018) assumes an additional lead shielding of the
detector, which has not been put in place at JSNS2. Cosmic backgrounds are
therefore underestimated in that work. However, the axion signal we consider
benefits from the high energy signal region, $E>227\,\mathrm{MeV}$, where
cosmic events are suppressed. Further reductions in the cosmic muon rate can
be achieved using coincident tagging with daughter Michel electrons Jeon
(2021). Since studies of cosmic background mitigation are ongoing, we plot
projections using a 5-event contour. Backgrounds will be smooth above
$227~{}{\rm MeV}$, and a sideband analysis can be used to estimate their size
in situ such that a search will always be statistics rather than systematics
limited. As we discuss below, if backgrounds are high at JSNS2 a search could
be performed using their second detector Ajimura  _et al._ (2020); Maruyama
(2022).

To compute the expected number of axion events, we take $L=24~{}{\rm m}$ as
the distance between the KDAR source and the JSNS2 detector, and the detector
volume $V=(17~{}{\rm tonnes})/(0.852~{}{\rm g/ml})=20.35~{}{\rm m}^{3}$
Ajimura  _et al._ (2021); Kim  _et al._ (2011). We assume $10^{23}$ POT,
corresponding to roughly three years of live time. This results in $\sim
3\times 10^{20}$ stopped $K^{+}$ in total, which is much larger than the
$10^{12}\,\mathchar 45\relax\mathchar 45\relax\,10^{13}$ $K^{+}$ decay events
at NA62 Cortina Gil  _et al._ (2021).

In Figure 1 (gluon dominance) and Figure 2 (co-dominance), we plot the
sensitivity of JSNS2 by the blue lines, together with the existing constraints
and future sensitivities. The figures demonstrate that JSNS2 has excellent
sensitivity to heavy axions. Note that the event number scales as
$f_{a}^{-4}$. Therefore, even if we instead require e.g. 50 events, the
sensitivity only weakens by $10^{1/4}\simeq 1.8$, which does not alter our
main conclusion.

In our analysis, we have focused on the JSNS2 near detector. However, the
JSNS2 collaboration recently installed another detector at a far location
Ajimura  _et al._ (2020); Maruyama (2022), $48\,\mathrm{m}$ away from the
source. This second detector is expected to start taking data soon and
contains $32$ tonnes of liquid scintillator as a fiducial volume; the larger
volume partially compensating for the longer baseline. Together with possibly
smaller backgrounds due to its location, the far detector may be better suited
for axion searches once it starts its operation.

Conclusions: KDAR provides a clean smoking gun signature of hadronically
coupled axions. A $K^{+}$ production target, coupled with a large volume
detector placed $\sim 10\text{--}100~{}{\rm m}$ away, allows for a powerful
probe of visibly decaying particles lighter than the kaon (e.g. dark scalars
or heavy neutral leptons Batell  _et al._ (2019); Kelly and Machado (2021)).
In the context of KDAR, axions are particularly compelling due to their well-
motivated hadronic couplings, which are necessary in any model that addresses
the strong CP problem.

We have focused on two benchmark scenarios (gluon dominance and co-dominance)
for ease of comparison with the literature. It is interesting to understand
how constraints vary with model-dependent coupling textures. The visible
decays we consider here are governed both by ${\rm BR}(K\rightarrow\pi a)$ and
the axion decay length and scale as $1/f_{a}^{4}$. By way of contrast, the
constraints from NA62 depend only on ${\rm BR}(K\rightarrow\pi a)$ and scale
as $1/f_{a}^{2}$. Stronger hadronic couplings will therefore favor NA62 over
JSNS2. Conversely, weaker hadronic couplings and/or a larger
$c_{\gamma\gamma}^{\rm eff}$ will favor JSNS2 over NA62. Beam dump searches
scale the same way as JSNS2.

We find that JSNS2 will have world-leading sensitivity to heavy axions. One
might imagine a competitive experimental landscape of modern high-intensity
low-energy proton beams. Notably, we find that JSNS2 is likely to provide
unsurpassed sensitivity, since other facilities suffer from low $K^{+}$ yields
(e.g. at a PIP-II beam dump Toups  _et al._ (2022a), LANSCE Aguilar-Arevalo
_et al._ (2022), or the SNS Akimov  _et al._ (2020)) or do not have
competitive intensity (e.g. the SBN beam dump concept Toups  _et al._
(2022b)). Experiments with detectors far downstream also suffer from large
$1/L^{2}$ geometric suppressions c.f. Eq. 5. Modified experimental designs,
e.g. a PIP-II beam dump with a proton beam energy $T_{p}\gtrsim 2~{}{\rm
GeV}$, could allow for competitive KDAR rates. A large volume detector placed
near the DUNE hadron absorber or coupled with a high-intensity 8 GeV beam for
a muon collider demonstrator would both offer promising future sensitivity.
Nevertheless, with JSNS2 already taking data Maruyama (2022), there is an
immediate opportunity to shed light on heavy axion models in currently
unprobed regions of parameter space. We encourage the JSNS2 collaboration to
incorporate axion searches into their central physics program.

Acknowledgments: This project arose from discussions at the Aspen Center for
Physics (supported by National Science Foundation grant PHY-2210452) and we
would like to collectively thank the center for their hospitality and support.
We thank Josh Spitz and Takasumi Maruyama for useful correspondence and
feedback, and Bertrand Echenard for suggestions involving future muon
facilities. YE and ZL are supported in part by the DOE grant DE-SC0011842. RP
is supported by DOE Award Number DE-SC0011632 and by the Walter Burke
Institute for Theoretical Physics. RP’s research is supported by the Neutrino
Theory Network Program Grant under Award Number DEAC02-07CHI11359 and the US
DOE under Award Number DE-SC0020250.

## References

  * Graner  _et al._ (2016) B. Graner, Y. Chen, E. G. Lindahl, and B. R. Heckel, “Reduced Limit on the Permanent Electric Dipole Moment of Hg199,” Phys. Rev. Lett. 116, 161601 (2016), [Erratum: Phys.Rev.Lett. 119, 119901 (2017)], arXiv:1601.04339 [physics.atom-ph] .
  * Abel  _et al._ (2020) C. Abel _et al._ , “Measurement of the Permanent Electric Dipole Moment of the Neutron,” Phys. Rev. Lett. 124, 081803 (2020), arXiv:2001.11966 [hep-ex] .
  * Peccei and Quinn (1977a) R. D. Peccei and Helen R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38, 1440–1443 (1977a).
  * Peccei and Quinn (1977b) R. D. Peccei and Helen R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D 16, 1791–1797 (1977b).
  * Weinberg (1978) Steven Weinberg, “A New Light Boson?” Phys. Rev. Lett. 40, 223–226 (1978).
  * Wilczek (1978) Frank Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40, 279–282 (1978).
  * Kamionkowski and March-Russell (1992) Marc Kamionkowski and John March-Russell, “Planck scale physics and the Peccei-Quinn mechanism,” Phys. Lett. B 282, 137–141 (1992), arXiv:hep-th/9202003 .
  * Barr and Seckel (1992) Stephen M. Barr and D. Seckel, “Planck scale corrections to axion models,” Phys. Rev. D 46, 539–549 (1992).
  * Ghigna  _et al._ (1992) S. Ghigna, Maurizio Lusignoli, and M. Roncadelli, “Instability of the invisible axion,” Phys. Lett. B 283, 278–281 (1992).
  * Holman  _et al._ (1992) Richard Holman, Stephen D. H. Hsu, Thomas W. Kephart, Edward W. Kolb, Richard Watkins, and Lawrence M. Widrow, “Solutions to the strong CP problem in a world with gravity,” Phys. Lett. B 282, 132–136 (1992), arXiv:hep-ph/9203206 .
  * Dimopoulos (1979) Savas Dimopoulos, “A Solution of the Strong CP Problem in Models With Scalars,” Phys. Lett. B 84, 435–439 (1979).
  * Tye (1981) S. H. H. Tye, “A Superstrong Force With a Heavy Axion,” Phys. Rev. Lett. 47, 1035 (1981).
  * Holdom and Peskin (1982) Bob Holdom and Michael E. Peskin, “Raising the Axion Mass,” Nucl. Phys. B 208, 397–412 (1982).
  * Flynn and Randall (1987) Jonathan M. Flynn and Lisa Randall, “A Computation of the Small Instanton Contribution to the Axion Potential,” Nucl. Phys. B 293, 731–739 (1987).
  * Rubakov (1997) V. A. Rubakov, “Grand unification and heavy axion,” JETP Lett. 65, 621–624 (1997), arXiv:hep-ph/9703409 .
  * Berezhiani  _et al._ (2001) Zurab Berezhiani, Leonida Gianfagna, and Maurizio Giannotti, “Strong CP problem and mirror world: The Weinberg-Wilczek axion revisited,” Phys. Lett. B 500, 286–296 (2001), arXiv:hep-ph/0009290 .
  * Hook (2015) Anson Hook, “Anomalous solutions to the strong CP problem,” Phys. Rev. Lett. 114, 141801 (2015), arXiv:1411.3325 [hep-ph] .
  * Fukuda  _et al._ (2015) Hajime Fukuda, Keisuke Harigaya, Masahiro Ibe, and Tsutomu T. Yanagida, “Model of visible QCD axion,” Phys. Rev. D 92, 015021 (2015), arXiv:1504.06084 [hep-ph] .
  * Gherghetta  _et al._ (2016) Tony Gherghetta, Natsumi Nagata, and Mikhail Shifman, “A Visible QCD Axion from an Enlarged Color Group,” Phys. Rev. D 93, 115010 (2016), arXiv:1604.01127 [hep-ph] .
  * Dimopoulos  _et al._ (2016) Savas Dimopoulos, Anson Hook, Junwu Huang, and Gustavo Marques-Tavares, “A collider observable QCD axion,” JHEP 11, 052 (2016), arXiv:1606.03097 [hep-ph] .
  * Agrawal and Howe (2018a) Prateek Agrawal and Kiel Howe, “Factoring the Strong CP Problem,” JHEP 12, 029 (2018a), arXiv:1710.04213 [hep-ph] .
  * Agrawal and Howe (2018b) Prateek Agrawal and Kiel Howe, “A Flavorful Factoring of the Strong CP Problem,” JHEP 12, 035 (2018b), arXiv:1712.05803 [hep-ph] .
  * Gaillard  _et al._ (2018) M. K. Gaillard, M. B. Gavela, R. Houtz, P. Quilez, and R. Del Rey, “Color unified dynamical axion,” Eur. Phys. J. C 78, 972 (2018), arXiv:1805.06465 [hep-ph] .
  * Lillard and Tait (2018) Benjamin Lillard and Tim M. P. Tait, “A High Quality Composite Axion,” JHEP 11, 199 (2018), arXiv:1811.03089 [hep-ph] .
  * Fuentes-Martín  _et al._ (2019) Javier Fuentes-Martín, Mario Reig, and Avelino Vicente, “Strong $CP$ problem with low-energy emergent QCD: The 4321 case,” Phys. Rev. D 100, 115028 (2019), arXiv:1907.02550 [hep-ph] .
  * Csáki  _et al._ (2020) Csaba Csáki, Maximilian Ruhdorfer, and Yuri Shirman, “UV Sensitivity of the Axion Mass from Instantons in Partially Broken Gauge Groups,” JHEP 04, 031 (2020), arXiv:1912.02197 [hep-ph] .
  * Hook  _et al._ (2020) Anson Hook, Soubhik Kumar, Zhen Liu, and Raman Sundrum, “High Quality QCD Axion and the LHC,” Phys. Rev. Lett. 124, 221801 (2020), arXiv:1911.12364 [hep-ph] .
  * Gherghetta  _et al._ (2020) Tony Gherghetta, Valentin V. Khoze, Alex Pomarol, and Yuri Shirman, “The Axion Mass from 5D Small Instantons,” JHEP 03, 063 (2020), arXiv:2001.05610 [hep-ph] .
  * Gherghetta and Nguyen (2020) Tony Gherghetta and Minh D. Nguyen, “A Composite Higgs with a Heavy Composite Axion,” JHEP 12, 094 (2020), arXiv:2007.10875 [hep-ph] .
  * Valenti  _et al._ (2022) Alessandro Valenti, Luca Vecchi, and Ling-Xiao Xu, “Grand Color axion,” JHEP 10, 025 (2022), arXiv:2206.04077 [hep-ph] .
  * Kivel  _et al._ (2022) Alexey Kivel, Julien Laux, and Felix Yu, “Supersizing axions with small size instantons,” JHEP 11, 088 (2022), arXiv:2207.08740 [hep-ph] .
  * Bauer  _et al._ (2017) Martin Bauer, Matthias Neubert, and Andrea Thamm, “Collider Probes of Axion-Like Particles,” JHEP 12, 044 (2017), arXiv:1708.00443 [hep-ph] .
  * Irastorza and Redondo (2018) Igor G. Irastorza and Javier Redondo, “New experimental approaches in the search for axion-like particles,” Prog. Part. Nucl. Phys. 102, 89–159 (2018), arXiv:1801.08127 [hep-ph] .
  * Arvanitaki  _et al._ (2010) Asimina Arvanitaki, Savas Dimopoulos, Sergei Dubovsky, Nemanja Kaloper, and John March-Russell, “String Axiverse,” Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Jaeckel and Ringwald (2010) Joerg Jaeckel and Andreas Ringwald, “The Low-Energy Frontier of Particle Physics,” Ann. Rev. Nucl. Part. Sci. 60, 405–437 (2010), arXiv:1002.0329 [hep-ph] .
  * Essig  _et al._ (2010) Rouven Essig, Roni Harnik, Jared Kaplan, and Natalia Toro, “Discovering New Light States at Neutrino Experiments,” Phys. Rev. D 82, 113008 (2010), arXiv:1008.0636 [hep-ph] .
  * Döbrich  _et al._ (2016) 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), arXiv:1512.03069 [hep-ph] .
  * Dolan  _et al._ (2017) 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), [Erratum: JHEP 03, 190 (2021)], arXiv:1709.00009 [hep-ph] .
  * Döbrich  _et al._ (2019) 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), [Erratum: JHEP 10, 046 (2020)], arXiv:1904.02091 [hep-ph] .
  * Harland-Lang  _et al._ (2019) Lucian Harland-Lang, Joerg Jaeckel, and Michael Spannowsky, “A fresh look at ALP searches in fixed target experiments,” Phys. Lett. B 793, 281–289 (2019), arXiv:1902.04878 [hep-ph] .
  * Dent  _et al._ (2020) James B. Dent, Bhaskar Dutta, Doojin Kim, Shu Liao, Rupak Mahapatra, Kuver Sinha, and Adrian Thompson, “New Directions for Axion Searches via Scattering at Reactor Neutrino Experiments,” Phys. Rev. Lett. 124, 211804 (2020), arXiv:1912.05733 [hep-ph] .
  * Brdar  _et al._ (2021) Vedran Brdar, Bhaskar Dutta, Wooyoung Jang, Doojin Kim, Ian M. Shoemaker, Zahra Tabrizi, Adrian Thompson, and Jaehoon Yu, “Axionlike Particles at Future Neutrino Experiments: Closing the Cosmological Triangle,” Phys. Rev. Lett. 126, 201801 (2021), arXiv:2011.07054 [hep-ph] .
  * Jaeckel  _et al._ (2013) Joerg Jaeckel, Martin Jankowiak, and Michael Spannowsky, “LHC probes the hidden sector,” Phys. Dark Univ. 2, 111–117 (2013), arXiv:1212.3620 [hep-ph] .
  * Mimasu and Sanz (2015) Ken Mimasu and Verónica Sanz, “ALPs at Colliders,” JHEP 06, 173 (2015), arXiv:1409.4792 [hep-ph] .
  * Jaeckel and Spannowsky (2016) Joerg Jaeckel and Michael Spannowsky, “Probing MeV to 90 GeV axion-like particles with LEP and LHC,” Phys. Lett. B 753, 482–487 (2016), arXiv:1509.00476 [hep-ph] .
  * Izaguirre  _et al._ (2017) Eder Izaguirre, Tongyan Lin, and Brian Shuve, “Searching for Axionlike Particles in Flavor-Changing Neutral Current Processes,” Phys. Rev. Lett. 118, 111802 (2017), arXiv:1611.09355 [hep-ph] .
  * Knapen  _et al._ (2017) Simon Knapen, Tongyan Lin, Hou Keong Lou, and Tom Melia, “Searching for Axionlike Particles with Ultraperipheral Heavy-Ion Collisions,” Phys. Rev. Lett. 118, 171801 (2017), arXiv:1607.06083 [hep-ph] .
  * Brivio  _et al._ (2017) I. Brivio, M. B. Gavela, L. Merlo, K. Mimasu, J. M. No, R. del Rey, and V. Sanz, “ALPs Effective Field Theory and Collider Signatures,” Eur. Phys. J. C 77, 572 (2017), arXiv:1701.05379 [hep-ph] .
  * Mariotti  _et al._ (2018) Alberto Mariotti, Diego Redigolo, Filippo Sala, and Kohsaku Tobioka, “New LHC bound on low-mass diphoton resonances,” Phys. Lett. B 783, 13–18 (2018), arXiv:1710.01743 [hep-ph] .
  * Cid Vidal  _et al._ (2019) Xabier Cid Vidal, Alberto Mariotti, Diego Redigolo, Filippo Sala, and Kohsaku Tobioka, “New Axion Searches at Flavor Factories,” JHEP 01, 113 (2019), [Erratum: JHEP 06, 141 (2020)], arXiv:1810.09452 [hep-ph] .
  * Beacham  _et al._ (2020) J. Beacham _et al._ , “Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report,” J. Phys. G 47, 010501 (2020), arXiv:1901.09966 [hep-ex] .
  * Alonso-Álvarez  _et al._ (2019) G. Alonso-Álvarez, M. B. Gavela, and P. Quilez, “Axion couplings to electroweak gauge bosons,” Eur. Phys. J. C 79, 223 (2019), arXiv:1811.05466 [hep-ph] .
  * Ebadi  _et al._ (2019) Javad Ebadi, Sara Khatibi, and Mojtaba Mohammadi Najafabadi, “New probes for axionlike particles at hadron colliders,” Phys. Rev. D 100, 015016 (2019), arXiv:1901.03061 [hep-ph] .
  * Gavela  _et al._ (2020) M. B. Gavela, J. M. No, V. Sanz, and J. F. de Trocóniz, “Nonresonant Searches for Axionlike Particles at the LHC,” Phys. Rev. Lett. 124, 051802 (2020), arXiv:1905.12953 [hep-ph] .
  * Altmannshofer  _et al._ (2020) Wolfgang Altmannshofer, Stefania Gori, and Dean J. Robinson, “Constraining axionlike particles from rare pion decays,” Phys. Rev. D 101, 075002 (2020), arXiv:1909.00005 [hep-ph] .
  * Gershtein  _et al._ (2021) Yuri Gershtein, Simon Knapen, and Diego Redigolo, “Probing naturally light singlets with a displaced vertex trigger,” Phys. Lett. B 823, 136758 (2021), arXiv:2012.07864 [hep-ph] .
  * Di Luzio  _et al._ (2020) Luca Di Luzio, Maurizio Giannotti, Enrico Nardi, and Luca Visinelli, “The landscape of QCD axion models,” Phys. Rept. 870, 1–117 (2020), arXiv:2003.01100 [hep-ph] .
  * Knapen  _et al._ (2022) Simon Knapen, Soubhik Kumar, and Diego Redigolo, “Searching for axionlike particles with data scouting at ATLAS and CMS,” Phys. Rev. D 105, 115012 (2022), arXiv:2112.07720 [hep-ph] .
  * Co  _et al._ (2023) Raymond T. Co, Soubhik Kumar, and Zhen Liu, “Searches for heavy QCD axions via dimuon final states,” JHEP 02, 111 (2023), arXiv:2210.02462 [hep-ph] .
  * Acciarri  _et al._ (2023) R. Acciarri _et al._ (ArgoNeuT), “First Constraints on Heavy QCD Axions with a Liquid Argon Time Projection Chamber Using the ArgoNeuT Experiment,” Phys. Rev. Lett. 130, 221802 (2023), arXiv:2207.08448 [hep-ex] .
  * Coloma  _et al._ (2022) Pilar Coloma, Pilar Hernández, and Salvador Urrea, “New bounds on axion-like particles from MicroBooNE,” JHEP 08, 025 (2022), arXiv:2202.03447 [hep-ph] .
  * Batell  _et al._ (2019) Brian Batell, Joshua Berger, and Ahmed Ismail, “Probing the Higgs Portal at the Fermilab Short-Baseline Neutrino Experiments,” Phys. Rev. D 100, 115039 (2019), arXiv:1909.11670 [hep-ph] .
  * Spitz (2014) J. Spitz, “Cross Section Measurements with Monoenergetic Muon Neutrinos,” Phys. Rev. D 89, 073007 (2014), arXiv:1402.2284 [physics.ins-det] .
  * Nikolakopoulos  _et al._ (2021) A. Nikolakopoulos, V. Pandey, J. Spitz, and Natalie Jachowicz, “Modeling quasielastic interactions of monoenergetic kaon decay-at-rest neutrinos,” Phys. Rev. C 103, 064603 (2021), arXiv:2010.05794 [nucl-th] .
  * Kelly  _et al._ (2021) Kevin J. Kelly, Soubhik Kumar, and Zhen Liu, “Heavy axion opportunities at the DUNE near detector,” Phys. Rev. D 103, 095002 (2021), arXiv:2011.05995 [hep-ph] .
  * Jerhot  _et al._ (2022) 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), arXiv:2201.05170 [hep-ph] .
  * Afik  _et al._ (2023) Yoav Afik, Babette Döbrich, Jan Jerhot, Yotam Soreq, and Kohsaku Tobioka, “Probing Long-lived Axions at the KOTO Experiment,” (2023), arXiv:2303.01521 [hep-ph] .
  * Bauer  _et al._ (2021a) Martin Bauer, Matthias Neubert, Sophie Renner, Marvin Schnubel, and Andrea Thamm, “Consistent Treatment of Axions in the Weak Chiral Lagrangian,” Phys. Rev. Lett. 127, 081803 (2021a), arXiv:2102.13112 [hep-ph] .
  * Abratenko  _et al._ (2021) P. Abratenko _et al._ (MicroBooNE), “Search for a Higgs Portal Scalar Decaying to Electron-Positron Pairs in the MicroBooNE Detector,” Phys. Rev. Lett. 127, 151803 (2021), arXiv:2106.00568 [hep-ex] .
  * Chang  _et al._ (2018) Jae Hyeok Chang, Rouven Essig, and Samuel D. McDermott, “Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle,” JHEP 09, 051 (2018), arXiv:1803.00993 [hep-ph] .
  * Ertas and Kahlhoefer (2020) Fatih Ertas and Felix Kahlhoefer, “On the interplay between astrophysical and laboratory probes of MeV-scale axion-like particles,” JHEP 07, 050 (2020), arXiv:2004.01193 [hep-ph] .
  * Depta  _et al._ (2020) Paul Frederik Depta, Marco Hufnagel, and Kai Schmidt-Hoberg, “Robust cosmological constraints on axion-like particles,” JCAP 05, 009 (2020), arXiv:2002.08370 [hep-ph] .
  * Artamonov  _et al._ (2009) A. V. Artamonov _et al._ (BNL-E949), “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), arXiv:0903.0030 [hep-ex] .
  * Cortina Gil  _et al._ (2021) Eduardo Cortina Gil _et al._ (NA62), “Measurement of the very rare K+→${\pi}^{+}\nu\overline{\nu}$ decay,” JHEP 06, 093 (2021), arXiv:2103.15389 [hep-ex] .
  * Bergsma  _et al._ (1985) F. Bergsma _et al._ (CHARM), “Search for Axion Like Particle Production in 400-GeV Proton - Copper Interactions,” Phys. Lett. B 157, 458–462 (1985).
  * Blumlein  _et al._ (1991) 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).
  * Aidala  _et al._ (2019) C. A. Aidala _et al._ (SeaQuest), “The SeaQuest Spectrometer at Fermilab,” Nucl. Instrum. Meth. A 930, 49–63 (2019), arXiv:1706.09990 [physics.ins-det] .
  * Berlin  _et al._ (2018) Asher Berlin, Stefania Gori, Philip Schuster, and Natalia Toro, “Dark Sectors at the Fermilab SeaQuest Experiment,” Phys. Rev. D 98, 035011 (2018), arXiv:1804.00661 [hep-ph] .
  * Blinov  _et al._ (2022) Nikita Blinov, Elizabeth Kowalczyk, and Margaret Wynne, “Axion-like particle searches at DarkQuest,” JHEP 02, 036 (2022), arXiv:2112.09814 [hep-ph] .
  * Ariga  _et al._ (2019) Akitaka Ariga _et al._ (FASER), “FASER’s physics reach for long-lived particles,” Phys. Rev. D 99, 095011 (2019), arXiv:1811.12522 [hep-ph] .
  * Yamanaka (2012) Taku Yamanaka (KOTO), “The J-PARC KOTO experiment,” PTEP 2012, 02B006 (2012).
  * Anelli  _et al._ (2015) M. Anelli _et al._ (SHiP), “A facility to Search for Hidden Particles (SHiP) at the CERN SPS,” (2015), arXiv:1504.04956 [physics.ins-det] .
  * Bauer  _et al._ (2021b) Martin Bauer, Matthias Neubert, Sophie Renner, Marvin Schnubel, and Andrea Thamm, “The Low-Energy Effective Theory of Axions and ALPs,” JHEP 04, 063 (2021b), arXiv:2012.12272 [hep-ph] .
  * Ajimura  _et al._ (2017) S. Ajimura _et al._ , “Technical Design Report (TDR): Searching for a Sterile Neutrino at J-PARC MLF (E56, JSNS2),” (2017), arXiv:1705.08629 [physics.ins-det] .
  * Ajimura  _et al._ (2021) S. Ajimura _et al._ (JSNS2), “The JSNS2 detector,” Nucl. Instrum. Meth. A 1014, 165742 (2021), arXiv:2104.13169 [physics.ins-det] .
  * Maruyama (2022) Takasumi Maruyama (JSNS2, JSNS2 -II,), “The status of JSNS2 and JSNS2-II,” PoS NuFact2021, 159 (2022).
  * Aguilar  _et al._ (2001) A. Aguilar _et al._ (LSND), “Evidence for neutrino oscillations from the observation of $\bar{\nu}_{e}$ appearance in a $\bar{\nu}_{\mu}$ beam,” Phys. Rev. D 64, 112007 (2001), arXiv:hep-ex/0104049 .
  * Gudima  _et al._ (2010) Konstantin K. Gudima, N. V. Mokhov, and S. I. Striganov, “Kaon Yields for 2 to 8 GeV Proton Beams,” in  _Workshop on Applications of High Intensity Proton Accelerators_ (2010) pp. 115–119.
  * Axani  _et al._ (2015) S Axani, G Collin, JM Conrad, MH Shaevitz, J Spitz, and T Wongjirad, “Decisive disappearance search at high $\Delta m^{2}$ with monoenergetic muon neutrinos,” Phys. Rev. D 92, 092010 (2015), arXiv:1506.05811 [physics.ins-det] .
  * Jordan  _et al._ (2018) Johnathon R. Jordan, Yonatan Kahn, Gordan Krnjaic, Matthew Moschella, and Joshua Spitz, “Signatures of Pseudo-Dirac Dark Matter at High-Intensity Neutrino Experiments,” Phys. Rev. D 98, 075020 (2018), arXiv:1806.05185 [hep-ph] .
  * Jeon (2021) HyoungKu Jeon, “Cosmic ray induced Background study at the JSNS2 experiment,”  (2021), for the JSNS2 collaboration.
  * Ajimura  _et al._ (2020) S. Ajimura _et al._ , “Proposal: JSNS2-II,” (2020), arXiv:2012.10807 [hep-ex] .
  * Kim  _et al._ (2011) Byoung Chan Kim, Seung Won Seo, and Kyung Kwang Joo, “Measurement of the density of liquid scintillator solvents for neutrino experiments,” (2011).
  * Kelly and Machado (2021) Kevin James Kelly and Pedro A. N. Machado, “MicroBooNE experiment, NuMI absorber, and heavy neutral leptons,” Phys. Rev. D 104, 055015 (2021), arXiv:2106.06548 [hep-ph] .
  * Toups  _et al._ (2022a) M. Toups _et al._ , “PIP2-BD: GeV Proton Beam Dump at Fermilab’s PIP-II Linac,” in  _Snowmass 2021_ (2022) arXiv:2203.08079 [hep-ex] .
  * Aguilar-Arevalo  _et al._ (2022) A. A. Aguilar-Arevalo _et al._ (CCM), “First Leptophobic Dark Matter Search from the Coherent–CAPTAIN-Mills Liquid Argon Detector,” Phys. Rev. Lett. 129, 021801 (2022), arXiv:2109.14146 [hep-ex] .
  * Akimov  _et al._ (2020) D. Akimov _et al._ (COHERENT), “Sensitivity of the COHERENT Experiment to Accelerator-Produced Dark Matter,” Phys. Rev. D 102, 052007 (2020), arXiv:1911.06422 [hep-ex] .
  * Toups  _et al._ (2022b) Matt Toups _et al._ , “SBN-BD: $\mathcal{O}$(10 GeV) Proton Beam Dump at Fermilab’s PIP-II Linac,” in  _Snowmass 2021_ (2022) arXiv:2203.08102 [hep-ex] .


# INTERPRETABLE DYSARTHRIC SPEAKER ADAPTATION BASED ON OPTIMAL-TRANSPORT

###### Abstract

This work addresses the mismatch problem between the distribution of training
data (source) and testing data (target), in the challenging context of
dysarthric speech recognition. We focus on Speaker Adaptation (SA) in command
speech recognition, where data from multiple sources (i.e., multiple speakers)
are available. Specifically, we propose an unsupervised Multi-Source Domain
Adaptation (MSDA) algorithm based on optimal-transport, called MSDA via
Weighted Joint Optimal Transport (MSDA-WJDOT). We achieve a Command Error Rate
relative reduction of 16$\%$ and 7$\%$ over the speaker-independent model and
the best competitor method, respectively. The strength of the proposed
approach is that, differently from any other existing SA method, it offers an
interpretable model that can also be exploited, in this context, to diagnose
dysarthria without any specific training. Indeed, it provides a closeness
measure between the target and the source speakers, reflecting their
similarity in terms of speech characteristics. Based on the similarity between
the target speaker and the healthy/dysarthric source speakers, we then define
the healthy/dysarthric score of the target speaker that we leverage to perform
dysarthria detection. This approach does not require any additional training
and achieves a 95$\%$ accuracy in the dysarthria diagnosis.

Index Terms: unsupervised speaker adaptation, dysarthric speech, optimal
transport, dysarthria detection

##  1 Introduction

Many Automatic Speech Recognition (ASR) models assume that test and training
datasets are sampled from the same distribution. However, new data can exhibit
a distribution change and drastically decrease the recognition performance.
This is due, for example, to different characteristics between speakers (e.g.,
different accents, different speaking style). Such a mismatch is even more
evident in speech-impaired people as the speech characteristics also depend on
the type and the severity of the disorder. This problem can be addressed
through supervised speaker adaptation (SA), where transcribed utterances of
the target speaker are used to adapt a speaker-independent (SI) model to
better recognize the speech of the target speaker.

Here we consider the more frequent unsupervised SA scenario, where the
transcriptions of the adaptation utterances are not available. Specifically,
we focus on unsupervised SA for spoken command recognition of dysarthric
speech. Dysarthria is a motor speech disorder consisting in the disruption of
the normal control of the vocal tract musculature, common in conditions such
as Parkinson Disease (PD), Amyotrophic Lateral Sclerosis (ALS) and post-stroke
motor impairments [1, 2]. Dysarthric patients tend to have mumbled speech,
acceleration or deceleration in the speaking rate, abnormal pitch and rhythm.
However, dysarthric speech is characterized by high variability as the
impairment varies based on the location and severity of the neurological
damage. This, beside the limited amount of available data, makes dysarthric
speech recognition very challenging. Indeed, although many studies improved
the ASR performance on dysarthric speech, they still poorly perform on low
intelligible speech [3, 4, 5, 6]. That points out the necessity of ad hoc SA
techniques for dysarthric speakers.

In [7], authors suggest that different adaptation techniques should be adopted
based on the impairment severity of the speaker. They employ two large
datasets to train two different SI models. The first dataset only includes
healthy speakers, while the second dataset contains additional recordings from
dysarthric speakers. Adaptation is then performed to better recognize the
utterances from a dysarthric corpus. Results show that adaptation methods for
dysarthric patients with mild or low dysarthria are more effective when the
initial model is trained on healthy speech, whereas the adaptation of severe
dysarthric patients improves when the SI model is trained on impaired-speech.
Following this direction, [8] introduces a two-stage adaptation approach in
which a SI model is first adapted to the general speaking style of multiple
dysarthric speakers, and then the adapted model is further adapted to the
speaker of interest. However, the described techniques require a large amount
of data for fine-tuning, that is usually not available. To overcome this
problem, [9] proposes to only fine-tune a subset of the network layers to
better adapt an ASR model to the dysarthric speech. Differently, [10] adopts
confusion matrix transducers within the weighted finite state transducer
decoding paradigm, to handle phonetic-articulatory errors made by the target
speakers. However, all the aforementioned studies are limited to the
supervised framework, whereas the unsupervised speaker adaptation is the most
common scenario.

We adopt an unsupervised SA approach based on Optimal Transport (OT) [11, 12],
named Multi-Source Domain Adaptation Weighted Joint Distribution Optimal
Transport (MSDA-WJDOT), firstly introduced in [13] for other applications. In
this work, we assume to have access to spoken commands of several speakers
(called source) and aim at finding a classification function for spoken
commands of a new speaker (called target) through a two-step procedure. The
first step learns a global feature extractor from a SI model, e.g., the last
hidden layer of a neural network. The second step learns a target-specific
classifier, from the feature space to the output space (i.e., the commands to
classify). Specifically, we define i) a weighted combination (with weights
$\boldsymbol{\alpha}$) of the joint probability distributions of the source
speakers on the feature-output product space; ii) a proxy joint probability
distribution of the target speaker on the feature-output product space, in
which the output space is represented by the predictions of the target
classifier. The algorithm is trained to jointly look for the optimal weights
$\boldsymbol{\alpha}$ and the optimal target classifier that minimize the
Wasserstein distance between the distributions defined in i) and ii).
Intuitively, at each iteration, this method simultaneously selects the sources
most similar to the target (by attributing them higher weights) and updates
the target classifier to increase this similarity. The advantage of MSDA-WJDOT
over the other adaptation approaches is that only the relevant source speakers
are used, avoiding redundant or misleading information. Mostly important, the
proposed approach offers an interpretable model as we found the source weights
$\boldsymbol{\alpha}$ reflecting a similarity between the source speakers and
the target speaker in terms of speech characteristics. Contrary to [7] in
which the sources are manually chosen based on the impairment severity of the
target speaker, MSDA-WJDOT is able to perform this selection automatically
considering also other important similarity aspects (e.g., the gender).

##  2 MSDA-WJDOT

In this section, we recall MSDA-WJDOT, firstly proposed in [13]. In the
classical machine learning problem, the aim is to learn a feature extractor
and a classifier, such that the combination of these two functions describes
the input-output relationship. MSDA-WJDOT assumes to have access to the
feature extraction function $g$ from the input space to an embedding space
$\mathcal{G}$ and looks for a target classifier defined on $\mathcal{G}$. This
is the case in which, for example, a pre-trained network is used to extract
meaningful features. Once $g$ is known, MSDA-WJDOT performs adaptation on the
feature-output product space. However, in the considered framework the feature
extractor is not available and, hence, the proposed method becomes a two-step
procedure, as described in the following.

###  2.1 Global feature extractor

The first step consists in estimating a global feature extractor, e.g., a
function that extracts features from the audio of a general speaker. To do
that, we learn a SI model, composed by a feature extraction function $g$ and a
classifier $f_{SI}$. Let us suppose to have $N_{j}$ utterances $(X_{j},Y_{j})$
of $J$ source speakers, representing the utterance acoustic features $X_{j}$
and the corresponding command label $Y_{j}$. We can train the model on the
source speakers $\\{(X_{j},Y_{j})\\}_{j=1}^{J}$ by

$\operatorname*{min}_{g,f_{SI}}\quad\sum_{j=1}^{J}\frac{1}{N_{j}}\sum_{i=1}^{N_{j}}\mathcal{L}(f_{SI}\circ
g(x_{j}^{i}),y_{j}^{i}),$ (1)

where $\mathcal{L}$ is the cross-entropy loss function. Note that in the
second step, $f_{SI}$ will not be used.

###  2.2 Adaptation

Once $g$ is given, we can define the joint distributions $p_{S_{j}}$, for each
source speaker $j$, with support on the product space
$\mathcal{G}\times\mathcal{Y}$, where $\mathcal{Y}$ is the label space. We
then define a convex combination of the source distributions

$p_{S}^{\alpha}=\sum_{j=1}^{J}\alpha_{j}p_{S_{j}}$

with $\boldsymbol{\alpha}$ element of the simplex $\Delta^{J}$ (i.e.,
$0\leq\alpha_{j}\leq 1$ for each $j$ and $\sum_{j=1}^{J}\alpha_{j}=1$).
Finally, we define the proxy target joint distribution $p_{T}^{f}$, with
support on $\mathcal{G}\times f(\mathcal{G})$, where $f$ is the target
classifier to learn. Here, the upper index $f$ indicates that the joint
distribution is over predicted labels and not reference labels. MSDA-WJDOT
aims at finding the $\boldsymbol{\alpha}$ weights and the classifier $f$ that
align the source combination of distributions $p_{S}^{\alpha}$ and the proxy
target distribution $p_{T}^{f}$ . This can be expressed as

$\operatorname*{min}_{\boldsymbol{\alpha},f}\quad
W_{D}\left(p_{T}^{f},p_{S}^{\alpha}\right)=\operatorname*{min}_{\boldsymbol{\alpha},f}\quad
W_{D}\left(p_{T}^{f},\sum_{j=1}^{J}\alpha_{j}p_{S,j}\right),$ (2)

where $W_{D}$ is the Wasserstein distance [11, 14]. Roughly speaking, the
source joint distributions can be seen as proxy joint distributions with
support on $\mathcal{G}\times f_{j}(\mathcal{G})$, where $f_{j}$ is the true
classifier (i.e., $f_{j}(X_{j})\equiv Y_{j}$). This approach dynamically
selects the most similar sources and re-adapts $f$ by “imitating” the best
$f_{j}$.

###  2.3 The role of $\boldsymbol{\alpha}$ and the similarity

We can consider $\boldsymbol{\alpha}$ as a measure of similarity between the
sources and the target, in terms of joint probability distribution. Indeed, in
order to minimize Eq. 2, the algorithm will attribute higher weights to the
source distributions closest to the proxy target distribution, and smaller (or
zero) weights to far source distributions. The zero-weighted distributions
will not contribute to the adaptation procedure. Consequently, misleading and
unrelated source domains are in practice not used and this allows to learn a
more accurate classifier. For example, if $\boldsymbol{\alpha}$ were a one-hot
vector with $\alpha_{j}=1$ , the method would consist in estimating $f$ such
that $p_{T}^{f}\simeq p_{S_{j}}$. Authors in [13] only offer a mathematical
explanation of the weight alpha. In this study, we explore the interpretation
of alpha in terms of speaker similarity and speech characteristics.

##  3 Experimental setup

###  3.1 Dataset

We employ the AllSpeak dataset [15], that consists of speech recordings from
29 Italian native speakers. Seventeen of these (thirteen males, four females)
are affected by Amyotrophic Lateral Sclerosis, while the remaining twelve (six
males, six females) are healthy control speakers. The dataset contains 25
commands in Italian, relative to basic needs such as “I am thirsty”, “I am
hungry”. This dataset is very challenging due to the small amount of
recordings. Indeed, only 2387 and 1857 examples have been recorded from
control and dysarthric speakers, respectively. Specifically, each command has
been recorded from 8 to 10 times by healthy speakers, while from 4 to 10 times
by patients – depending on the patient medical condition. In the following we
will adopt a speaker coding according to which F/M stands for female/male, and
the letter C is added to refer to control subjects.

We perform speaker adaptation of each dysarthric speaker (target speaker) by
using all the remaining speakers as source speakers. We trained the SI model
on the source speakers, by splitting the dataset into training and validation
sets. Specifically, the validation set consists of one example of each command
of each speaker and it has been used to perform early stopping. The unlabelled
target speaker data is split into adaptation set (80$\%$) and testing (20$\%$)
set, which contains one example of each command. To train MSDA-WJDOT, we
simultaneously employ the source training dataset and the adaptation dataset
of the target speaker.

###  3.2 Model details

The feature extraction function $g$ is represented by a BLSTM with five hidden
layers, each containing 250 memory blocks. Both the SI classifier $f_{SI}$ and
the target classifier $f$ are represented by a softmax layer. All weights are
initialized with Xavier initialization. Both SI model training and adaptation
are performed with Adam optimizer with 0.9 momentum and $\epsilon=e^{-8}$.
Learning rate is fixed to 0.001.

##  4 Results

###  4.1 Competitor models

We compare the proposed approach with an unadapted speaker- independent (SI)
model, in which the target speaker is directly tested on a model trained on
all source speakers. Notice that the literature lacks unsupervised MSDA
specific for dysarthric speaker adaptation. Hence, we implemented two MSDA-
WJDOT approaches - CJDOT and MJDOT (details can be found in [13]) - which are
different extensions of a state-of-the-art unsupervised domain adaptation
method - JDOT [16] \- adapted to the multi-speaker adaptation case. In
addition, we trained a supervised SA (SSA) model, in which we added a feed-
forward linear layer atop the input to the SI model and trained it on the
target speaker [17]. As this method uses the target labels, it must be
considered as the lower bound of the unsupervised SA performance error.

###  4.2 Command speech recognition

Table 1: Command Error Rate (CER) for each dysarthric target speaker provided
by the SI, SSA, MSDA-WJDOT and the competitor models.

Speaker | SI |  CJDOT [16] |  MJDOT [16] | MSDA-WJDOT |  SSA [17]
---|---|---|---|---|---
M01 | 35.79 | 32.77 | 31.93 | 31.93 | 25.00
M02 | 34.26 | 36.75 | 36.75 | 37.71 | 54.17
F01 | 63.16 | 52.63 | 49.12 | 49.12 | 60.00
F02 | 48.50 | 40.00 | 40.00 | 40.00 | 36.00
M03 | 64.44 | 68.89 | 71.11 | 57.78 | 55.56
M04 | 30.00 | 31.20 | 32.00 | 30.40 | 24.00
M05 | 18.62 | 17.46 | 15.08 | 14.29 | 17.39
F03 | 68.33 | 61.74 | 70.43 | 62.61 | 56.52
M06 | 48.67 | 34.78 | 33.91 | 35.65 | 26.09
M07 | 11.00 | 7.20 | 8.00 | 8.80 | 20.00
M08 | 39.50 | 36.00 | 41.60 | 33.60 | 32.00
M09 | 24.79 | 16.81 | 18.49 | 19.33 | 13.04
F04 | 48.07 | 38.60 | 38.60 | 38.60 | 40.91
M10 | 18.00 | 12.80 | 12.00 | 12.80 | 16.00
M11 | 56.50 | 47.20 | 48.80 | 45.60 | 40.00
M12 | 7.50 | 5.60 | 7.20 | 4.80 | 4.00
M13 | 30.91 | 45.45 | 43.64 | 21.82 | 9.09
A. CER | 38.11 | 34.46 | 34.37 | 32.05 | 31.16

Table 1 reports the results in terms of Command Error Rate (CER). A first
remark is that, although the SI model always achieved a CER between $15\%$ and
$20\%$ on the source validation set, it often had low accuracy on the target
speaker. Once again, this emphasizes the difficulty of an ASR system to
generalize to a new dysarthric speaker and the importance of the speaker
adaptation in this context.
The unsupervised SA carried out by MSDA- WJDOT outperforms all the methods by
providing the best Average CER. Indeed, it reduces the CER of $16\%$ and $7\%$
over the SI system and the best competitor model, respectively. Surprisingly,
MSDA-WDJOT achieves an Average CER similar to the SSA approach, in which the
labels are used.

It is essential to recall that the proposed approach provides a measure of
similarity between the source speakers and the target speaker, in terms of
probability distribution. We found that the recovered $\boldsymbol{\alpha}$
always attributes highest similarity scores to dysarthric speakers rather than
healthy ones, when the target speaker is dysarthric, and vice versa for
healthy target speakers. This means that $\boldsymbol{\alpha}$ weights can
catch the similarity in terms of speech characteristics and that this approach
can realistically estimate speakers similarity. Fig. 1 shows the recovered
weights for a dysarthric target speaker (M13) during the training.
Specifically, we report on the $x$-axis the epochs of training, while the
$y$-axis lists the healthy (in green) and dysarthric (in red) source speakers.
As we can see, after few epochs $\alpha_{j}$ becomes sparse, using in practice
only the data of the most similar speakers. Interestingly, all the healthy
speakers have zero weights whereas the $\boldsymbol{\alpha}$ mass is spread
along dysarthric speakers. Also, the highest weights are attributed to M03 and
M08, which are both males, while a lower importance is attributed to F01. This
suggests that MSDA-WJDOT also incorporates similarity information about the
gender.

Figure 1: $\boldsymbol{\alpha}$ coefficients recovered during MSDA-WJDOT
training for dysarthric target speaker M13. The $\alpha_{j}$ coefficients are
close to zero for healthy speakers (in green), while the highest weights are
attributed to dysarthric speakers (in red). Figure 2: HS and DS computed for
healthy (in green) and dysarthric (in red) speakers.

###  4.3 The $\boldsymbol{\alpha}$ weight and the dysarthria detection

We here propose a completely novel approach to perform dysarthia detection
based on the $\boldsymbol{\alpha}$ weights. As shown in Fig. 1, MSDA-WJDOT
associates speakers with similar voice characteristics. We investigated the
possibility of leveraging the $\boldsymbol{\alpha}$ weights of the command
classifier to infer if the target speaker is affected by dysarthria.
Specifically, we propose to classify a speaker as healthy or dysarthric based
on his/her similarity with the other subjects.

Let define $I_{c}$ as the set indexing the control speakers and $I_{d}$ as the
set of indices related to dysarthric speakers. We then define the Healthy
Score (HS) and the Dysarthric Score (DS) as follow

$HS=\sum_{j\in I_{c}}\alpha_{j},\quad DS=\sum_{j\in I_{d}}\alpha_{j}.$ (3)

As $\boldsymbol{\alpha}\in\Delta^{J}$, we have that $HS,DS\in[0,1]$ such that
$HS+DS=1$. Fig. 2 reports the computed scores for all dysarthric speakers and
for 5 control speakers. The control subjects always show a higher $HS$, while
for the patients (except for F02), we have a higher $DS$. Moving a step
forward, we can use these scores to perform dysarthria detection by stating
that

A speaker is affected by dysarthria if $DS>HS$.

This results in a final accuracy of 95$\%$. It is crucial to emphasize that,
differently from standard dysarthria detection approaches [18, 19, 20], any
specific training has been performed: we directly infer the speaker health
status from the similarity coefficients $\boldsymbol{\alpha}$ that provide us
an interpretable model. Note that none of existing SA methods are able to also
perform dysarthria detection. Hence, we cannot fairly compare our approach to
any competitor model.

##  5 Conclusions

In this work, we addressed usupervised speaker adaptation in the challenging
context of pathological speech by employing MSDA-WJDOT algorithm. Our approach
significantly improves the speaker-independent system, reducing the CER of
$16\%$, and it also outperforms all the competitor models. The advantage of
the proposed method is that it does not require the labels for the target
speaker and, interestingly, it also provides source-target probability
distribution similarity coefficients $\boldsymbol{\alpha}$, reflecting the
source-target speakers characteristics similarity. From these coefficients, we
derived the Healthy (HI) and Dysarthric (DI) Index of a target speaker that
allowed us to detect dysarthria with an accuracy of 95$\%$.

The novelty of our work relies on several aspects. Firstly, this work brings a
contribution in the field of dysarthric speaker adaptation, in which
unsupervised methods dramatically lack while they are essential for a good
recognition accuracy. The second type of novelty is the interpretation of the
method. To the best of our knowledge, this is the first SA study in which the
model provides a measure of similarity between speakers. Note that this
approach may be transferred to other SA contexts in which, for instance, the
algorithm may exploit speakers accents and provide a similarity in terms of
region of origin of the speakers. Last but not least, we moved a step forward
and we introduced new ranking indices that we used to perform dysarthria
detection. This approach is completely new and can achieve high diagnosis
performance, without requiring a specific model training.

Future directions could delve into the dysarthria assessment by individuating
intervals of values in which the DI corresponds to dysarthria severity levels
(e.g., mild, moderate, severe). This may bring to very efficient ASR systems
that simultaneously improve their performance via SA, and compute the DI
warning the subject when the index is close to the right endpoint of its
interval. Such a device could perform an early detection of the disorder or
predict the disease degeneration and allow the patient to act in time in order
to prevent it. Furthermore, the dysarthric index may be used to track the
speaker changes during a treatment and monitoring the benefits of the
rehabilitation.

Finally, our MSDA-WJDOT approach was applied on a limited resource scenario
where the initial unadapted speaker independent model was trained on the
AllSpeak training set. In future work we will experiment with speaker
independent models whose hidden layers are pre-trained on much larger datasets
(e.g., datasets used for ASR).

## References

  * [1] M. R. McNeil, _Clinical management of sensorimotor speech disorders_. Thieme, 2009.
  * [2] F. L. Darley, A. E. Aronson, and J. R. Brown, “Differential diagnostic patterns of dysarthria,” _Journal of speech and hearing research_ , vol. 12, no. 2, pp. 246–269, 1969.
  * [3] E.-B. C. and F. J. A. R., “Automatic speech recognition with deep neural networks for impaired speech,” in _International Conference on Advances in Speech and Language Technologies for Iberian Languages_. Springer, 2016, pp. 97–107.
  * [4] B. Vachhani, C. Bhat, B. Das, and S. K. Kopparapu, “Deep autoencoder based speech features for improved dysarthric speech recognition.” in _Interspeech_ , 2017, pp. 1854–1858.
  * [5] J. N. M. and U. S., “Improving acoustic models in torgo dysarthric speech database,” _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , vol. 26, no. 3, pp. 637–645, 2018.
  * [6] X. Xie, R. Ruzi, X. Liu, and L. Wang, “Variational auto-encoder based variability encoding for dysarthric speech recognition,” _Proc. Interspeech 2021_ , pp. 4808–4812, 2021.
  * [7] M. B. Mustafa, S. S. Salim, N. Mohamed, B. Al-Qatab, and C. E. Siong, “Severity-based adaptation with limited data for asr to aid dysarthric speakers,” _PloS one_ , vol. 9, no. 1, p. e86285, 2014.
  * [8] R. Takashima, T. Takiguchi, and Y. Ariki, “Two-step acoustic model adaptation for dysarthric speech recognition,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6104–6108.
  * [9] J. Shor, D. Emanuel, O. Lang, O. Tuval, M. Brenner, J. Cattiau, F. Vieira, M. McNally, T. Charbonneau, M. Nollstadt _et al._ , “Personalizing asr for dysarthric and accented speech with limited data,” _arXiv preprint arXiv:1907.13511_ , 2019.
  * [10] S. O. C. Morales and S. J. Cox, “Modelling errors in automatic speech recognition for dysarthric speakers,” _EURASIP Journal on Advances in Signal Processing_ , vol. 2009, no. 1, p. 308340, 2009.
  * [11] G. Monge, “Mémoire sur la théorie des déblais et de remblais,” in _Histoire de l’Académie Royale des Sciences de Paris, avec les Mémoires de Mathématique et de Physique pour la même année_ , 1781.
  * [12] L. V. Kantorovich, “On the translocation of masses,” in _Journal of Mathematical Sciences_ , 2006.
  * [13] R. Turrisi, R. Flamary, A. Rakotomamonjy, and M. Pontil, “Multi-source domain adaptation via weighted joint distributions optimal transport,” _arXiv preprint arXiv:2006.12938_ , 2020.
  * [14] A. Genevay, G. Peyré, and M. Cuturi, “Learning generative models with sinkhorn divergences,” _arXiv preprint arXiv:1706.00292_ , 2017.
  * [15] C. Di Nardi, R. Turrisi, A. Inuggi, N. Riva, I. Mauri, and L. Badino, “An automatic speech recognition android app for als patients.”
  * [16] N. Courty, R. Flamary, A. Habrard, and A. Rakotomamonjy, “Joint distribution optimal transportation for domain adaptation,” in _Advances in Neural Information Processing Systems 30_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds. Curran Associates, Inc., 2017, pp. 3730–3739. [Online]. Available: http://papers.nips.cc/paper/6963-joint-distribution-optimal-transportation-for-domain-adaptation.pdf
  * [17] J. Neto, L. Almeida, M. Hochberg, C. Martins, L. Nunes, S. Renals, and T. Robinson, “Speaker-adaptation for hybrid hmm-ann continuous speech recognition system,” 1995.
  * [18] J. R. Williamson, T. F. Quatieri, B. S. Helfer, J. Perricone, S. S. Ghosh, G. Ciccarelli, and D. D. Mehta, “Segment-dependent dynamics in predicting parkinson’s disease,” in _Sixteenth annual conference of the international speech communication association_ , 2015.
  * [19] M. Tu, V. Berisha, and J. Liss, “Interpretable objective assessment of dysarthric speech based on deep neural networks.” in _INTERSPEECH_ , 2017, pp. 1849–1853.
  * [20] A. Hernandez, E. J. Yeo, S. Kim, and M. Chung, “Dysarthria detection and severity assessment using rhythm-based metrics,” in _Proceedings of the Annual Conference of the International Speech Communication Association (INTERSPEECH), Shanghai, China_ , 2020, pp. 25–29.