IFT-UAM/CSIC-23-145, MI-HET-820

# Heavy Neutral Leptons via Axion-Like Particles at Neutrino Facilities

Asli M Abdullahi  Instituto de Física Teórica UAM-CSIC, Universidad Autónoma
de Madrid, Cantoblanco, 28049 Madrid, Spain  Particle Theory Department, Fermi
National Accelerator Laboratory, Batavia, IL, 60510, USA  André de Gouvêa
Northwestern University, Department of Physics & Astronomy, 2145 Sheridan
Road, Evanston, IL 60208, USA  Bhaskar Dutta  Mitchell Institute for
Fundamental Physics and Astronomy, Department of Physics and Astronomy, Texas
A&M University, College Station, TX 77843, USA  Ian M. Shoemaker  Center for
Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, VA 24061,
USA  Zahra Tabrizi  Northwestern University, Department of Physics &
Astronomy, 2145 Sheridan Road, Evanston, IL 60208, USA

###### Abstract

Heavy neutral leptons (HNLs) are often among the hypothetical ingredients
behind nonzero neutrino masses. If sufficiently light, they can be produced
and detected in fixed-target-like experiments. We show that if the HNLs belong
to a richer – but rather generic – dark sector, their production mechanism can
deviate dramatically from expectations associated to the standard-model weak
interactions. In more detail, we postulate that the dark sector contains an
axion-like particle (ALP) that naturally decays into HNLs. Since ALPs mix with
the pseudoscalar hadrons, the HNL flux might be predominantly associated to
the production of neutral mesons (e.g., $\pi^{0}$, $\eta$) as opposed to
charge hadrons (e.g., $\pi^{\pm}$, $K^{\pm}$). In this case, the physics
responsible for HNL production and decay are not directly related and
experiments like DUNE might be sensitive to HNLs that are too weakly coupled
to the standard model to be produced via weak interactions, as is generically
the case of HNLs that play a direct role in the type-I seesaw mechanism.

Introduction. The discovery of neutrino oscillation Fukuda  _et al._ (1998,
1999); Ahmad  _et al._ (2002), and the confirmation of non-zero neutrino
masses, implies the existence of fields beyond those of the standard model
(SM) of particle physics. While little is known about these new degrees of
freedom, a plethora of extensions have been proposed to explain the origin of
the neutrino mass. Many of these extensions invoke the existence of SM-singlet
fermions, commonly referred to as right-handed, or sterile, neutrinos. Popular
examples of such models include the Type-1 seesaw mechanism Minkowski (1977);
Mohapatra and Senjanovic (1980); Gell-Mann  _et al._ (1979); Yanagida (1979);
Lazarides  _et al._ (1981); Mohapatra and Senjanovic (1981); Schechter and
Valle (1980); Cheng and Li (1980); Foot  _et al._ (1989) and other seesaw
variants Schechter and Valle (1980); Konetschny and Kummer (1977); Cheng and
Li (1980); Mohapatra and Senjanovic (1981); Foot  _et al._ (1989); Ma (1998);
Bajc and Senjanovic (2007); Dorsner and Fileviez Perez (2007). Depending on
their properties, these singlet fermions may help shed light on other
outstanding problems in particle physics, including the abundance of dark
matter Boehm and Fayet (2004); Boehm  _et al._ (2004); Pospelov  _et al._
(2008); Pospelov (2009) and the baryon asymmetry of the universe Fukugita and
Yanagida (1986); Davidson  _et al._ (2008).

SM-singlet fermions can interact with SM particles through Yukawa interactions
involving the Higgs-doublet $H$ and the lepton-doublets $L$, in what is known
as the neutrino portal. Such a coupling induces a mixing between the active
neutrinos, $\nu_{e,\mu,\tau}$, and the SM-singlet fermions, $N_{s_{j}}$,
$j=1,\ldots,n_{s}$, generating a Dirac mass for the neutrinos after
electroweak symmetry breaking (EWSB). The singlet fermions may or may not have
a Majorana mass, although such a mass term is always permitted by symmetry. We
are left with mass eigenstates $\nu_{i}$, with masses $m_{i}$, that are linear
superpositions of the active and singlet states, $\nu_{i}=U^{\dagger}_{\alpha
i}\nu_{\alpha}$, where $\alpha$ now runs over
$\alpha=e,\mu,\tau,s_{1},s_{2},\ldots,s_{n_{s}}$. The transformation is
parametrized by the unitary $(3+n_{s})\times(3+n_{s})$ matrix, with components
$U_{\alpha i}$. Henceforth, we will restrict our discussions to $n_{s}=1$ and
refer to the new fermion as $N$. We are interested in mass scales of order 1
MeV to 1 GeV for the new particle, where $N$ is often referred to as a heavy
neutral lepton (HNL). In this region of parameter space, the $U_{\alpha 4}$
elements of the mixing matrix often govern the production of $N$ in fixed-
target experiments, along with their scattering cross sections inside
detectors, and their decay properties.

The existence of SM-singet fermions invites one to consider a richer “dark
sector” with its own particle content and interactions Pospelov (2011); Harnik
_et al._ (2012); Batell  _et al._ (2016); Farzan and Heeck (2016); De Romeri
_et al._ (2017); Magill  _et al._ (2018); Bertuzzo  _et al._ (2019, 2018);
Ballett  _et al._ (2019, 2020a); Coloma (2019); Fischer  _et al._ (2020);
Cline  _et al._ (2020); Berbig  _et al._ (2020). Generically, it is useful to
allow for the possibility that $N$, while a SM-singlet, is charged under a
“dark” gauge group and interacts with other particles that are also SM-
singlets.

We will concentrate on dark sectors that also contain a new pseudoscalar, $a$,
with non-zero mass, sometimes referred to as an axion-like particle (ALP). We
will further assume that $a$ couples to $N$ and that its mass is such that it
decays predominantly to pairs of HNL. This is not an especially extraordinary
hypothesis as it mimics the SM. In the SM, the strong interactions confine and
the lightest propagating hadronic degrees of freedom are pseudoscalars (pions,
kaons, etc). The lightest among these are restricted to decaying weakly into
pairs of leptons (e.g., $\pi^{+}\to\mu^{+}\nu$) or electromagnetically into
photons (e.g., $\pi^{0}\to\gamma\gamma$). If the dark sector is anything like
the SM, it is easy to imagine scenarios where $a\to NN$ is similarly
inevitable. A complete chiral dark sector model with these characteristics was
proposed and explored in Ref. Berryman  _et al._ (2017). There, the ALPs were
referred to as ‘dark pions’, and the HNLs ‘dark neutrinos’. We will refer to
Ref. Berryman  _et al._ (2017) when discussing concrete realizations of the
scenario of interest, but emphasize that our results remain general and not
tied to any specific model.

In the absence of the ALP, in fixed-target-like experiments and for HNL masses
below a GeV, the HNLs are typically produced in the decay of charged mesons
(e.g., $\pi^{+}\rightarrow\mu^{+}N$). In this case, the production rate is
proportional to $|U_{\alpha 4}|^{2}$. Further downstream, the $N$ can be
detected via its decays to SM particles, with the $N$ partial widths also
proportional to $|U_{\alpha 4}|^{2}$. If a kinematically accessible ALP is now
added to the picture, we allow for a new HNL production mechanism in fixed-
target-like experiments. The ALPs mix with neutral pseudoscalar mesons
($\pi^{0},\eta$, etc.) and are produced in tandem with light hadrons at the
target. This production rate is, of course, independent of $U_{\alpha 4}$. The
ALPs then decay to HNL pairs with large branching ratios and HNL decays
downstream can lead to an observable signal in the detector.

While generic ALP-HNL couplings have been studied in the literature, the focus
has typically been on the constraint posed by experiment on the coupling
itself de Giorgi  _et al._ (2023); Alves  _et al._ (2020), or the contribution
to existing processes coming from such a coupling Bonilla  _et al._ (2023). In
this work, we focus on the impact of an ALP-HNL coupling on the discovery
potential of HNLs in fixed-target-like experiments. We explore the
consequences of this new HNL production mechanism at the DUNE near-detector
complex, taken as a particularly relevant experimental setup. Fig. 1 provides
a schematic of the production–decay–detection process. We find that, depending
on the ALP properties, one can explore previously inaccessible regions of the
HNL parameter space, including the region preferred by the type-I seesaw model
for neutrino masses.

Figure 1:  Schematic description of the production of axion-like particles
(ALP, $a$) followed by the production of heavy neutral leptons (HNL, $N$), in
a fixed-target experimental setup. In the DUNE near-detector complex, a beam
of 120 GeV protons impinges on a target, producing charged and neutral mesons.
ALP production is proportional to that of neutral mesons: the crosses
represent $a$–neutral-meson mixing. The ALP travels towards the detector and
decays into a pair of HNLs a distance $h$ from the target. The HNLs that reach
the detector (width $\Delta L$, a distance $L$ from the target) may decay
inside the detector volume into SM particles, yielding an observable signal.
Not to scale.

Heavy ALP production. The ALP $a$ is a pseudo-scalar field and mixes with the
SM pseudoscalar mesons, including the $\pi^{0}$ and the $\eta$. Following
Refs. Krauss and Wise (1986); Alves and Weiner (2018); Altmannshofer  _et al._
(2020) and especially Ref. Kelly  _et al._ (2021), the relevant mixing
contributions of the SM mesons to the ALP can be defined as

$\displaystyle\pi^{0}$ $\displaystyle\rightarrow\pi^{0}+g_{\pi a}a\,,$ (1)
$\displaystyle\eta$ $\displaystyle\rightarrow\eta^{0}+g_{\eta a}a\,,$ (2)

where the meson mixing parameters are Kelly  _et al._ (2021)

$\displaystyle g_{\pi a}$
$\displaystyle=\frac{1}{6}\frac{f_{\pi}}{f_{a}}\frac{m_{a}^{2}}{m_{a}^{2}-m_{\pi^{0}}^{2}}\,,$
(3) $\displaystyle g_{\eta a}$
$\displaystyle=\frac{1}{\sqrt{6}}\frac{f_{\pi}}{f_{a}}\bigg{(}\frac{m_{a}^{2}-\frac{4}{9}m_{\pi^{0}}^{2}}{m_{a}^{2}-m_{\eta}^{2}}\bigg{)}\,.$
(4)

In order to estimate the ALP production rate at a fixed-target facility, we
assume that the probability of producing an $a$ particle with three-momentum
$\vec{p}_{a}$ is the same as that to produce a pseudoscalar $\mathbf{m}$ with
the same three-momentum, up to a scaling factor $g_{\mathbf{m}a}^{2}$. This is
a good approximation as long as the mass of the ALP is similar to the hadronic
pseudoscalar masses. With this in mind, the ALP differential flux is
proportional to that of the pseudoscalar $\mathbf{m}$:

$\displaystyle\frac{d^{2}\phi_{a}}{dE_{a}d\theta_{a}}=g_{\mathbf{m}a}^{2}\frac{E_{a}}{E_{\mathbf{m}}}\frac{d^{2}\phi_{\mathbf{m}}}{dE_{\mathbf{m}}d\theta_{\mathbf{m}}}\,,$
(5)

where $\phi_{a,\mathbf{m}}$ are the ALP and meson fluxes, $E_{a,\mathbf{m}}$
their respective energies, and $\theta_{a}=\theta_{\mathbf{m}}$ is the angle
defined by the outgoing ALP/meson three-momentum and the direction of the
incoming proton beam. For meson fluxes at the DUNE near-detector facility,
assuming the SM, we used the results from Refs. Brdar  _et al._ (2021, 2023),
obtained with GEANT4 Agostinelli  _et al._ (2003)

Integrating Eq. (5) over $\theta_{a}$ and $E_{a}$, Fig. 2 depicts the ALP flux
as a function of the ALP mass $m_{a}$, considering the contributions from
ALP–$\pi^{0}$ (dotted, orange curve) and ALP–$\eta$ (dashed, blue curve)
mixing. We assume a proton-on-target (POT) rate of $1.1\times 10^{21}$ per
year for a period of 10 years Abi  _et al._ (2020), and chose $f_{a}=1$ TeV.
The ALP flux is proportional to $f_{a}^{-2}$ (see Eqs. (3), (4) and (5)). We
restrict the plot to $m_{a}$ values under a few GeV since the mixing formalism
discussed above fails for larger values of $m_{a}$.

Figure 2:  Expected ALP flux at the DUNE near-detector facility as a function
of the ALP mass $m_{a}$. The dotted, orange (dashed, blue) curve corresponds
to the contribution from ALP–$\pi^{0}$ (ALP–$\eta$) mixing while the total
flux is depicted with the solid, red curve. We assume a total exposure of
$1.1\times 10^{22}$ protons on target and $f_{a}=1$ TeV. The ALP flux is
proportional to $f_{a}^{-2}$), see text for details.

HNL production. Once produced, we assume that the ALP decays exclusively into
a pair of HNLs. For the decay rate, we use the model discussed in Ref.
Berryman  _et al._ (2017), where ALP decays into HNLs parallel the decays of
charged-pions into lepton pairs, mediated by $W$-boson exchange. The partial
width is given by,

$\displaystyle\Gamma(a\to 2N)=\frac{1}{4\pi
v_{D}^{4}}f_{a}^{2}m_{a}m_{N}^{2}.$ (6)

Here, $v_{D}$ is the vacuum expectation value of the dark sector Higgs-like
scalar and $m_{N}$ is the HNL mass. Here we assume $v_{D}$ to be much larger
than the weak scale Berryman  _et al._ (2017). Like SM pion decay, $a\to 2N$
is chirality suppressed. For very small $m_{N}$ values, the four-body decay,
$a\to 4N$ is significant or even dominant (one SM parallel to this decay mode
is $\pi^{0}\to Z^{*}Z^{*}\to 4\nu$). See Ref. Berryman  _et al._ (2017) for
details. Here, we restrict our discussion to regions of parameter space where
the two-body final state dominates.

The probability that the ALP decays inside an infinitesimal interval $\Delta
h$ a distance $h$ from the production target is

$\displaystyle\left(\frac{{\rm d}P^{a}_{\rm{decay}}}{{\rm d}h}\right)\Delta
h=\left(\frac{\Gamma}{(\beta\gamma)_{a}}e^{-\Gamma_{a}h/(\beta\gamma)_{a}}\right)\Delta
h\,,$ (7)

where $\Gamma_{a}$ is the ALP decay width and
$(\beta\gamma)_{a}=|\vec{p}_{a}|/m_{a}$ is the boost-velocity factor. We are
interested in decays that occur before the detector, a distance $L$ away:
$0\leqslant h\leqslant L$. In the case of DUNE, $L=574$ m.

In the ALP rest-frame, the two HNLs are produced back-to-back and the decay
is, of course, isotropic. In order to compute the HNL flux, we simulate ALP
decays in the lab frame and assume the detector to be a cylinder with radius
$r$, aligned with the beam direction (“on-axis”). Since the HNLs are produced
at a distance $h$ from the target (see Fig. 1), HNLs reach the detector as
long as their production angles (relative to the beam direction) in the lab
frame are smaller than the opening angle of the detector $\theta_{\rm
det}=\arctan(r/L-h)$.

The total flux of HNLs arriving at the detector can be computed by combining
the differential ALP flux and the decay rate. It is

$\displaystyle\phi_{N}$
$\displaystyle=\int\sin\theta_{a}\frac{d^{2}\phi_{a}}{dE_{a}d\theta_{a}}\frac{{\rm
d}P^{a}_{\rm{decay}}}{{\rm d}h}dE_{a}d\theta_{a}$
$\displaystyle\times\Theta(\theta_{\text{det}}-\theta)P^{N}_{\rm{sur}}\frac{d\Omega}{4\pi}dh$
(8)

where $d\Omega$ ad $\theta$ refer to the direction of the HNL and
$P^{N}_{\rm{sur}}$ is the probability that an HNL reaches the detector before
decaying,

$P^{N}_{\rm{sur}}=e^{-\Gamma_{N}(L-h)/(\beta\gamma)_{N}}.$ (9)

$\Gamma_{N},(\beta\gamma)_{N}$ are, respectively, the decay width and boost-
velocity factor of the HNL.

Fig. 3 depicts the HNL flux at the detector as a function of the HNL mass
$m_{N}$, for different values of $f_{a}$ and $v_{D}$, for
$m_{a}=0.2~{}(0.8)~{}$GeV [solid (dashed) curves] and $|U_{\alpha 4}|^{2}$
equals one for all $\alpha=e,\mu,\tau$. The HNL flux gets enhanced for
$m_{a}=0.5$ due to the resonance-like behavior in Fig 2. The ALP production
rate is proportional to $1/f_{a}^{2}$ while the ALP decay rate into HNLs is
proportional to
$(f_{a}/v_{D}^{2})^{2}\exp{-(f_{a}m_{a}m_{N}/v_{D}^{2})^{2}h/|p_{a}|}$. For
our choices of parameters, ALP always decays between the source and the
detector.

Figure 3: HNL flux at the detector as a function of the HNL mass $m_{N}$,
assuming a total exposure of $1.1\times 10^{22}$ protons on target. Different
colors correspond to different values of $f_{a}$ and $v_{D}$ (labelled), while
the solid (dashed) curves correspond to $m_{a}=0.2~{}(0.8)~{}$GeV.

HNL event rates. HNLs decay via SM interactions through their mixing with the
active neutrinos; we assume the dark sector parameters are such that there are
no other allowed decay modes. The total and partial decay widths are governed
by the HNL mass, $m_{N}$, and the elements of the mixing matrix $U_{\alpha
4}$, $\alpha=e,\mu,\tau$ Berryman  _et al._ (2017); Pal and Wolfenstein
(1982); Gorbunov and Shaposhnikov (2007); Ballett  _et al._ (2020b). HNL
decays inside the detector are potentially observable and constitute our
signal.

For a given exposure time $\mathcal{T}$, the total event rate from HNL decays
into a final state $X$, $N_{\rm signal}(X)$, at the detector is given by a
convolution of the decay probability of the HNL into $X$ and the HNL flux, Eq.
(Heavy Neutral Leptons via Axion-Like Particles at Neutrino Facilities):

$\displaystyle N_{\rm signal}(X)$
$\displaystyle=\mathcal{T}\int\sin\theta_{a}\frac{d^{2}\phi_{a}}{dE_{a}d\theta_{a}}\frac{{\rm
d}P^{a}_{\rm{decay}}}{{\rm d}h}dE_{a}d\theta_{a}$
$\displaystyle\times\Theta(\theta_{\text{det}}-\theta)P^{N}_{\rm{sur}}\frac{d\Omega}{4\pi}dh$
$\displaystyle\times P^{N}_{\rm{decay}}\,\mathbf{Br}(N\to X)$ (10)

where $\mathbf{Br}(N\to X)$ is the branching ratio of the HNL decay into $X$
and the probability the decay happens inside the volume of the detector, with
the length $\Delta L$, is

$P^{N}_{\rm
decay}=e^{-\Gamma_{N}^{\rm{Total}}(L-h)/(\beta\gamma)_{N}}\left[1-e^{-\Gamma_{N}^{\rm{Total}}\Delta
L/(\beta\gamma)_{N}}\right]\,.$ (11)

The HNL lifetimes and the various branching fractions can be found in Appendix
A of ref. Berryman  _et al._ (2017); we depict different branching ratios as a
function of $m_{N}$ in Fig. 5, in the supplemental materials. We consider the
charged-current mediated channels $X=e^{\pm}\pi^{\mp},\mu^{\pm}\pi^{\mp},\nu
e^{\pm}\mu^{\mp}$, the neutral-current mediated channel $X\nu\pi^{0}$, the
“mixed” channels $X=\nu\ell^{+}\ell^{-}$ ($\ell=e,\mu$), and the one-loop
channel $X=\nu\gamma$. We pay special attention to $X=\nu e^{+}e^{-}$ since it
has the largest branching fraction for $m_{N}$ below $m_{\pi}$ and contributes
appreciably for larger masses. At small $m_{N}$, $X=\nu\gamma$ is also of
interest, even though it has a much smaller branching smaller.

DUNE Analysis and Results. In order to estimate the sensitivity of the DUNE
near-detector facility to HNLs produced via ALP decays, we assume a
cylindrical (radius $r=2.6$ m and length $\Delta L=10$ m), on-axis, $50$ ton
liquid argon near-detector, located $L=574$ m from the target. We assume ten
years of data taking with an exposure of $1.1\times 10^{21}$ POT per year Abi
_et al._ (2020). We estimate the sensitivity to the HNL parameter space using
a naive $\chi^{2}$ analysis based on the expected HNL rate $N_{\rm{HNL}}$,
where for $90\%$ CL and 2-degrees of freedom we use $N_{\rm{HNL}}=4.61$. The
model of interest contains eight, assumed-to-be-uncorrelated, parameters:
$(f_{a},v_{D},m_{a},m_{N},|U_{\alpha 4}|^{2})$, $\alpha=e,\mu,\tau$. For each
final state, we perform an independent analysis and “turn on” only one
$|U_{\alpha 4}|^{2}$ at a time (i.e., only one among $|U_{eN}|^{2}$, $|U_{\mu
N}|^{2}$, and $|U_{\tau N}|^{2}$ are nonzero at a given analysis). Unless
otherwise noted, we fix $m_{\pi_{D}}=3m_{N}$, for concreteness.

For simplicity, we do not include backgrounds in our analyses. We expect this
to be a reasonable assumption for some HNL-decay final states, including
$e^{+}e^{-}\nu$, but are quite confident this is not the case for other final
states, where significant beam-induced backgrounds are expected. A dedicated
study is necessary for properly estimating the different backgrounds. We leave
such background studies for future work.

Our results in the $m_{N}\times|U_{\alpha 4}|^{2}$–plane, for fixed
$f_{a}=1~{}$TeV and $v_{D}=10^{2}~{}$TeV are depicted in Fig. 4. The top,
middle, and bottom panels correspond $\nu_{e}$-coupled, $\nu_{\mu}$-coupled,
and $\nu_{\tau}$-coupled HNLs, respectively. For very small $|U_{\alpha 4}|$,
HNL production is severely suppressed since the HNLs are too long-lived to
decay inside the detector. There are regions of enhanced sensitivy at
$m_{N}\sim 40~{}$MeV and $m_{N}\sim 200~{}$MeV. These correspond to the
resonant-production peaks observed in Fig. 2 when $m_{a}\sim 120$ MeV and
$\sim 600$ MeV, respectively. For $m_{N}$ smaller than, roughly, the pion
mass, only the $\nu e^{+}e^{-}$ and $\nu\gamma$ channels are open. These
remain competitive for larger values of $m_{N}$ but are superceded by decay
modes with larger branching ratios. Results for $v_{D}=10^{3}$ TeV are
depicted in Fig. 6, in the supplemental materials.

Figure 4: The sensitivity of the DUNE near detector to the HNLs produced from
the dark pion decay, assuming a total of $1.1\times 10^{22}$ protons-on-target
(POT). The panels correspond to different HNL mixing structures, namely mixing
only with $\nu_{e}$ (top panel), $\nu_{\mu}$ (middle panel) and $\nu_{\tau}$
(bottom panel), respectively.

Regardless of whether the HNL is purely $\nu_{e}$-coupled (only $|U_{e4}|\neq
0$), $\nu_{\mu}$-coupled (only $|U_{\mu 4}|\neq 0$), or $\nu_{\tau}$-coupled
(only $|U_{\tau 4}|\neq 0$), the overall sensitivities are approximately the
same. This is to be contrasted with the standard scenario in which the HNLs
are the product of charged-meson decays. There, due to the comparative sizes
of the respective parent meson fluxes, the sensitivity to $\nu_{\mu}$-coupled
HNLs is superior to that of $\nu_{e}$-coupled HNLs, which is in turn superior
to that of $\nu_{\tau}$-coupled HNLs. This is especially relevant for
$\tau$-coupled HNLs whose production is severely suppressed by the flux of
charm mesons (e.g. $D_{s}$) at the target, or, in the case of lower-energy
facilities, completely absent. The enhanced production of $\tau$-coupled HNLs
makes it possible to better probe the $|U_{\tau 4}|$ coupling, as well as
search for the potential detector signatures of the $\tau$-lepton Conrad  _et
al._ (2010); De Gouvêa  _et al._ (2019); Machado  _et al._ (2020); Dev  _et
al._ (2023); Meighen-Berger  _et al._ (2023); Mammen Abraham  _et al._ (2022).

For the ALP parameters of choice, in this work we are sensitive to
significantly smaller $|U_{\alpha 4}|$ values. This is expected. For charged-
meson-produced HNLs, the signal rate is naively proportional to $|U_{\alpha
4}|^{4}$, while for neutral-meson-produced HNLs, the signal rate is naively
proportional to $|U_{\alpha 4}|^{2}$. For the scenario considered here, we
find that DUNE is sensitive to HNLs whose parameters naturally agree with
expectation from the naive type-I seesaw. These regions of parameter space are
highlighted in grey in Fig. 4 and correspond to $|U_{\alpha 4}|^{2}\leq
0.1~{}{\rm eV}/m_{N}$ Ballett  _et al._ (2020b).

Conclusion. SM-singlet fermions are often among the ingredients behind nonzero
neutrino masses. If this is the case, and if these new particles are light
enough, they can be produced and detected in fixed-target-like experiments,
like the DUNE near-detector facility. Here, we explore the fact that the
production and detection of SM-singlet fermions is model dependent. We show
that if these belong to a richer – but rather generic – dark sector, it is
possible that an experiment like DUNE can explore regions of parameter space
relevant for the neutrino mass puzzle.

In more detail, After electroweak symmetry breaking, SM-singlet fermions, also
referred to as heavy neutral leptons (HNL), naturally mix with SM neutrinos.
If these are light enough, they can be produced via the weak interactions in
the decay of charged mesons (e.g., $\pi^{+}\to\mu^{+}N$) and decay via the
weak interactions into leptons and other light states. On the other hand,
axion-like particles (ALP) will naturally decay, often with large branching
ratio, into HNLs if these belong to the same dark sector. Since ALPs mix with
the SM neutral pseudoscalar hadrons, we explored for the first time the
consequences of the fact that the HNL flux at a fixed-target-like facility
might be predominantly associated to the production of neutral mesons,
including the $\pi^{0}$ and the $\eta$. If this is the case, the physics
responsible for HNL production and decay is not directly related and the
production rate is not suppressed by the active–sterile mixing parameters.
Furthermore, experiments like DUNE might be sensitive to HNLs that are too
weakly coupled, directly, to the SM to be produced in the decay of charged
mesons. This is naively expected, for example, of HNLs that play a direct part
in the type-I seesaw mechanism.

Our analyses are confined to the DUNE experiment but we anticipate that the
HNL–ALP connection will lead to interesting consequences across the spectrum
of neutrino facilities, including ICARUS, SBND, $\mu$BooNE, COHERENT, CCM, and
FASER$\nu$ and provide a new avenue for future research and experimentation in
the quest to uncover the mysteries of HNLs and their role in particle physics.

###### Acknowledgements.

The work of AdG is supported in part by the U.S. Department of Energy grant
DE-SC0010143 and in part by the NSF grant PHY-1630782. The work of BD is
supported by the U.S. Department of Energy Grant DE-SC0010813. The work of ZT
is supported by the Neutrino Theory Network Program Grant No. DE-
AC02-07CHI11359 and the U.S. Department of Energy under the award number DE-
SC0020250. She also appreciate the hospitality of Aspen Physics Center where
this work was partially developed. The work of IMS is supported by the U.S.
Department of Energy grant DE-SC0020262. AA acknowledges partial support from
the European Union’s Horizon 2020 Research and Innovation Programme under the
Marie Sklodowska-Curie grant agreement no. 860881-HIDDeN and the Marie
Sklodowska-Curie Staff Exchange grant agreement no. 101086085–ASYMMETRY; from
Fermi National Accelerator Laboratory, managed and operated by Fermi Research
Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of
Energy; and from the Spanish Research Agency (Agencia Estatal de
Investigación) through the grant IFT Centro de Excelencia Severo Ochoa no.
CEX2020-001007-S and the grant RYC2018-024240-I, funded by
MCIN/AEI/10.13039/501100011033 and by “ESF Investing in your future”.

##  Appendix A Supplemental Material

For the expressions of the partial HNL decay widths we have used the
expressions in appendix A of Berryman  _et al._ (2017), as well as Ref. Pal
and Wolfenstein (1982); Gorbunov and Shaposhnikov (2007); Ballett  _et al._
(2020b). In general with all the decay channels that are kinetically allowed
for the mass range we consider in this work, the total decay width reads:

$\displaystyle\Gamma_{N}^{\rm{Total}}=\,\,$ $\displaystyle
2\sum_{\alpha=e,\mu,\tau}\Big{[}\Gamma(N\to\nu_{\alpha}\bar{\nu}\nu)+\Gamma(N\to\nu_{\alpha}\gamma)$
$\displaystyle+\Gamma(N\to\nu_{\alpha}\pi^{0})\Big{]}+2\sum_{\alpha=e,\mu}\Big{[}\Gamma(N\to\pi^{+}\ell_{\alpha})$
$\displaystyle+\Gamma(N\to\ell^{-}_{\alpha}\ell^{+}\nu)+\Gamma(N\to\ell^{-}_{\alpha}\ell^{+}\nu)$
$\displaystyle+\Gamma(N\to\nu_{\alpha}\ell^{-}\ell^{+})\Big{]}\,,$ (12)

where the overall factor of 2 accounts for the fact that the HNLs are Majorana
fermions. We have shown the branching fractions in Fig. 5, where we have
assumed all the mixing angles $|U_{\alpha 4}|^{2}$ are equal, and can be
cancelled from the expressions of the branching fractions.

Figure 5: Branching fractions for different decay modes of the heavy neutral
lepton (HNL), as a function of the HNL mass $m_{N}$.  Figure 6: Same as Fig.
4, for $v_{D}=10^{3}~{}$TeV.

To demonstrate the dependence of our final results on the value of the dark
vev $v_{D}$, we have made Fig. 6, for which all the parameters are the same as
in Fig. 4, but $v_{D}=10^{3}~{}$TeV.

## References

  * Fukuda  _et al._ (1998) Y. Fukuda _et al._ (Super-Kamiokande), Phys. Rev. Lett. 81, 1562 (1998), arXiv:hep-ex/9807003 .
  * Fukuda  _et al._ (1999) Y. Fukuda _et al._ (Super-Kamiokande), Phys. Rev. Lett. 82, 2644 (1999), arXiv:hep-ex/9812014 .
  * Ahmad  _et al._ (2002) Q. R. Ahmad _et al._ (SNO), Phys. Rev. Lett. 89, 011301 (2002), arXiv:nucl-ex/0204008 .
  * Minkowski (1977) P. Minkowski, Phys. Lett. B 67, 421 (1977).
  * Mohapatra and Senjanovic (1980) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).
  * Gell-Mann  _et al._ (1979) M. Gell-Mann, P. Ramond, and R. Slansky, Conf. Proc. C 790927, 315 (1979), arXiv:1306.4669 [hep-th] .
  * Yanagida (1979) T. Yanagida, Conf. Proc. C 7902131, 95 (1979).
  * Lazarides  _et al._ (1981) G. Lazarides, Q. Shafi, and C. Wetterich, Nucl. Phys. B 181, 287 (1981).
  * Mohapatra and Senjanovic (1981) R. N. Mohapatra and G. Senjanovic, Phys. Rev. D 23, 165 (1981).
  * Schechter and Valle (1980) J. Schechter and J. W. F. Valle, Phys. Rev. D 22, 2227 (1980).
  * Cheng and Li (1980) T. P. Cheng and L.-F. Li, Phys. Rev. D 22, 2860 (1980).
  * Foot  _et al._ (1989) R. Foot, H. Lew, X. G. He, and G. C. Joshi, Z. Phys. C 44, 441 (1989).
  * Konetschny and Kummer (1977) W. Konetschny and W. Kummer, Phys. Lett. B 70, 433 (1977).
  * Ma (1998) E. Ma, Phys. Rev. Lett. 81, 1171 (1998), arXiv:hep-ph/9805219 .
  * Bajc and Senjanovic (2007) B. Bajc and G. Senjanovic, JHEP 08, 014 (2007), arXiv:hep-ph/0612029 .
  * Dorsner and Fileviez Perez (2007) I. Dorsner and P. Fileviez Perez, JHEP 06, 029 (2007), arXiv:hep-ph/0612216 .
  * Boehm and Fayet (2004) C. Boehm and P. Fayet, Nucl. Phys. B 683, 219 (2004), arXiv:hep-ph/0305261 .
  * Boehm  _et al._ (2004) C. Boehm, P. Fayet, and J. Silk, Phys. Rev. D 69, 101302 (2004), arXiv:hep-ph/0311143 .
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Lett. B 662, 53 (2008), arXiv:0711.4866 [hep-ph] .
  * Pospelov (2009) M. Pospelov, Phys. Rev. D 80, 095002 (2009), arXiv:0811.1030 [hep-ph] .
  * Fukugita and Yanagida (1986) M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986).
  * Davidson  _et al._ (2008) S. Davidson, E. Nardi, and Y. Nir, Phys. Rept. 466, 105 (2008), arXiv:0802.2962 [hep-ph] .
  * Pospelov (2011) M. Pospelov, Phys. Rev. D 84, 085008 (2011), arXiv:1103.3261 [hep-ph] .
  * Harnik  _et al._ (2012) R. Harnik, J. Kopp, and P. A. N. Machado, JCAP 07, 026 (2012), arXiv:1202.6073 [hep-ph] .
  * Batell  _et al._ (2016) B. Batell, M. Pospelov, and B. Shuve, JHEP 08, 052 (2016), arXiv:1604.06099 [hep-ph] .
  * Farzan and Heeck (2016) Y. Farzan and J. Heeck, Phys. Rev. D 94, 053010 (2016), arXiv:1607.07616 [hep-ph] .
  * De Romeri  _et al._ (2017) V. De Romeri, E. Fernandez-Martinez, J. Gehrlein, P. A. N. Machado, and V. Niro, JHEP 10, 169 (2017), arXiv:1707.08606 [hep-ph] .
  * Magill  _et al._ (2018) G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Phys. Rev. D 98, 115015 (2018), arXiv:1803.03262 [hep-ph] .
  * Bertuzzo  _et al._ (2019) E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, Phys. Lett. B 791, 210 (2019), arXiv:1808.02500 [hep-ph] .
  * Bertuzzo  _et al._ (2018) E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, Phys. Rev. Lett. 121, 241801 (2018), arXiv:1807.09877 [hep-ph] .
  * Ballett  _et al._ (2019) P. Ballett, M. Hostert, and S. Pascoli, Phys. Rev. D 99, 091701 (2019), arXiv:1903.07590 [hep-ph] .
  * Ballett  _et al._ (2020a) P. Ballett, M. Hostert, and S. Pascoli, Phys. Rev. D 101, 115025 (2020a), arXiv:1903.07589 [hep-ph] .
  * Coloma (2019) P. Coloma, Eur. Phys. J. C 79, 748 (2019), arXiv:1906.02106 [hep-ph] .
  * Fischer  _et al._ (2020) O. Fischer, A. Hernández-Cabezudo, and T. Schwetz, Phys. Rev. D 101, 075045 (2020), arXiv:1909.09561 [hep-ph] .
  * Cline  _et al._ (2020) J. M. Cline, M. Puel, and T. Toma, JHEP 05, 039 (2020), arXiv:2001.11505 [hep-ph] .
  * Berbig  _et al._ (2020) M. Berbig, S. Jana, and A. Trautner, Phys. Rev. D 102, 115008 (2020), arXiv:2004.13039 [hep-ph] .
  * Berryman  _et al._ (2017) J. M. Berryman, A. de Gouvêa, K. J. Kelly, and Y. Zhang, Phys. Rev. D 96, 075010 (2017), arXiv:1706.02722 [hep-ph] .
  * de Giorgi  _et al._ (2023) A. de Giorgi, L. Merlo, and J.-L. Tastet, Fortsch. Phys. 71, 2300027 (2023), arXiv:2212.11290 [hep-ph] .
  * Alves  _et al._ (2020) A. Alves, A. G. Dias, and D. D. Lopes, JHEP 08, 074 (2020), arXiv:1911.12394 [hep-ph] .
  * Bonilla  _et al._ (2023) J. Bonilla, B. Gavela, and J. Machado-Rodríguez, (2023), arXiv:2309.15910 [hep-ph] .
  * Krauss and Wise (1986) L. M. Krauss and M. B. Wise, Phys. Lett. B 176, 483 (1986).
  * Alves and Weiner (2018) D. S. M. Alves and N. Weiner, JHEP 07, 092 (2018), arXiv:1710.03764 [hep-ph] .
  * Altmannshofer  _et al._ (2020) W. Altmannshofer, S. Gori, and D. J. Robinson, Phys. Rev. D 101, 075002 (2020), arXiv:1909.00005 [hep-ph] .
  * Kelly  _et al._ (2021) K. J. Kelly, S. Kumar, and Z. Liu, Phys. Rev. D 103, 095002 (2021), arXiv:2011.05995 [hep-ph] .
  * Brdar  _et al._ (2021) V. Brdar, B. Dutta, W. Jang, D. Kim, I. M. Shoemaker, Z. Tabrizi, A. Thompson, and J. Yu, Phys. Rev. Lett. 126, 201801 (2021), arXiv:2011.07054 [hep-ph] .
  * Brdar  _et al._ (2023) V. Brdar, B. Dutta, W. Jang, D. Kim, I. M. Shoemaker, Z. Tabrizi, A. Thompson, and J. Yu, Phys. Rev. D 107, 055043 (2023), arXiv:2206.06380 [hep-ph] .
  * Agostinelli  _et al._ (2003) S. Agostinelli _et al._ (GEANT4), Nucl. Instrum. Meth. A 506, 250 (2003).
  * Abi  _et al._ (2020) B. Abi _et al._ (DUNE), (2020), arXiv:2002.03005 [hep-ex] .
  * Pal and Wolfenstein (1982) P. B. Pal and L. Wolfenstein, Phys. Rev. D 25, 766 (1982).
  * Gorbunov and Shaposhnikov (2007) D. Gorbunov and M. Shaposhnikov, JHEP 10, 015 (2007), [Erratum: JHEP 11, 101 (2013)], arXiv:0705.1729 [hep-ph] .
  * Ballett  _et al._ (2020b) P. Ballett, T. Boschi, and S. Pascoli, JHEP 03, 111 (2020b), arXiv:1905.00284 [hep-ph] .
  * Conrad  _et al._ (2010) J. Conrad, A. de Gouvea, S. Shalgar, and J. Spitz, Phys. Rev. D 82, 093012 (2010), arXiv:1008.2984 [hep-ph] .
  * De Gouvêa  _et al._ (2019) A. De Gouvêa, K. J. Kelly, G. V. Stenico, and P. Pasquini, Phys. Rev. D 100, 016004 (2019), arXiv:1904.07265 [hep-ph] .
  * Machado  _et al._ (2020) P. Machado, H. Schulz, and J. Turner, Phys. Rev. D 102, 053010 (2020), arXiv:2007.00015 [hep-ph] .
  * Dev  _et al._ (2023) P. S. B. Dev, B. Dutta, T. Han, and D. Kim, (2023), arXiv:2304.02031 [hep-ph] .
  * Meighen-Berger  _et al._ (2023) S. A. Meighen-Berger, J. F. Beacom, N. F. Bell, and M. J. Dolan, (2023), arXiv:2311.01667 [hep-ph] .
  * Mammen Abraham  _et al._ (2022) R. Mammen Abraham _et al._ , J. Phys. G 49, 110501 (2022), arXiv:2203.05591 [hep-ph] .