# Up-scattering production of a sterile fermion at DUNE: complementarity with
spallation source and direct detection experiments

Pablo M. Candela  0009-0009-8416-9295<EMAIL_ADDRESS>Instituto de Física
Corpuscular (CSIC-Universitat de València), Parc Científic UV C/ Catedrático
José Beltrán, 2 E-46980 Paterna (Valencia) - Spain  Valentina De Romeri
0000-0003-3585-7437<EMAIL_ADDRESS>Instituto de Física Corpuscular (CSIC-
Universitat de València), Parc Científic UV C/ Catedrático José Beltrán, 2
E-46980 Paterna (Valencia) - Spain  Pantelis Melas  0009-0000-2901-3444
<EMAIL_ADDRESS>Department of Physics, National and Kapodistrian University of
Athens, Zografou Campus GR-15772 Athens, Greece  Dimitrios K. Papoulias
0000-0003-0453-8492<EMAIL_ADDRESS>Department of Physics, National
and Kapodistrian University of Athens, Zografou Campus GR-15772 Athens, Greece
Niki Saoulidou  0000-0001-6958-4196<EMAIL_ADDRESS>Department of
Physics, National and Kapodistrian University of Athens, Zografou Campus
GR-15772 Athens, Greece

###### Abstract

We investigate the possible production of a MeV-scale sterile fermion through
the up-scattering of neutrinos on nuclei and atomic electrons at different
facilities. We consider a phenomenological model that adds a new fermion to
the particle content of the Standard Model and we allow for all possible
Lorentz-invariant non-derivative interactions (scalar, pseudoscalar, vector,
axial-vector and tensor) of neutrinos with electrons and first-generation
quarks. We first explore the sensitivity of the DUNE experiment to this
scenario, by simulating elastic neutrino-electron scattering events in the
near detector. We consider both options of a standard and a tau-optimized
neutrino beams, and investigate the impact of a mobile detector that can be
moved off-axis with respect to the beam. Next, we infer constraints on the
typical coupling, new fermion and mediator masses from elastic neutrino-
electron scattering events induced by solar neutrinos in two current dark
matter direct detection experiments, XENONnT and LZ. Under the assumption that
the new mediators couple also to first-generation quarks, we further set
constraints on the up-scattering production of the sterile fermion using
coherent elastic neutrino-nucleus scattering data from the COHERENT
experiment. Moreover, we set additional constraints assuming that the sterile
fermion may decay within the detector. We finally compare our results and
discuss how these facilities are sensitive to different regions of the
relevant parameter space due to kinematics arguments and can hence provide
complementary information on the up-scattering production of a sterile
fermion.

neutrinos, COHERENT, LZ, XENONnT, DUNE, light mediators, sterile fermion

##  I Introduction

The experimental evidence for neutrino oscillations [1, 2] provides strong
motivation for the existence of new physics beyond the Standard Model (SM). On
the theoretical side, many SM extensions meant to accommodate neutrino masses
call upon the introduction of new sterile states that mix with active
neutrinos (see e.g. [3]). Such sterile neutrinos, whose masses can in
principle range over many orders of magnitude depending on the actual
mechanism that originates neutrino masses, have received increasing attention
in the recent years. Despite their extended phenomenological implications,
comprising cosmological and astrophysical ones, no positive evidence of heavy
sterile neutrinos has been found to date [4, 5, 6]. Nevertheless, searches for
new sterile fermions (and in particular leptons) in the MeV-GeV range are
theoretically motivated and particularly appealing, as they are in principle
accessible at a range of laboratory experiments, including high-intensity and
high-energy facilities (see [6] for a recent review).

A crucial role in these searches may be played by neutrino oscillation
experiments, performing beam-dump style measurements [7]. The Deep Underground
Neutrino Experiment (DUNE) is a next-generation long-baseline neutrino
oscillation experiment with manifold applications besides its primary goal of
a precise determination of neutrino oscillation patterns [8, 9, 10]. Its Near
Detector (ND) complex, located at the Fermi National Accelerator Laboratory
(FNAL), is characterized by strong particle reconstruction capabilities
typical of Liquid Argon (LAr) detectors, which, together with a high-power,
broadband neutrino beam anticipate a great potential in searches for new
physics [11, 12]. The ND will be located approximately 574 m from the neutrino
source and will consist of three primary detector components with the
capability for two of them to move off the beam axis [11]. This feature is
called the DUNE Precision Reaction-Independent Spectrum Measurement (PRISM)
and is primarily designed to produce data-driven predictions of the oscillated
neutrino event rate spectrum at the far detector [13, 11], and to help
constraining neutrino cross sections and flux-related uncertainties [14]. The
liquid argon time-projection chamber (LArTPC) is expected to have a 67 tonne
fiducial volume and to be mobile up to 30.5 m off-axis [11]. With these
characteristics, the DUNE-ND — specially if equipped with the off-axis option
— will offer excellent opportunities to search for low-energy new physics
[15], for instance in the form of heavy neutral leptons [16, 17, 18, 19] or
light DM [20, 21, 22].

Besides “traditional” neutrino oscillation experiments [23, 24], the recent
observation of coherent elastic neutrino-nucleus scattering (CE$\nu$NS) [25,
26, 27] has proven to be a powerful tool to probe the neutrino sector and
search for new physics [28, 29, 30, 31]. Experiments using neutrinos from
stopped-pion beams, like COHERENT [32], Coherent CAPTAIN-Mills [33] and the
planned European Spallation Source [34] are specially suited to study new
physics in the form of neutrino non-standard interactions (NSI) [35] or
neutrino generalized interactions (NGI) [36, 37, 38, 31]. These experiments
can also probe the existence of new light [39, 40, 41, 42, 31] or heavy states
[43, 44], the inelastic transition into a new, heavy, sterile neutrino through
the presence of a active-sterile transition magnetic moment [45, 42, 46, 31]
or the up-scattering production of a new fermion at the MeV scale through new
scalar or vector mediators [47, 48, 49, 50, 51].

Apace with neutrino experiments, there is another vast array of facilities
that have currently reached impressive sensitivities: dark matter (DM) direct
detection (DD) experiments [52, 53]. These underground detectors looking for
ghostly particles of cosmological or astrophysical origin share many common
features with neutrino experiments. Despite their searches for DM in the form
of weakly interacting massive particles (WIMPs) have given null results till
now [54], the incredibly low thresholds, large volumes and background
reduction of these facilities allow them to observe astrophysical neutrinos.
While constituting an irreducible background for WIMP searches [55, 56, 57,
58, 59, 60], the observation of solar neutrinos at DM DD experiments may
provide new unique opportunities to improve the neutrino physics picture [61,
62, 63, 64, 65, 66, 67, 68, 69, 70, 71]. The latest generation of dual-phase
liquid xenon (LXe) detectors like those employed by the XENON [72, 73, 74] and
LUX-ZEPLIN (LZ) [75, 76, 77, 78] experiments, have reached world-leading
sensitivities not only on WIMPs but also on new-physics scenarios beyond the
SM. Solar (and other sources, e.g., atmospheric and diffuse supernova
background) neutrinos can induce not only CE$\nu$NS but also elastic neutrino-
electron scattering (E$\nu$ES) events in DM DD experiments, both having a
strong degeneracy with the corresponding DM-electron [64, 79, 80] and DM-
nucleus signals [58, 60, 81]. However, due to energy-threshold limitations,
solar neutrinos at DM DD experiments are currently detected mainly through
E$\nu$ES [82, 69, 83, 84, 71].

In this paper we explore the complementarity of different facilities in
probing the up-scattering production of a new sterile fermion $\chi$ with a
mass in the MeV range. We consider a simple, phenomenological extension of the
SM with the addition of $\chi$, a generic sterile fermion, and new neutrino-
quark and neutrino-electron interactions (scalar, pseudoscalar, vector, axial-
vector or tensor). In this regard, we extend the NGI scenario [85, 36, 38, 31,
86] by allowing for the inelastic transition of an active neutrino into a new
fermion. For the sake of completeness, we work in a light mediator regime thus
exploring a wider region of parameter space. We are mainly motivated to
explore the sensitivity to this scenario through E$\nu$ES events at the DUNE-
ND. However, in order to provide a more complete picture, we also analyze
current data from DM DD experiments, namely XENONnT [74] and LZ [77], hence
testing the up-scattering production of $\chi$ in E$\nu$ES events induced by
solar (mainly $pp$ and 7Be) neutrinos. Moreover, and under the assumption that
the new mediator couples not only to electrons but also to quarks, we
investigate the up-scattering production of $\chi$ through CE$\nu$NS by
analyzing the most recent COHERENT CsI dataset [27]. Given the different
neutrino sources, solar neutrinos for XENONnT and LZ, neutrinos from pions
decaying at rest for COHERENT and a proton beam for DUNE, these facilities
will be able to explore different regions in parameter space, being sensitive
to different new fermion masses.

One similar framework is the dipole portal, that consists in the possible
production of a heavy neutral lepton due to the presence of active-sterile
neutrino transition magnetic moments [87, 88, 89, 90]. The implications of
this scenario have been studied for DUNE [91, 92, 93], COHERENT [31, 46, 42]
and DM DD experiments [94, 51, 95]. Here, we go one step beyond this scenario
by considering not electromagnetic-type but generalized interactions that may
lead to the up-scattering production of $\chi$. Searches for new generalized
interactions have also been extensively done in the literature, for example at
DUNE [96, 86], neutrino-electron scattering experiments [97, 98], DM DD
experiments [99, 100] and CE$\nu$NS experiments [38, 36].

The remainder of the work is organized as follows: in Section II we discuss
the theoretical scenario under consideration and the relevant E$\nu$ES and
CE$\nu$NS cross sections for the new fermion production. In Section III we
present the details of the analyses of current DM DD and COHERENT data,
together with the simulation of DUNE sensitivity. In Section IV we show the
results of our analyses and we present our final discussion in Section V.

##  II Neutrino up-scattering into a sterile fermion

We consider the possible production of a sterile fermion (SF) $\chi$ from the
up-scattering of neutrinos on atomic electrons or nuclei:

$\nu_{\ell}e\rightarrow\chi
e,\qquad\nu_{\ell}\mathcal{N}\rightarrow\chi\mathcal{N},$ (1)

where $\nu_{\ell}$ is a neutrino or antineutrino of any flavor
$\ell=e,\mu,\tau$, while the nature of the sterile fermion is that of a
spin-1/2 field. In the case of an incoming antineutrino, the corresponding
produced SF is its antiparticle.

Working in the framework of a simplified phenomenological scenario, we extend
the SM Lagrangian at energies below the electroweak scale with

$\mathcal{L}_{\mathrm{SF}}^{a}\supseteq\dfrac{G_{F}}{\sqrt{2}}\,\varepsilon_{\ell}^{a}\left(\overline{\chi}\,\Gamma^{a}P_{L}\,\nu_{\ell}\right)\left(\overline{f}\,\Gamma_{a}f\right)+\mathrm{H.c.}$
(2)

This Lagrangian encodes new neutrino-$\chi$ interactions via one of the
following generalized interactions: scalar ($S$), pseudoscalar ($P$), vector
($V$), axial-vector ($A$) or tensor ($T$). These are summarized through the
subscript $a\equiv\left\\{S,\,P,\,V,\,A,\,T\right\\}$ as
$\Gamma^{a}\equiv\left\\{I,\,i\gamma^{5},\,\gamma^{\mu},\,\gamma^{\mu}\gamma^{5},\,\sigma^{\mu\nu}\right\\}$,
with $\sigma^{\mu\nu}\equiv\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$. Here,
$G_{F}$ is the Fermi coupling constant, $\varepsilon_{\ell}^{a}$ is a coupling
that quantifies the strength of the interaction of type $a$ for a neutrino
$\nu_{\ell}$, $\chi$ is the SF field, $P_{L}\equiv(1-\gamma^{5})/2$ is the
left-handed projector and $f$ is the fermion field that corresponds to
electrons in E$\nu$ES or up and down quarks in CE$\nu$NS, i.e.,
$f=u,\,d,\,e^{-}$.

Given the Lagrangian (2), one can compute the cross section for each
interaction in terms of the target’s recoil energy. Concerning the latter, DM
DD experiments like LZ and XENONnT involve a typical momentum transfer of
$\mathcal{O}$(100) keV, while for COHERENT it is of the order of
$\mathcal{O}$(10) MeV and finally for the case of DUNE it is in the range of a
few GeV. The magnitude of the three-momentum transfer is defined as
$|\mathbf{q}|\approx\sqrt{2m_{e}T_{e}}$ for scattering on electrons or
$|\mathbf{q}|\approx\sqrt{2m_{\mathcal{N}}T_{\mathcal{N}}}$ for scattering on
nuclei. It is not unreasonable to think about a light mediator with a mass
comparable to the typical momentum transfer of these experiments. Thus, we
extend our discussion by including the possibility of having light mediators
and we translate the effective-interaction couplings through the following
substitution

$G_{F}^{2}|\varepsilon_{\ell}^{a}|^{2}\rightarrow\dfrac{2g_{a}^{4}}{(m_{a}^{2}+|\mathbf{q}|^{2})^{2}},$
(3)

where $m_{a}$ is the mass of the mediator of type
$a\equiv\left\\{S,\,P,\,V,\,A,\,T\right\\}$ and $g_{a}$ is the new coupling
that quantifies the strength of the interaction. In Eq. (3) we have re-scaled
the coupling constant and we have explicitly indicated the propagator
dependence on the mediator mass. In addition, we have assumed that the
coupling $\varepsilon_{\ell}^{a}$ is the same for every neutrino flavor
$\nu_{\ell}$. After this transformation, the resulting expressions for the
E$\nu$ES cross section are

$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}}\right|_{\mathrm{\mathrm{E\nu
ES}}}^{\mathrm{S}}(E_{\nu},T_{e})=$
$\displaystyle\,\dfrac{m_{e}g_{S}^{4}}{4\pi(m_{S}^{2}+2m_{e}T_{e})^{2}}Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)\left(1+\dfrac{T_{e}}{2m_{e}}\right)\left(\dfrac{m_{e}T_{e}}{E_{\nu}^{2}}+\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\right),$
(4)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}}\right|_{\mathrm{\mathrm{E\nu
ES}}}^{\mathrm{P}}(E_{\nu},T_{e})=$
$\displaystyle\,\dfrac{m_{e}g_{P}^{4}}{4\pi(m_{P}^{2}+2m_{e}T_{e})^{2}}Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)\dfrac{T_{e}}{2m_{e}}\left(\dfrac{m_{e}T_{e}}{E_{\nu}^{2}}+\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\right),$
(5)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}}\right|_{\mathrm{\mathrm{E\nu
ES}}}^{\mathrm{V}}(E_{\nu},T_{e})=$
$\displaystyle\,\dfrac{m_{e}g_{V}^{4}}{2\pi(m_{V}^{2}+2m_{e}T_{e})^{2}}Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)$
$\displaystyle\times\left[\left(1-\dfrac{m_{e}T_{e}}{2E_{\nu}^{2}}-\dfrac{T_{e}}{E_{\nu}}+\dfrac{T_{e}^{2}}{2E_{\nu}^{2}}\right)-\dfrac{m_{\chi}^{2}}{4E_{\nu}^{2}}\left(1+\dfrac{2E_{\nu}}{m_{e}}-\dfrac{T_{e}}{m_{e}}\right)\right],$
(6)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}}\right|_{\mathrm{\mathrm{E\nu
ES}}}^{\mathrm{A}}(E_{\nu},T_{e})=$
$\displaystyle\,\dfrac{m_{e}g_{A}^{4}}{2\pi(m_{A}^{2}+2m_{e}T_{e})^{2}}Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)$
$\displaystyle\times\left[\left(1+\dfrac{m_{e}T_{e}}{2E_{\nu}^{2}}-\dfrac{T_{e}}{E_{\nu}}+\dfrac{T_{e}^{2}}{2E_{\nu}^{2}}\right)+\dfrac{m_{\chi}^{2}}{4E_{\nu}^{2}}\left(1-\dfrac{2E_{\nu}}{m_{e}}+\dfrac{T_{e}}{m_{e}}\right)\right],$
(7)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}}\right|_{\mathrm{\mathrm{E\nu
ES}}}^{\mathrm{T}}(E_{\nu},T_{e})=$
$\displaystyle\,\dfrac{4m_{e}g_{T}^{4}}{\pi(m_{T}^{2}+2m_{e}T_{e})}Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)$
$\displaystyle\times\left[\left(1-\dfrac{m_{e}T_{e}}{4E_{\nu}^{2}}-\dfrac{T_{e}}{E_{\nu}}+\dfrac{T_{e}^{2}}{4E_{\nu}^{2}}\right)-\dfrac{m_{\chi}^{2}}{4E_{\nu}^{2}}\left(\dfrac{1}{2}+\dfrac{2E_{\nu}}{m_{e}}-\dfrac{T_{e}}{2m_{e}}\right)\right],$
(8)

where $T_{e}$ is the electron recoil energy, $E_{\nu}$ is the incident
neutrino energy, $m_{e}$ is the electron mass and $m_{\chi}$ is the SF mass.
Notably, the results for the vector and scalar cases are consistent with the
cross sections presented in Ref. [50] (and [98] in the regime of effective
interactions), while in the limit of $m_{\chi}\to 0$ the NGI E$\nu$ES cross
sections presented in Ref. [86] are recovered. Tabulated values of the
effective charges $Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)$,
indicating the amount of atomic electrons that can be ionized for a given
energy deposition $T_{e}$ and different isotopes, are given in Appendix A. Let
us clarify that, based upon the NGI definition given in Lagrangian (2) that
denotes all possible four-fermion non-derivative Lorentz structures, we have
included the tensor case also in the light-mediator regime, in the spirit of
performing the most comprehensive phenomenological analysis.

In addition to the study of E$\nu$ES, we consider the production of the SF via
up-scattering on nuclei, since we are interested in the recent COHERENT data,
involving CE$\nu$NS events. To obtain these cross sections, one can replace in
Eqs. (4–8) the electron mass with the nuclear mass, $m_{e}\rightarrow
m_{\mathcal{N}}$, the effective charge with the weak nuclear form factor,
$Z_{\mathrm{eff}}^{\mathcal{A}}\left(T_{e}\right)\rightarrow
F_{W}^{2}(\left|\mathbf{q}\right|^{2})$, and the quark-level couplings $g_{a}$
with the nuclear-level ones, $C_{a}$, through the typical procedure also
adopted in DM DD searches and described in [101]. Following these steps, we
obtain

$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{N}}}{\mathrm{d}T_{\mathcal{N}}}\right|_{\mathrm{CE\nu
NS}}^{\mathrm{S}}(E_{\nu},T_{\mathcal{N}})=$
$\displaystyle\,\dfrac{m_{\mathcal{N}}C_{S}^{4}}{4\pi(m_{S}^{2}+2m_{\mathcal{N}}T_{\mathcal{N}})^{2}}F_{W}^{2}(\left|\mathbf{q}\right|^{2})\left(1+\dfrac{T_{\mathcal{N}}}{2m_{\mathcal{N}}}\right)\left(\dfrac{m_{\mathcal{N}}T_{\mathcal{N}}}{E_{\nu}^{2}}+\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\right),$
(9)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{N}}}{\mathrm{d}T_{\mathcal{N}}}\right|_{\mathrm{CE\nu
NS}}^{\mathrm{P}}(E_{\nu},T_{\mathcal{N}})=$
$\displaystyle\,\dfrac{m_{\mathcal{N}}C_{P}^{4}}{4\pi(m_{P}^{2}+2m_{\mathcal{N}}T_{\mathcal{N}})^{2}}F_{W}^{2}(\left|\mathbf{q}\right|^{2})\dfrac{T_{\mathcal{N}}}{2m_{\mathcal{N}}}\left(\dfrac{m_{\mathcal{N}}T_{\mathcal{N}}}{E_{\nu}^{2}}+\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\right),$
(10)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{N}}}{\mathrm{d}T_{\mathcal{N}}}\right|_{\mathrm{CE\nu
NS}}^{\mathrm{V}}(E_{\nu},T_{\mathcal{N}})=$
$\displaystyle\,\dfrac{m_{\mathcal{N}}C_{V}^{4}}{2\pi(m_{V}^{2}+2m_{\mathcal{N}}T_{\mathcal{N}})^{2}}F_{W}^{2}(\left|\mathbf{q}\right|^{2})$
$\displaystyle\times\left[\left(1-\dfrac{m_{\mathcal{N}}T_{\mathcal{N}}}{2E_{\nu}^{2}}-\dfrac{T_{\mathcal{N}}}{E_{\nu}}+\dfrac{T_{\mathcal{N}}^{2}}{2E_{\nu}^{2}}\right)-\dfrac{m_{\chi}^{2}}{4E_{\nu}^{2}}\left(1+\dfrac{2E_{\nu}}{m_{\mathcal{N}}}-\dfrac{T_{\mathcal{N}}}{m_{\mathcal{N}}}\right)\right],$
(11)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{N}}}{\mathrm{d}T_{\mathcal{N}}}\right|_{\mathrm{CE\nu
NS}}^{\mathrm{A}}(E_{\nu},T_{\mathcal{N}})=$
$\displaystyle\,\dfrac{m_{\mathcal{N}}C_{A}^{4}}{2\pi(m_{A}^{2}+2m_{\mathcal{N}}T_{\mathcal{N}})^{2}}F_{W}^{2}(\left|\mathbf{q}\right|^{2})$
$\displaystyle\times\left[\left(1+\dfrac{m_{\mathcal{N}}T_{\mathcal{N}}}{2E_{\nu}^{2}}-\dfrac{T_{\mathcal{N}}}{E_{\nu}}+\dfrac{T_{\mathcal{N}}^{2}}{2E_{\nu}^{2}}\right)+\dfrac{m_{\chi}^{2}}{4E_{\nu}^{2}}\left(1-\dfrac{2E_{\nu}}{m_{\mathcal{N}}}+\dfrac{T_{\mathcal{N}}}{m_{\mathcal{N}}}\right)\right],$
(12)
$\displaystyle\left.\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{N}}}{\mathrm{d}T_{\mathcal{N}}}\right|_{\mathrm{CE\nu
NS}}^{\mathrm{T}}(E_{\nu},T_{\mathcal{N}})=$
$\displaystyle\,\dfrac{4m_{\mathcal{N}}C_{T}^{4}}{\pi(m_{T}^{2}+2m_{\mathcal{N}}T_{\mathcal{N}})}F_{W}^{2}(\left|\mathbf{q}\right|^{2})$
$\displaystyle\times\left[\left(1-\dfrac{m_{\mathcal{N}}T_{\mathcal{N}}}{4E_{\nu}^{2}}-\dfrac{T_{\mathcal{N}}}{E_{\nu}}+\dfrac{T_{\mathcal{N}}^{2}}{4E_{\nu}^{2}}\right)-\dfrac{m_{\chi}^{2}}{4E_{\nu}^{2}}\left(\dfrac{1}{2}+\dfrac{2E_{\nu}}{m_{\mathcal{N}}}-\dfrac{T_{\mathcal{N}}}{2m_{\mathcal{N}}}\right)\right].$
(13)

The nuclear form factor accounts for nuclear-physics effects. We adopt the
Klein-Nystrand parametrization111We assume that this form factor is the same
for protons and neutrons. This is a good approximation for the typical
momentum transfer involved at the COHERENT experiment., defined as

$F_{W}(\left|\mathbf{q}\right|^{2})=\dfrac{3\,j_{1}(\left|\mathbf{q}\right|R_{A})}{\left|\mathbf{q}\right|R_{A}(1+a_{k}^{2}\left|\mathbf{q}\right|^{2})}\,,$
(14)

where $j_{1}(x)=\sin(x)/x^{2}-\cos(x)/x$ is the spherical Bessel function of
order one, $R_{A}=1.23\,A^{1/3}~{}\mathrm{fm}$ is the nuclear radius, $A$ is
the mass number and $a_{k}=0.7~{}\mathrm{fm}$ is the Yukawa potential range.
The couplings $C_{a}$ are obtained following the procedure of [101, 102] and
take the form

$\displaystyle C_{S}^{2}$ $\displaystyle\equiv
g_{S}^{2}\left(Z\sum_{q=u,d}\dfrac{m_{p}}{m_{q}}f_{T_{q}}^{(p)}+N\sum_{q=u,d}\dfrac{m_{n}}{m_{q}}f_{T_{q}}^{(n)}\right),$
(15) $\displaystyle C_{P}^{2}$ $\displaystyle\equiv
g_{P}^{2}\left[Z\sum_{q=u,d}\dfrac{m_{p}}{m_{q}}\Delta_{q}^{(p)}+N\sum_{q=u,d}\dfrac{m_{n}}{m_{q}}\Delta_{q}^{(n)}\right]\left(1-\sum_{q^{\prime}=u,d}\dfrac{\overline{m}}{m_{q^{\prime}}}\right),$
(16) $\displaystyle C_{V}^{2}$ $\displaystyle\equiv 3Ag_{V}^{2},$ (17)
$\displaystyle C_{A}^{2}$ $\displaystyle\equiv
g_{A}^{2}\left(Z\sum_{q=u,d}\Delta_{q}^{(p)}+N\sum_{q=u,d}\Delta_{q}^{(n)}\right),$
(18) $\displaystyle C_{T}^{2}$ $\displaystyle\equiv
g_{T}^{2}\left(Z\sum_{q=u,d}\delta_{q}^{(p)}+N\sum_{q=u,d}\delta_{q}^{(n)}\right),$
(19)

where $Z$ is the number of protons, $N=A-Z$ is the number of neutrons, $m_{p}$
and $m_{n}$ are the proton and neutron masses, respectively, and $m_{q}$ are
the quark $q$ masses. $f_{T_{q}}^{(p)}$ and $f_{T_{q}}^{(n)}$ represent the
quark mass contributions to the nucleon (proton and neutron) mass,
$\Delta_{q}^{(p)}$ and $\Delta_{q}^{(n)}$ parametrize the quark spin content
of the nucleon and $\delta_{q}^{(p)}$ and $\delta_{q}^{(n)}$ are tensor
charges that characterize the difference between the spin of quarks and anti-
quarks inside the nucleon. Finally, the constant $\overline{m}$ is defined as
$1/\overline{m}\equiv 1/m_{u}+1/m_{d}+1/m_{s}$. Even though the determination
of some of these parameters is still uncertain, in this work we fix the
following values [102]

$\displaystyle f_{T_{u}}^{(p)}$ $\displaystyle=0.026,$ $\displaystyle
f_{T_{u}}^{(n)}$ $\displaystyle=0.018\,,$ (20) $\displaystyle f_{T_{d}}^{(p)}$
$\displaystyle=0.038,$ $\displaystyle f_{T_{d}}^{(n)}$
$\displaystyle=0.056\,,$ (21) $\displaystyle\Delta_{u}^{(p)}$
$\displaystyle=\Delta_{d}^{(n)}=0.777,$ $\displaystyle\Delta_{d}^{(p)}$
$\displaystyle=\Delta_{u}^{(n)}=-0.438\,,$ (22)
$\displaystyle\delta_{u}^{(p)}$ $\displaystyle=\delta_{d}^{(n)}=0.784,$
$\displaystyle\delta_{d}^{(p)}$ $\displaystyle=\delta_{u}^{(n)}=-0.204\,.$
(23)

Note that in the pseudoscalar and tensor parameters (last two rows) the
isospin symmetry limit has been assumed, i.e.,
$\delta_{u}^{(p)}=\delta_{d}^{(n)}$, $\delta_{d}^{(p)}=\delta_{u}^{(n)}$,
$\Delta_{u}^{(p)}=\Delta_{d}^{(n)}$ and $\Delta_{d}^{(p)}=\Delta_{u}^{(n)}$.
Finally, let us comment that using different values [103, 104] for these
charges would imply, at most, a difference of $27\%$ in the scalar and tensor
couplings, of $10\%$ in the axial-vector coupling $g_{A}$ and of $12\%$ in the
pseudoscalar coupling $g_{P}$.

Before proceeding to our numerical analysis, a final comment should be made,
specially regarding the axial-vector interaction. In this work we first
compute the differential cross sections at the quark level and then follow the
procedure of [101, 102] to obtain the expressions at the nuclear level. A
rigorous procedure would instead require to perform a multipole decomposition
of the hadronic current [105]. Such a treatment calls for extensive nuclear
structure calculations and goes well beyond the scope of our present study,
see e.g. [106]. Compared to conventional studies that rely upon
phenomenological form factors to account for the finite nuclear size, the use
of different nuclear models may lead to differences in the expected event
rates of up to 10% [106, 107]. It is worth noting, however, that nuclear
model-dependent uncertainties as large as $\sim 5\%$ also exist [108]. In our
present work, the latter effects are effectively accounted for through the
introduction of relevant nuisance parameters in the statistical analysis (see
Sec. III).

##  III Data analysis

In this section, we provide a description of the procedure followed to analyze
the most recent COHERENT, LZ and XENONnT datasets as well as the simulation of
E$\nu$ES events at the DUNE-ND.

Given a neutrino flux with an energy endpoint $E_{\nu}^{\mathrm{max}}$, the
maximum and minimum kinetic recoil energies that the atomic electrons in the
detector can acquire, depend on the energy flux endpoint as follows

$T^{\prime\begin{subarray}{c}\scriptscriptstyle\mathrm{min}\\\
\scriptscriptstyle\mathrm{(max)}\end{subarray}}_{e}=\dfrac{2m_{e}{(E_{\nu}^{\mathrm{max}})}^{2}-m_{\chi}^{2}\left(E_{\nu}^{\mathrm{max}}+m_{e}\right){\begin{subarray}{c}-\\\
(+)\end{subarray}}E_{\nu}^{\mathrm{max}}\sqrt{4m_{e}^{2}{(E_{\nu}^{\mathrm{max}})}^{2}-4m_{e}m_{\chi}^{2}\left(E_{\nu}^{\mathrm{max}}+m_{e}\right)+m_{\chi}^{4}}}{2m_{e}\left(2E_{\nu}^{\mathrm{max}}+m_{e}\right)}.$
(24)

Likewise, the minimum neutrino energy needed for the production of the SF is
given, in terms of the true recoil energy of the electron, as

$E_{\nu}^{\mathrm{min}}(T^{\prime}_{e})=\dfrac{1}{2}\left(T^{\prime}_{e}+\sqrt{2m_{e}T^{\prime}_{e}+(T^{\prime}_{e})^{2}}\right)\left(1+\dfrac{m_{\chi}^{2}}{2m_{e}T^{\prime}_{e}}\right).$
(25)

In the case of nuclear recoil these expressions hold under the changes
$m_{e}\rightarrow m_{\mathcal{N}}$ and $T_{e}\rightarrow T_{\mathcal{N}}$.
Besides, in the limit $m_{\chi}\rightarrow 0$, the usual SM expressions are
recovered.

###  III.1 DUNE-ND

For the case of DUNE-ND we simulate the number of E$\nu$ES events in the
$E_{e}\theta_{e}^{2}$ space where $E_{e}=m_{e}+T_{e}$ and $\theta_{e}$ are the
total reconstructed electron energy and scattering angle, respectively. This
way it becomes possible to remove a large part of the background which is
mainly coming from charged-current quasielastic (CCQE) scattering and
misidentified $\pi^{0}$, as explained in Ref. [109, 14]. From the kinematics
of the process one has $1-\cos\theta_{e}=m_{e}(1-y)/E_{e}$ where $y$ denotes
the inelasticity parameter, defined as $y\equiv T_{e}/E_{e}$ for SM E$\nu$ES.
For the case of up-scattering into a SF it instead takes the form
$y\equiv(T_{e}/E_{e})\left[1+m_{\chi}^{2}/(2m_{e}T_{e})\right]$. Then, the
number of events as a function of $E_{e}\theta_{e}^{2}$ can be calculated as

$\left[\dfrac{\mathrm{d}N}{\mathrm{d}E_{e}\theta_{e}^{2}}\right]_{\lambda}=t_{\mathrm{run}}\,N_{e}\,N_{\mathrm{POT}}\sum_{\ell}\int_{E_{\nu}^{\mathrm{min}}}^{E_{\nu}^{\mathrm{max}}}\mathrm{d}E_{\nu}\dfrac{\mathrm{d}\Phi_{\nu_{\ell}}(E_{\nu})}{\mathrm{d}E_{\nu}}\left[\dfrac{\mathrm{d}\sigma_{\nu_{\ell}}}{\mathrm{d}E_{e}\theta_{e}^{2}}\right]_{\lambda}\,,$
(26)

where the differential cross section reads

$\dfrac{\mathrm{d}\sigma_{\nu_{\ell}}}{\mathrm{d}E_{e}\theta_{e}^{2}}=\left.\dfrac{E_{\nu}}{2m_{e}}\dfrac{\mathrm{d}\sigma_{\nu_{\ell}}}{\mathrm{d}T_{e}}\right|_{T_{e}=E_{\nu}\left(1-\frac{E_{e}\theta_{e}^{2}}{2m_{e}}-\Delta\right)}\,.$
(27)

Here, $\Delta=0$ when $\lambda=\mathrm{SM}$ and
$\Delta=m_{\chi}^{2}/\left(2m_{e}E_{\nu}\right)$ when $\lambda=\mathrm{SF}$.
In Eq. (26), $\mathrm{d}\Phi_{\nu_{\ell}}(E_{\nu})/\mathrm{d}E_{\nu}$ denotes
the Long Baseline Neutrino Facility (LBNF) neutrino flux for the different
neutrino flavors
$\nu_{\ell}=\\{\nu_{e},\nu_{\mu},\bar{\nu}_{e},\bar{\nu}_{\mu}\\}$ and the
various on-axis or off-axis locations [110, 111]. We assume the DUNE-ND to be
mobile up to 30 m off-axis [11]. Moreover, when simulating the neutrino
spectra, for each neutrino flavor we consider both the nominal neutrino mode,
namely the _CP-optimized_ beam and a _$\tau$ -optimized_ configuration [112,
113], the latter allowing for a shift of the neutrino spectra to higher
energies.

Regarding the integration limits, the minimum neutrino energy to generate a
detectable recoil is given by Eq. (25). The maximum neutrino energy
$E_{\nu}^{\text{max}}$ is taken by the endpoint of the LBNF beam. Following
[114], the signal spectra are evaluated assuming $50~{}\mathrm{MeV}\leq
E_{e}\leq 20~{}\mathrm{GeV}$ (as discussed in Ref. [17], concerning electron
detection a threshold as low as 30 MeV is possible). Let us stress that a
higher threshold of — for example — 600 MeV would modify our extracted
sensitivities, at most, by a factor of 3 in the tensor case, whereas for the
other interactions this change is milder. It should be kept in mind that the
FNAL beam peaks at 2–3 GeV. Finally, the number of protons on target (POT) per
year assuming a 120 GeV proton beam is taken to be $N_{\mathrm{POT}}=1.1\times
20^{21}$ [115], while $t_{\text{run}}$ is the exposure time and $N_{e}$ is the
number of electron targets contained in the LArTPC of the DUNE-ND, for which
we assume a fiducial volume of 67 tonnes.

We evaluate the number of events by performing a Monte Carlo integration
assuming an angular uncertainty of $\sigma_{\theta}=1^{\mathrm{o}}$, following
the method described in Ref. [86]. Before closing this discussion let us
stress that we have verified that an uncertainty in the reconstruction of
electron energy as large as $\sigma_{E}/E_{e}=10\%/\sqrt{E_{e}/\mathrm{GeV}}$
has negligible impact on the event rates, see e.g. [116, 96, 86].

To extract the sensitivity at the DUNE-ND, we rely on the analysis method of
Ref. [86]. A simultaneous $\chi^{2}$ fit is performed taking into account on-
axis (loc = 0 m) and off-axis (loc = 6, 12, 18, 24 and 30 m) spectra, using
shape and normalization information, through the following $\chi^{2}$ function

$\chi^{2}=2\sum_{k=\nu,\,\bar{\nu}}\,\sum_{j=\mathrm{loc}}\,\sum_{i=1}^{20}\left[N^{\mathrm{th}}_{ijk}-N^{\mathrm{exp}}_{ijk}+N^{\mathrm{exp}}_{ijk}\,\ln{\left(\dfrac{N^{\mathrm{exp}}_{ijk}}{N^{\mathrm{th}}_{ijk}}\right)}\right]+\left(\frac{\alpha_{1}}{\sigma_{\alpha_{1}}}\right)^{2}+\left(\frac{\alpha_{2}}{\sigma_{\alpha_{2}}}\right)^{2}\,.$
(28)

Here, the index $i$ runs over the reconstructed $E_{e}\theta_{e}^{2}$ bins for
which we assume 20 bins in the range $[0,10~{}m_{e}]$ as done in Ref. [86].
The index $j$ runs over the different on-axis and off-axis locations, while
the index $k$ refers to the neutrino (3.5 years) and antineutrino mode (3.5
years). The expected theoretical number of events depends on the nuisance
parameters $\alpha_{1}$ and $\alpha_{2}$ which account for the uncertainty on
the neutrino flux ($\sigma_{\alpha_{1}}=5\%$) and background
($\sigma_{\alpha_{2}}=10\%$) normalizations, and reads
$N_{\mathrm{th}}=(1+\alpha_{1})N_{\mathrm{SM}}+(1+\alpha_{2})N_{\mathrm{bkg}}$.
On the other hand, $N_{\mathrm{exp}}$ accounts for the simulated data that
DUNE-ND would detect including the new-physics SF events. Finally, the
background events are taken to be the sum of CCQE and misidentified $\pi^{0}$
events, i.e., $N_{\text{bkg}}=N_{\pi^{0}}^{\text{missID}}+N_{\text{CCQE}}$
(see Ref. [114]). The latter events are weighted accordingly to match the
exposure time at the various off-axis locations. We also note that although
the shape of the background events is expected to vary for the different on-
axis or off-axis locations, this effect is foreseen to be small.

###  III.2 XENONnT and LUX-ZEPLIN

Solar neutrinos induce low-energy electron recoils via E$\nu$ES events that
are detectable at DM DD experiments. We are interested in extracting limits on
the up-scattering production of the SF using available data from the state-of-
the-art experiments XENONnT and LZ. The signal can be obtained from the
expression

$N^{\mathrm{E\nu
ES}}_{i,k}=\mathcal{E}N_{T}~{}\int_{T_{e}^{i}}^{T_{e}^{i+1}}\mathrm{d}T_{e}~{}\epsilon_{E}(T_{e})~{}\int_{T_{e}^{{}^{\prime}{\mathrm{min}}}}^{T_{e}^{{}^{\prime}\mathrm{max}}}\mathrm{d}T_{e}^{\prime}~{}R(T_{e},T_{e}^{\prime})~{}\sum_{k=pp,\,^{7}\mathrm{Be}}\int_{E_{\nu}^{\mathrm{min}}}^{E_{\nu,k}^{\mathrm{max}}}\mathrm{d}E_{\nu}~{}\sum_{\ell}~{}\dfrac{\mathrm{d}\Phi_{\nu_{\ell}}^{k}}{\mathrm{d}E_{\nu}}~{}\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}^{\prime}}\,,$
(29)

where $T_{e}$ and $T_{e}^{\prime}$ are the reconstructed and true electron
recoil energies, respectively, and $E_{\nu}$ is the neutrino energy. Here,
$\mathcal{E}$ represents the exposure given in units of $\mathrm{tonnes\times
years}$, while $N_{T}$ corresponds to the number of target nuclei (LXe in the
case of XENONnT and LZ) per $\mathrm{tonne}$. For XENONnT (LZ) the exposure
time is $\mathcal{E}=1.16~{}(0.90)~{}\mathrm{tonne\times year}$ (see Refs.
[74, 77]). The number of electron targets depends on the value of the energy-
dependent quantity, $Z_{\mathrm{eff}}^{\mathcal{A}}$, which is taken into
account in the cross section.

The lower integration limit for the neutrino energy is given by Eq. (25). On
the other hand, the upper limit $E_{\nu,k}^{\mathrm{max}}$ in this case
corresponds to the endpoints of the solar neutrino energy spectra and in
particular for $k=pp,\,^{7}\mathrm{Be}(0.861)$, the two relevant flux
components given the current detector thresholds. Likewise, the true recoil
energy limits are given by Eq. (24).

The fluxes $\mathrm{d}\Phi_{\nu_{\ell}}^{k}/\mathrm{d}E_{\nu}$ are given by

$\Phi_{\nu_{e}}^{k}=\Phi_{\nu_{e}}^{k\,\odot}P_{ee},\quad\Phi_{\nu_{\mu}}^{k}=\Phi_{\nu_{e}}^{k\,\odot}\left(1-P_{ee}\right)\cos^{2}\theta_{23},\quad\Phi_{\nu_{\tau}}^{k}=\Phi_{\nu_{e}}^{k\,\odot}\left(1-P_{ee}\right)\sin^{2}\theta_{23},$
(30)

where $\Phi_{\nu_{e}}^{k\,\odot}$ denote the unoscillated fluxes of $\nu_{e}$
produced in the interior of the Sun [117, 118, 119, 120], while the index $k$
corresponds to the various production mechanisms. As anticipated, for XENONnT
and LZ only $k=pp,\,^{7}\text{Be}$ are relevant, whose corresponding flux
normalizations are taken from Ref. [53]. In Eq. (30), $P_{ee}$ represents the
$\nu_{e}$-survival probability at the terrestrial detector accounting for
neutrino oscillations which we take from Ref. [121], while $\theta_{23}$
stands for the atmospheric mixing angle, see e.g. [1].

Finally, the quantities $R(T_{e},T_{e}^{\prime})$ and $\epsilon_{E}(T_{e})$
denote the detector resolution and efficiency, respectively. For the case of
XENONnT we employ the resolution function from Ref. [122] and the efficiency
from Ref. [74]. Similarly, for the case of LZ the resolution is taken from
Refs. [69, 123] and the efficiency from Ref. [78].

Next, let us analyze the most recent data reported by LZ and XENONnT
collaborations. Regarding the LZ experiment, our statistical analysis is based
on the following Poissonian $\chi^{2}$ function [71]

$\chi^{2}=2\left[\sum_{i=1}^{51}N^{\mathrm{th}}_{i}-N^{\mathrm{exp}}_{i}+N^{\mathrm{exp}}_{i}\ln\left(\dfrac{N^{\mathrm{exp}}_{i}}{N^{\mathrm{th}}_{i}}\right)\right]+\sum_{k}\left(\dfrac{\alpha_{k}}{\sigma_{\alpha_{k}}}\right)^{2}+\sum_{k}\left(\dfrac{\beta_{k}}{\sigma_{\beta_{k}}}\right)^{2}\,,$
(31)

where $N^{\mathrm{exp}}_{i}$ is the data reported in Ref. [77] in the $i$-th
recoil-energy bin, $\alpha_{k}$ and $\beta_{k}$ are nuisance parameters of the
backgrounds and neutrino fluxes, respectively, and $\sigma_{\alpha_{k}}$,
$\sigma_{\beta_{k}}$ are their associated uncertainties. Such uncertainties
are given in Ref. [124] for $\sigma_{\alpha_{k}}$ and in Ref. [53] for
$\sigma_{\beta_{k}}$. The predicted number of events, $N^{\mathrm{th}}_{i}$,
for the energy bin $i$ is given by

$N^{\mathrm{th}}_{i}=\sum_{k=pp,\,^{7}\mathrm{Be}}(1+\beta_{k})N^{\mathrm{E\nu
ES}}_{i,k}+\sum_{k}(1+\alpha_{k})N_{i,k}^{\mathrm{bkg}},$ (32)

where $N^{\mathrm{E\nu ES}}_{i,k}$ is given by Eq. (29) for each flux $k$ and
$N_{i,k}^{\mathrm{bkg}}$ is every background described in Table VI of Ref.
[124]. At this point we should stress that although we are varying all the
background components individually according to Eq. (31) as done in Ref. [71],
our present results would differ by less than 5% if we followed the method of
Ref. [123] where all the backgrounds were summed together and a single
nuisance parameter was used. We however, stick to the more complete method of
Ref. [71] despite the computational cost. For the case of XENONnT experiment,
we instead rely on the Gaussian $\chi^{2}$ function

$\chi^{2}=\sum_{i=1}^{30}\left(\dfrac{N^{\mathrm{th}}_{i}-N^{\mathrm{exp}}_{i}}{\sigma_{i}}\right)^{2}+\sum_{k}\left(\dfrac{\beta_{k}}{\sigma_{\beta_{k}}}\right)^{2}\,,$
(33)

where $N^{\mathrm{exp}}_{i}$ and $\sigma_{i}$ denote the experimental data and
the corresponding uncertainties, respectively, taken from Ref. [74]. The
nuisance parameters $\beta_{k}$ and their uncertainties $\sigma_{\beta_{k}}$
are the same as those of Eq. (31). Moreover,

$N^{\mathrm{th}}_{i}=\sum_{k}(1+\beta_{k})N^{\mathrm{E\nu
ES}}_{i,k}+B_{i}^{0}\,.$ (34)

Here, $B_{i}^{0}$ represents the modeled background reported in [74] from
which the SM E$\nu$ES contribution has been removed, and $N^{\mathrm{E\nu
ES}}_{i,k}$ is obtained with Eq. (29). While in principle the more
sophisticated method of Ref. [71] could be adopted like we did for the case of
LZ, here we choose the other one for simplicity. Indeed, the two methods give
anyway identical results, compare, e.g., the results of Refs. [123, 71].

###  III.3 COHERENT

In the case of the COHERENT experiment we analyze the most recent data
released by the COHERENT collaboration from the CsI measurement [27]222The
COHERENT collaboration has also released results using a LAr target [26].
However, we do not analyze this data set here as it currently has less
statistics and hence leads to less stringent bounds than the CsI-2021 data set
(see e.g. Refs. [31, 50]).. We perform a comprehensive analysis exploiting
both the energy and time data from the COHERENT-CsI-2021 measurement following
the fitting procedure described in Refs. [31, 50].

First, we compute the number of SF or SM events expected at the CsI detector,
of mass $m_{\mathrm{det}}=14.6~{}\mathrm{kg}$, and located at a distance
$L=19.3~{}\mathrm{m}$ from the neutrino source, i.e., the Spallation Neutron
Source (SNS). The expected number of events on a nuclear target,
$\mathcal{N}$, per neutrino flavor, $\nu_{\ell}$, in each nuclear recoil
energy bin, $i$, is given by

$\displaystyle N_{i,\nu_{\ell}}^{\mathrm{CE\nu NS}}(\mathcal{N})=$
$\displaystyle\,N_{\mathrm{target}}\int_{T_{\mathcal{N}}^{i}}^{T_{\mathcal{N}}^{i+1}}\mathrm{d}T_{\mathcal{N}}\,\epsilon_{E}(T_{\mathcal{N}})\int_{T_{\mathcal{N}}^{{}^{\prime}{\mathrm{min}}}}^{T_{\mathcal{N}}^{{}^{\prime}\mathrm{max}}}\mathrm{d}T_{\mathcal{N}}^{\prime}\,R(T_{\mathcal{N}},T_{\mathcal{N}}^{\prime})$
$\displaystyle\times\int_{E_{\nu}^{\mathrm{min}}(T_{\mathcal{N}}^{\prime})}^{E_{\nu}^{\mathrm{max}}}\mathrm{d}E_{\nu}\,\dfrac{\mathrm{d}\Phi_{\nu_{\ell}}}{\mathrm{d}E_{\nu}}(E_{\nu})\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{N}}}{\mathrm{d}T_{\mathcal{N}}^{\prime}}(E_{\nu},T_{\mathcal{N}}^{\prime})\,,$
(35)

where
$N_{\mathrm{target}}=N_{\mathrm{A}}m_{\mathrm{det}}/M_{\mathrm{\mathrm{target}}}$
is the number of target nuclei in the detector. The integration limits for
$T_{\mathcal{N}}^{\prime}$ and $E_{\nu}$ are given by Eqs. (24) and (25).
Moreover, $\epsilon_{E}(T_{\mathcal{N}})$ is the energy-dependent detector
efficiency, $R(T_{\mathcal{N}},T_{\mathcal{N}}^{\prime})$ is the energy
resolution function and $\mathrm{d}\Phi_{\nu_{\ell}}/\mathrm{d}E_{\nu}$ denote
the prompt ($\nu_{\mu}$) and delayed ($\nu_{e}$ and $\overline{\nu}_{\mu}$)
neutrino fluxes produced at the SNS, for which we employ the well-known Michel
spectra. The nuclear-recoil spectrum given in Eq. (III.3) is then converted to
an electron-equivalent spectrum using the scintillation response function
provided in Ref. [27] via the quenching factor (QF) measurements reported in
[125]. In order to be able to directly compare our results with the
experimental data, in a final step we need to eventually convert the electron-
recoil spectrum into a photoelectron (PE) spectrum using the light yield
information, as explained in Ref. [27].

Since the CsI detector cannot distinguish between a nuclear or an electron
recoil, we also include in the analysis the corresponding expected number of
E$\nu$ES events on an atomic nucleus $\mathcal{A}$, estimated as

$\displaystyle N_{i,\nu_{\ell}}^{\mathrm{E\nu ES}}(\mathcal{A})=$
$\displaystyle\,N_{\mathrm{target}}\int_{T_{e}^{i}}^{T_{e}^{i+1}}\mathrm{d}T_{e}\,\epsilon_{E}(T_{e})\int_{T_{e}^{{}^{\prime}{\mathrm{min}}}}^{T_{e}^{{}^{\prime}\mathrm{max}}}\mathrm{d}T_{e}^{\prime}\,R(T_{e},T_{e}^{\prime})$
$\displaystyle\times\int_{E_{\nu}^{\mathrm{min}}(T_{e}^{\prime})}^{E_{\nu}^{\mathrm{max}}}\mathrm{d}E_{\nu}\,\dfrac{\mathrm{d}\Phi_{\nu_{\ell}}}{\mathrm{d}E_{\nu}}(E_{\nu})\dfrac{\mathrm{d}\sigma_{\nu_{\ell}\mathcal{A}}}{\mathrm{d}T_{e}^{\prime}}(E_{\nu},T_{e}^{\prime})\,.$
(36)

Similar to the CE$\nu$NS case, the computed E$\nu$ES spectrum is eventually
converted into a PE-spectrum to be directly confronted to the experimental
data.

In addition, we include timing information in the analysis by distributing the
predicted number of events $N_{i,\nu_{\ell}}^{\mathrm{CE\nu NS}}(\mathcal{N})$
in each time bin $j$. We adopt the time distributions
$\mathcal{P}^{\nu_{\ell}}_{T}(t_{\mathrm{rec}})$ of
$\nu_{\ell}=\nu_{\mu},\,\overline{\nu}_{\mu},\,\nu_{e}$ given in [126, 27],
normalized to 6 $\mathrm{\mu s}$. Thus, the expected number of events per each
observed nuclear recoil energy bin $i$ and time bin $j$ is given by

$N_{ij}^{\mathrm{CE\nu
NS}}(\mathcal{N})=\sum_{\nu_{\ell}=\nu_{e},\nu_{\mu},\bar{\nu}_{\mu}}\int_{t_{\mathrm{rec}}^{j}}^{t_{\mathrm{rec}}^{j+1}}\mathrm{d}t_{\mathrm{rec}}\,\mathcal{P}^{\nu_{\ell}}_{T}(t_{\mathrm{rec}},\alpha_{6})\,\epsilon_{T}(t_{\mathrm{rec}})N_{i,\nu_{\ell}}^{\mathrm{CE\nu
NS}}(\mathcal{N}),$ (37)

where $t_{\mathrm{rec}}$ is the reconstructed time and
$\epsilon_{T}(t_{\mathrm{rec})}$ is the time-dependent efficiency. We have
included an additional nuisance parameter on the beam timing, $\alpha_{6}$, to
account for the uncertainty on the beam arrival time (further details are
provided at the end of this discussion). For the E$\nu$ES rate, the timing
information is also included using Eq. (37) but by exchanging the predicted
number of CE$\nu$NS events with the predicted number of E$\nu$ES events given
in Eq. (III.3). Finally, the SF and SM contributions must be summed, and a
weighted sum for the Cs and I nuclei must be performed to obtain the total
number of expected events ($N_{ij}^{\mathrm{CE\nu NS}}$ and
$N_{ij}^{\mathrm{E\nu ES}}$ for CE$\nu$NS and E$\nu$ES, respectively).

For the statistical analysis of the COHERENT-CsI-2021 data we employ the
following Poissonian least-squares function [31]

$\chi^{2}_{\mathrm{CsI}}\Big{|}_{\mathrm{CE\nu NS}+\mathrm{E\nu
ES}}=2\sum_{i=1}^{9}\sum_{j=1}^{11}\left[N^{\mathrm{th}}_{ij}-N_{ij}^{\mathrm{exp}}+N_{ij}^{\mathrm{exp}}\ln\left(\frac{N_{ij}^{\mathrm{exp}}}{N^{\mathrm{th}}_{ij}}\right)\right]\\\
+\sum_{k=0}^{5}\left(\dfrac{\alpha_{k}}{\sigma_{k}}\right)^{2},$ (38)

where $N_{ij}^{\mathrm{exp}}$ is the binned coincidence data given in [27] and
$N^{\mathrm{th}}_{ij}$ is the predicted number of events, given by

$\displaystyle N^{\mathrm{th}}_{ij}=$
$\displaystyle(1+\alpha_{0}+\alpha_{5})N_{ij}^{\mathrm{CE\nu
NS}}(\alpha_{4},\alpha_{6},\alpha_{7})+(1+\alpha_{0})N_{ij}^{\mathrm{E\nu
ES}}(\alpha_{6},\alpha_{7})$
$\displaystyle+(1+\alpha_{1})N_{ij}^{\mathrm{BRN}}(\alpha_{6})+(1+\alpha_{2})N_{ij}^{\mathrm{NIN}}(\alpha_{6})+(1+\alpha_{3})N_{ij}^{\mathrm{SSB}}\,.$
(39)

Here, the number of events for each background are $N_{ij}^{\mathrm{BRN}}$,
$N_{ij}^{\mathrm{NIN}}$ and $N_{ij}^{\mathrm{SSB}}$, which stand for beam
related neutrons (BRN), neutrino induced neutrons (NIN) and steady state
backgrounds (SSB), respectively, all taken from [27]. Equations (38) and (39)
involve several nuisance parameters ($\alpha_{i}$) together with their
associated uncertainties ($\sigma_{i}$). In particular, $\sigma_{0}=11\%$
encodes efficiency and flux uncertainties; $\sigma_{1}=25\%$,
$\sigma_{2}=35\%$ and $\sigma_{3}=2.1\%$ are related to BRN, NIN and SSB
backgrounds, respectively. On the other hand, $\sigma_{5}=3.8\%$ is associated
to the QF. Notice that the number of events $N_{ij}^{\mathrm{CE\nu NS}}$ and
$N_{ij}^{\mathrm{ES}}$ further include: $\alpha_{4}$, that enters the nuclear
form factor and thus affects only the CE$\nu$NS number of events 333This is
done by introducing the nuisance parameter $\alpha_{4}$ in the nuclear radius
entering Eq. (14) such that $R_{A}=1.23\,A^{1/3}(1+\alpha_{4})$., with
$\sigma_{4}=5\%$; $\alpha_{6}$, that accommodates the uncertainty in beam
timing with no prior assigned; and $\alpha_{7}$, that allows for deviations of
the uncertainty in the CE$\nu$NS efficiency. We refer the reader to Ref. [31]
for further details regarding the statistical analysis.

##  IV Results and discussion

In this section, we summarize the results of our analyses, starting with the
prospects for DUNE-ND and then comparing the constraints inferred using
current data from the COHERENT, XENONnT and LZ experiments.

###  IV.1 DUNE-ND sensitivities

In Fig. 1 we present the 90% C.L. sensitivities on the up-scattering
production of a SF at the DUNE-ND for various on-axis and off-axis locations.
We assume 7 years of total running time (3.5 years in neutrino mode and 3.5
years in antineutrino mode). The calculation takes into account the actual
exposure time in each location i.e., 50% of running time is assumed to be on-
axis and 10% to each of the off-axis locations (6, 12, 18, 24 and 30) m. For
each location the corresponding sensitivity is obtained (for a combined
analysis, see below). To avoid overcrowding the plots, here we depict only the
cases correspondent to on-axis, 6 m and 30 m. The results are presented for
the different interaction channels
$a\equiv\left\\{S,\,P,\,V,\,A,\,T\right\\}$. For illustration purposes three
benchmark choices $m_{\chi}=\\{0.1,\,1,\,10\\}\times\,m_{a}$ are depicted. The
left panels are plotted by projecting the contours onto the relevant coupling
[see Eqs. (4) through (8)] and mediator mass $m_{a}$, while the right panels
are plotted with respect to the SF mass $m_{\chi}$. For small masses of the
mediator, $m_{a}$, or of the SF, $m_{\chi}$, the constraints are reaching a
saturation point, regardless on whether the contours are projected with
respect to $m_{a}$ or $m_{\chi}$. On the other hand, for relatively large
$m_{a}$ or $m_{\chi}$, the contours behave differently. Specifically, the
sensitivity loss occurs at different mediator masses $m_{a}$, depending on the
relation between $m_{a}$ and $m_{\chi}$ (see left panels). However, when the
results are projected with respect to $m_{\chi}$, the sensitivity loss occurs
always at the same SF mass, due to the kinematics of the process.

Figure 1: 90% C.L. projected sensitivity of DUNE-ND on the SF parameter space,
after 7 years (3.5 in neutrino mode and 3.5 in antineutrino mode) of data
taking time. Individual sensitivities corresponding to some on-axis and off-
axis locations are shown. The results are presented for the various possible
interactions: scalar/pseudoscalar, vector/axial-vector and tensor. Dotted,
dashed and plain curves account for the three benchmark scenarios:
$m_{\chi}=\\{0.1,\,1,\,10\\}\times\,m_{a}$. Left (Right) panels depict the
sensitivity regions projected on the mediator (SF) mass.

Figure 2: Comparison of the expected DUNE-ND sensitivities at 90% C.L. after 7
years of data-taking time. The results are combined over all possible
locations, accounting for the different exposure times at each location. The
exclusion regions are presented in terms of the mediator mass $m_{a}$ (left)
and the SF mass $m_{\chi}$ (right) for the three benchmark scenarios (see the
main text for details).

As can be seen in Fig. 1, the on-axis location yields the most stringent
results in all cases, while the remaining off-axis locations yield more or
less similar results which become successively less stringent as we move to
further off-axis locations. This was expected, since more neutrinos per unit
time reach the DUNE-ND when the latter is placed on-axis than off-axis. The
advantage of going off-axis is to reduce backgrounds, although at the price of
also reducing the signal. For completeness, although not shown here, we have
also examined the impact of the different neutrino flux shape distributions
expected at the various locations of DUNE-ND. We did so by normalizing them
appropriately by hand such that they all have the same flux normalizations. In
this case we concluded that for all DUNE-ND locations the extracted
constraints become identical, hence rendering the neutrino flux shape not
relevant at the scope of this analysis.

We are now interested in estimating the combined DUNE-ND sensitivities after 7
years of running time. Our result is depicted in Fig. 2 where we also compare
the sensitivity reach of the different interaction channels. As previously, we
consider the three benchmark choices
$m_{\chi}=\\{0.1,\,1,\,10\\}\times\,m_{a}$, while left and right panels
illustrate the contours with respect to $m_{a}$ and $m_{\chi}$, respectively,
with the same behavior as discussed in Fig. 1. From the plots it can be
deduced that the scalar/pseudoscalar (tensor) interactions will be constrained
with the least (most) sensitivity. Indeed, at the leading order, the
scalar/pseudoscalar cross section is suppressed by a factor
$T_{e}\left(m_{\chi}^{2}+2m_{e}T_{e}\right)/\left(4E_{\nu}^{2}m_{e}\right)$
compared to the vector/axial-vector and tensor cross sections. On the other
hand, the tensor cross section is by a factor 8 larger compared to the
vector/axial-vector one. Moreover, as can be seen, the results for
scalar/pseudoscalar or vector/axial-vector are presented with the same contour
curves. This is due to the following reasons. The cross sections for scalar
and pseudoscalar cases have equivalent forms with the only difference being
that the scalar is proportional to $1+T_{e}/(2m_{e})$ while the pseudoscalar
one is proportional to $T_{e}/(2m_{e})$. However, for the GeV recoil energies
involved at DUNE-ND it holds that $T_{e}/(2m_{e})\gg 1$, rendering the two
cross sections practically equivalent. Note that this is not the case for the
analysis of DM DD experiments where keV recoil energies are involved. The
expressions of vector and axial-vector cross sections differ by relative signs
in terms that contribute negligibly to the event rates. In particular the term
$m_{e}T_{e}/(2E_{\nu}^{2})$ appearing in the first parenthesis of Eqs. (6) and
(7) is $\ll 1$, while in the second parenthesis both $E_{\nu}/m_{e}$ and
$T_{e}/m_{e}$ are $\gg 1$. Finally, let us comment on the marginal improvement
of the combined analysis presented in Fig. 2, compared to Fig. 1. As already
explained, the best sensitivity is obtained when the DUNE-ND is on-axis and
the neutrino flux is maximal. Indeed, the combined analysis is mainly driven
by the on-axis events with a small improvement coming from the contributions
corresponding to off-axis events.

###  IV.2 Comparing with COHERENT, XENONnT and LZ

We now proceed to explore the complementarity of a future DUNE-ND measurement
with existing data from CE$\nu$NS and DM DD experiments. In particular, we
analyze the recent CE$\nu$NS data from the COHERENT CsI-2021 measurement as
well as the E$\nu$ES data from the DM DD experiments XENONnT and LZ. To this
purpose we reproduce our DUNE-ND sensitivities discussed previously in Fig. 2
and we superimpose the corresponding exclusion regions in the light of
existing COHERENT, XENONnT and LZ data.

The results are summarized in Figs. 3 to 7, where a comparison of the
respective 90% C.L. sensitivities for scalar, pseudoscalar, vector, axial-
vector and tensor interactions, is presented. As before, the analysis has been
performed for the three benchmark scenarios
$m_{\chi}=\\{0.1,\,1,\,10\\}\times\,m_{a}$. The complementarity of our
extracted sensitivities is evident. As can be seen, the different experiments
dominate the sensitivities in different regions of mediator or SF mass. For
instance, the XENONnT- and LZ-driven constraints are the most stringent for
very low mediator or SF mass. This is due to the very low detection thresholds
achievable in LXe DM DD detectors, combined with the inverse recoil-energy
dependence in the cross sections. Then, for $m_{\chi}\geq 300$ keV the recent
COHERENT CsI-2021 data are improving the constraints placed by DM DD
experiments. At this point, however, it should be noted that the full
parameter space probed by DM DD experiments as well as a large part of it
probed with CE$\nu$NS data from COHERENT fall in the region $\lesssim$ few
MeV, and may be in conflict with cosmological and astrophysical observations.
The comprehensive estimation of these bounds will eventually depend on the
specific model and its thermal history in the early Universe, thus lying
beyond the scope of this work. However, we show in Figs. 3, 4, 5, and 6 a
vertical dot-dashed gray line corresponding to a bound imposed with big bang
nucleosynthesis (BBN) data (see, e.g., [127]) at the scope of general warning.
On the other hand, for sufficiently large mediator or SF mass, i.e., $\sim 2$
MeV, the current sensitivities will be improved drastically in view of future
measurements at DUNE-ND. This is especially true for the pseudoscalar and
axial-vector interactions, which in the case of CE$\nu$NS lead to a spin-
dependent and hence suppressed signal, and therefore translate into less
stringent bounds by COHERENT. By exploiting the GeV neutrino beam offered at
FNAL it becomes feasible to extend the region of the parameter space that is
currently excluded by the COHERENT measurement. Interestingly, in all cases
except scalar interactions, DUNE-ND will start dominating the constraints
being sensitive to mediator and SF masses in a region not in conflict with
cosmological observations.

Figure 3: 90% C.L. exclusion regions for the case of scalar-mediated SF
production. A comparison is given between the constraints coming from the
analysis of existing data by XENONnT, LZ and COHERENT with the projected
sensitivity from future measurements at the DUNE-ND. Left (Right) panels show
the exclusion regions projecting on the mediator (SF) mass, for the three
benchmark scenarios $m_{\chi}=\\{0.1,\,1,\,10\\}\times\,m_{a}$ considered.

Figure 4: Same as Fig. 3, but for pseudoscalar interactions.

Figure 5: Same as Fig. 3, but for vector interactions.

Figure 6: Same as Fig. 3, but for axial-vector interactions.

Figure 7: Same as Fig. 3, but for tensor interactions.

###  IV.3 Impact of the $\tau$-optimized beam at DUNE-ND

DUNE will benefit from operating with different horn configurations. So far,
we have only assumed the default configuration, which is widely considered in
the literature, that is, the CP-optimized configuration. At this point, we
find it useful to explore the impact of the so-called $\tau$-optimized
configuration on the SF sensitivities. The latter results into a neutrino
energy distribution peaking at higher neutrino energies (see e.g. Refs. [112,
113]). This in turn will translate into an increase of SF events with larger
masses compared to the CP-optimized flux.

In Fig. 8 we compare the SF sensitivities obtained utilizing the CP- vs.
$\tau$-optimized fluxes. Here, we only consider the on-axis event rates and 1
year of exposure for both cases. This explains the slight differences of the
sensitivities corresponding to the CP-optimized flux shown here vs. Figs. 1
and 2. As can be seen from the plots, using the $\tau$-optimized flux leads to
slight improvements on the SF sensitivity in the region of large mediator
$m_{a}$ and $m_{\chi}$ masses, for the various possible interaction channels.
Specially for the scalar/pseudoscalar interaction channel one sees that in the
low mass region the CP-optimized flux is more constraining444This is also the
case for the vector/axial-vector and tensor cases as well, but barely
visible., and depending on the $m_{\chi}/m_{a}$ ratio there is eventually a
turning point where the constraint driven by the $\tau$-optimized flux starts
becoming dominant. Notably, this feature also occurred in heavy sterile
neutrino searches in Ref. [93].

Figure 8: Comparison of the expected 90% C.L. sensitivities using the CP-
optimized versus the $\tau$-optimized flux at DUNE-ND. The results depict the
projected sensitivity assuming 1 yr of exposure for the on-axis location only.

###  IV.4 Sterile fermion decays

We finally turn our attention to potential signals that could trigger the
detector due to events arising from SF decays. Indeed, SFs, after being
produced via the up-scattering on nuclei and electrons inside the detector,
could eventually decay within the detector to SM particles, thus leading to
additional detectable signatures. For the sake of simplicity and illustration,
we focus only on the vector interaction. Following Ref. [50], in our framework
there are only two decay modes of $\chi$ at tree-level. Specifically, for
$m_{\chi}>m_{V}$ one has the two-body decay $\chi\to V\nu_{\alpha}$, while for
$m_{\chi}<m_{V}$ the three-body decay $\chi\to\nu_{\alpha}f\bar{f}$ is
possible via an off-shell mediator. The latter decay widths have been
previously computed to be [50] 555Note that in this work we start from an
effective Lagrangian (see Eq. (2)). These expressions are obtained in the
light-mediator case discussed in Ref. [50] and are equivalent under the
assumption $g_{V}\equiv g_{\chi_{L}}=g_{f}$.

$\displaystyle\Gamma(\chi\to V\nu_{\ell})$
$\displaystyle=\dfrac{g_{V}^{2}}{32\pi}m_{\chi}\left(1-\dfrac{m_{V}^{2}}{m_{\chi}^{2}}\right)\left(1+\dfrac{m_{\chi}^{2}}{m_{V}^{2}}-2\dfrac{m_{V}^{2}}{m_{\chi}^{2}}\right),$
$\displaystyle\Gamma(\chi\to\nu_{\ell}\,\ell^{-}\ell^{+})$
$\displaystyle=\dfrac{g_{V}^{4}}{768\pi^{3}}\,m_{\chi}^{5}\,I(m_{\chi},m_{V},\mu_{2})\,,$
(40)

where $I(m_{\chi},m_{V},\mu_{2})$ is expressed as

$\displaystyle
I(m_{\chi},m_{V},\mu_{2})\equiv\int_{0}^{1-4\mu_{2}}\mathrm{d}x_{1}\,$
$\displaystyle\dfrac{\lambda^{1/2}(1-x_{1},\mu_{2},\mu_{2})}{\left[m_{\chi}^{2}(1-x_{1})-m_{V}^{2}\right]^{2}}\dfrac{x_{1}^{2}}{(1-x_{1})^{3}}\bigg{\\{}6(1+2\mu_{2})\bigg{.}$
$\displaystyle\bigg{.}-x_{1}\left[15+\lambda(1-x_{1},\mu_{2},\mu_{2})+24\mu_{2}+3x_{1}^{2}-12x_{1}(1+\mu_{2})\right]\bigg{\\}}\,,$
(41)

with $\mu_{2}\equiv m_{\ell}^{2}/m_{\chi}^{2}$ and
$\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$ being the Källén function.
Equation (40) is valid in both the heavy- and light-mediator regimes. Notice
that the latter decay mode is kinematically open for $m_{\chi}>2m_{\ell}$
666As explained in Ref. [50] radiative decays such as $\chi\to\nu\gamma$,
$\chi\to\nu\gamma\gamma\gamma$ or $\chi\to\nu\nu\nu$ are possible. However,
here we are interested in the analysis of nuclear and electron recoil
signals.. Finally, if also $m_{V}>2m_{\ell}$, the dilepton decay channel is
open and its decay width reads

$\Gamma(V\to\ell^{-}\ell^{+})=\frac{g_{V}^{2}}{12\pi}m_{V}\left(1+2\frac{m_{\ell}^{2}}{m_{V}^{2}}\right)\sqrt{1-4\frac{m_{\ell}^{2}}{m_{V}^{2}}}\,.$
(42)

Notice that decay channels involving mesons could also be allowed, depending
on the SF mass. While not kinematically accessible at LZ, XENONnT or COHERENT,
the decay mode with a $\pi^{0}$ in the final state would be — in principle —
allowed at DUNE-ND, but only in a marginal corner of parameter space. For this
reason, we did not consider it and we instead focus our discussion on the
decays given in Eq. (40).

We estimate the number of decay events inside the detector, as

$N_{\mathrm{decay}}=N_{0}\int\langle\frac{\mathrm{d}\sigma}{\mathrm{d}E_{\chi}}\Big{|}^{V}_{\mathrm{\xi}}\rangle
P_{\mathrm{decay}}\,\mathcal{B}(\chi\to...)\,\epsilon_{E}\,\,\mathrm{d}E_{\chi}\,,$
(43)

where $N_{0}$ incorporates the neutrino flux normalization, the number of
targets and the total time of exposure for a given experiment. Moreover,
$\langle\frac{\mathrm{d}\sigma}{\mathrm{d}E_{\chi}}\big{|}^{V}_{\mathrm{\xi}}\rangle$
represents the differential cross section with respect to the energy of the
SF, averaged over the neutrino flux and $\xi=$ {CE$\nu$NS, E$\nu$ES}. In the
previous expression, $\epsilon_{E}$ denotes the efficiency of electron recoil
reconstruction. We may assume that the final state electrons and positrons
have on average the same energy, $T_{e}\sim m_{\chi}/2$. As a consequence, for
the COHERENT experiment the efficiency has a negligible impact. For DUNE-ND,
an ideal efficiency has been assumed to compute the contours in Fig. 9.
However, we have checked that our present results are not affected
significantly even when a flat, conservative, 15% efficiency is assumed in the
range $50<T_{e}<m_{\chi}^{\mathrm{max}}/2\simeq 150$ MeV. The probability for
$m_{\chi}$ to decay inside the detector is defined as
$P_{\mathrm{decay}}=1-\exp[-L_{\mathrm{det}}\Gamma/(\beta\gamma)]$, with
$L_{\mathrm{det}}$ being the detector length along the beam direction (see
Ref. [93]) and $\beta\gamma$ being the boost factor to go from the center of
mass to the laboratory reference frame. For the DUNE-ND LArTPC detector we set
$L_{\mathrm{det}}=4$ m [12], while for the COHERENT CsI detector we assume
$L_{\mathrm{det}}=11$ cm [128].

In our analysis, when two-body decays are involved, we furthermore assume that
the produced mediator $V$ decays promptly. We do so by introducing a second
probability factor in Eq. (43) and requiring that the sum of the two decay
lengths does not exceed the size of the detector.

A few comments are in order. First, given the SF masses, $m_{\chi}$,
accessible at DM DD experiments, only the two-body decay $\chi\to V\nu_{\ell}$
is kinematically allowed. However, a subsequent decay of $V\to e^{-}e^{+}$ is
not possible since the up-scattered $m_{\chi}$ at XENONnT and LZ are always
smaller than $2m_{e}$. Hence, we conclude that SF decays are not expected to
leave any extra electronic signal at DM DD experiments. Next, concerning
COHERENT, both the two-body and three-body decays leading to a final state
$e^{-}e^{+}$ pair are kinematically open, as explained in [50]. Finally, for
the case of DUNE-ND two-body and three-body decays that lead to either
$e^{-}e^{+}$ or $\mu^{-}\mu^{+}$ final state pairs are possible.

Figure 9: Left: comparison of the 90% C.L. exclusion regions extracted from SF
decays at COHERENT and DUNE-ND. Right: effect of the CP- versus
$\tau$-optimized flux on the sensitivity at DUNE-ND assuming 1 yr of exposure
and the on-axis position. In both figures $r\equiv m_{\chi}/m_{V}$.

Before proceeding to the discussion of our main results, let us stress that,
here, our goal is a simplistic order of magnitude estimation, a detailed Monte
Carlo analysis being beyond the scope of our present study. To this purpose,
we evaluate the 90% C.L. sensitivity contours within a background-free
framework that is based only on the number of decay events, i.e., we require
$N_{\mathrm{decay}}\geq 2.3$ per year. Finally, in our calculation the
decaying particles are assumed to be in the forward direction, which is
nonetheless a reasonable assumption as explained in Ref. [93].

Our results are summarized in Fig. 9 where the various 90% C.L. sensitivity
contours correspond to different values of the ratio $r=m_{V}/m_{\chi}$. In
the left panel of the figure we illustrate a comparison of the exclusion
regions for the cases of COHERENT and DUNE-ND experiments. We observe that
DUNE-ND has the prospect to probe a larger region of the parameter space,
covering SF masses up to about 280 MeV, while in the case of COHERENT the
sensitivity is lost at $\sim$52 MeV, as expected, given the nature of the
neutrino source. It is interesting, however, to note the impact of the
different detection channels, i.e. neutrino-nucleus vs. neutrino-electron
scattering together with the typical energies involved in the two experiments.
The COHERENT experiment is exploiting the SNS neutrino beam for which the
endpoint occurs at $m_{\mu}/2\simeq 52.8$ MeV, while in the case of DUNE the
FNAL beam endpoint is peaking at $\sim$2–3 GeV. (The FNAL neutrino beam
extends up to 125 GeV, however the relevant part of the neutrino spectrum lies
at energies $E_{\nu}\leq 40$ GeV.) The COHERENT experiment is optimized to
detect nuclear scatterings of the heavy CsI target. As a consequence, the
recoil energies involved are tiny, and thus most of the incident neutrino
energy is available to generate SFs777The conservation of energy implies
$E_{\nu}=E_{\chi}-T_{\mathcal{N}}$, with the nuclear recoil energy being
tiny.. On the other hand, for the case of DUNE where the signal is due to
scattering on electrons, most of the incident neutrino energy is carried away
through electron recoils, leaving a small amount of energy for SF production.
This explains that while GeV-order neutrino energies are involved at DUNE, SFs
with a maximum energy of only up to 280 MeV can be probed.

In the right panel of the figure we explore the effect of the obtained
sensitivity when the $\tau$-optimized flux is assumed. By comparing the
endpoint of the contours where a sharp loss of sensitivity occurs, one can
immediately conclude that by using the $\tau$-optimized horn it becomes
feasible to probe slightly larger SF masses $m_{\chi}$. We furthermore find
that there is a crossing point for the CP- vs. $\tau$-optimized driven
contours. Notice that the same behavior of the contours was found for the
scattering-based analysis discussed previously in Fig. 8 as well as in the
analysis of the dipole portal performed in Ref. [93].

Before closing this discussion we find it useful to elaborate on the shape of
the obtained contours. For the cases with $r>1$ (or, equivalently, when
$m_{V}>m_{\chi}$) that correspond to events coming from the three-body decay
$\chi\to\nu_{\ell}e^{-}e^{+}$, the left leg of the contours starts at exactly
$m_{\chi}=2m_{e}$, regardless of the specific value of $r$. On the other hand,
cases with $r<1$ that follow from the two-body decay $\chi\to V\nu_{\ell}$ and
the subsequent decay of the vector mediator $V\to\nu_{\ell}e^{-}e+$, depend
upon the condition $m_{V}\geq 2m_{e}$. Hence the left leg of the respective
contours depends on the value of the ratio $r$.

##  V Conclusions

In this work we have addressed the possible up-scattering of neutrinos into a
new, sterile fermion with a mass in the MeV range. We have considered a
simplified phenomenological scenario that allows for all possible Lorentz-
invariant non-derivative interactions (scalar, pseudoscalar, vector, axial-
vector and tensor) of neutrinos with electrons and first-generation quarks.

Motivated by the rich physics opportunities offered by the next-generation
long-baseline neutrino oscillation experiment DUNE, we have explored the
sensitivity of the DUNE near detector to this scenario. We have simulated
E$\nu$ES data considering both the nominal neutrino mode and a high-energy
configuration. We have found that, while the $\tau$-optimized configuration
allows to reach slightly larger sterile fermion masses, the standard
configuration actually leads to better sensitivities for the scalar and
pseudoscalar interactions. For all other types of interactions, the two
configurations produce similar results. Additionally, we have investigated how
sensitivities change when moving the DUNE near detector off-axis with respect
to the neutrino beam. Although going off-axis allows to reduce backgrounds, as
expected the neutrino flux is also reduced and hence the signal, in turn
translating into poorer sensitivities the farther the detector is moved, for
all interactions. Eventually, the best projected sensitivity is obtained when
the DUNE near detector is on-axis and independently of the beam configuration,
allowing to reach couplings as small as $g_{a}\sim 10^{-4},~{}2\times
10^{-5},~{}10^{-5}$ for $a=$ scalar/pseudoscalar, vector/axial-vector and
tensor interactions respectively, and $m_{\chi}\lesssim$ few/tens of MeV,
after 7 years of data taking. The possible decay of the sterile fermion within
the same detector may provide additional bounds on the parameter space.

Next, we have compared DUNE sensitivities to bounds obtained with existing
data. In this scope, we have analyzed E$\nu$ES data from solar neutrinos at
current DM DD experiments, LZ and XENONnT and CE$\nu$NS data from the most
recent result of the COHERENT-CsI experiment. We have found that all these
facilities dominate the sensitivities in different regions of mediator or
sterile fermion mass. Namely, XENONnT and LZ dominate the constraints at very
small masses; COHERENT bounds exclude smaller couplings in the intermediate
($\sim$ tens of MeV) mass range, while DUNE definitely will provide the most
stringent constraints at heavier masses, being able to probe sterile fermions
as heavy as $\sim 280$ MeV. These considerations hold especially for spin-
dependent interactions, for which CE$\nu$NS data are suppressed. We have
finally explored the possibility of detecting additional signatures from
sterile fermion decays inside the detectors under consideration. While such
decays are not relevant for DM DD experiments, for COHERENT and DUNE-ND our
results lead to the same general conclusions as in the case of up-scattering.

To conclude, our work is meant to show the complementarity of independent
facilities involving different neutrino sources to a phenomenological scenario
in which a neutrino can up-scatter into a sterile fermion. As a follow-up, it
would be interesting to pursue further theoretical work to explore realistic
UV-completions of this phenomenological framework.

## Acknowledgments

This work has been supported by the Spanish grants PID2020-113775GB-I00
(MCIN/AEI/ 10.13039/501100011033) and CNS2023-144124
(MCIN/AEI/10.13039/501100011033 and “Next Generation EU”/PRTR). VDR
acknowledges financial support by the CIDEXG/2022/20 grant (project
“D’AMAGAT”) funded by Generalitat Valenciana. PMC is supported by the grant
CIACIF/2021/281 also funded by Generalitat Valenciana. The work of PM, DKP and
NS was supported by the Hellenic Foundation for Research and Innovation
(H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Post-
Doctoral Researchers” (Project Number: 7036).

##  Appendix A Effective electron charge for Cs, I and Xe

In this Appendix we provide the effective electron charge values for Cs, I and
Xe (see Table 1).

$\boldsymbol{Z^{\mathrm{Cs}}_{\mathrm{eff}}}$ |  $\boldsymbol{T_{e}~{}(}$keV$\mathbf{)}$
---|---
55 |  $T_{e}>$ 35.99
53 | 5.71 $<T_{e}\leq$ | 35.99
51 | 5.36 $<T_{e}\leq$ | 5.71
49 | 5.01 $<T_{e}\leq$ | 5.36
45 | 1.21 $<T_{e}\leq$ | 5.01
43 | 1.07 $<T_{e}\leq$ | 1.21
41 | 1.00 $<T_{e}\leq$ | 1.07
37 | 0.74 $<T_{e}\leq$ | 1.00
33 | 0.73 $<T_{e}\leq$ | 0.74
27 | 0.23 $<T_{e}\leq$ | 0.73
25 | 0.17 $<T_{e}\leq$ | 0.23
23 | 0.16 $<T_{e}\leq$ | 0.17
19 |  $T_{e}\leq$ 0.16

$\boldsymbol{Z^{\mathrm{I}}_{\mathrm{eff}}}$ |  $\boldsymbol{T_{e}~{}(}$keV$\mathbf{)}$
---|---
53 |  $T_{e}>$ 33.17
51 | 5.19 $<T_{e}\leq$ | 33.17
49 | 4.86 $<T_{e}\leq$ | 5.19
47 | 4.56 $<T_{e}\leq$ | 4.86
43 | 1.07 $<T_{e}\leq$ | 4.56
41 | 0.93 $<T_{e}\leq$ | 1.07
39 | 0.88 $<T_{e}\leq$ | 0.93
35 | 0.63 $<T_{e}\leq$ | 0.88
31 | 0.62 $<T_{e}\leq$ | 0.63
25 | 0.19 $<T_{e}\leq$ | 0.62
23 | 0.124 $<T_{e}\leq$ | 0.19
21 | 0.123 $<T_{e}\leq$ | 0.124
17 |  $T_{e}\leq$ 0.123

$\boldsymbol{Z^{\mathrm{Xe}}_{\mathrm{eff}}}$ |  $\boldsymbol{T_{e}~{}(}$keV$\mathbf{)}$
---|---
54 |  $T_{e}>$ 34.76
52 | 5.51 $<T_{e}\leq$ | 34.76
50 | 5.16 $<T_{e}\leq$ | 5.51
48 | 4.84 $<T_{e}\leq$ | 5.16
44 | 1.17 $<T_{e}\leq$ | 4.84
42 | 1.02 $<T_{e}\leq$ | 1.17
40 | 0.96 $<T_{e}\leq$ | 1.02
36 | 0.71 $<T_{e}\leq$ | 0.96
32 | 0.69 $<T_{e}\leq$ | 0.71
26 | 0.23 $<T_{e}\leq$ | 0.69
24 | 0.18 $<T_{e}\leq$ | 0.23
22 | 0.16 $<T_{e}\leq$ | 0.18
18 | 0.074 $<T_{e}\leq$ | 0.16
14 | 0.072 $<T_{e}\leq$ | 0.074
10 | 0.028 $<T_{e}\leq$ | 0.072
4 | 0.013 $<T_{e}\leq$ | 0.028
2 | 0.012 $<T_{e}\leq$ | 0.013
0 |  $T_{e}\leq$ 0.012

Table 1: Effective electron charge for Cs, I [129] and Xe [130] as a function
of the energy deposition $T_{e}$.

## References

  * [1] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Martínez-Miravé, O. Mena, C. A. Ternes, M. Tórtola, and J. W. F. Valle JHEP 02 (2021) 071, arXiv:2006.11237 [hep-ph].
  * [2] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou JHEP 09 (2020) 178, arXiv:2007.14792 [hep-ph].
  * [3] R. N. Mohapatra et al. Rept. Prog. Phys. 70 (2007) 1757–1867, arXiv:hep-ph/0510213.
  * [4] A. Atre, T. Han, S. Pascoli, and B. Zhang JHEP 05 (2009) 030, arXiv:0901.3589 [hep-ph].
  * [5] M. Drewes and B. Garbrecht Nucl. Phys. B 921 (2017) 250–315, arXiv:1502.00477 [hep-ph].
  * [6] A. M. Abdullahi et al. J. Phys. G 50 no. 2, (2023) 020501, arXiv:2203.08039 [hep-ph].
  * [7] C. A. Argüelles et al. Rept. Prog. Phys. 83 no. 12, (2020) 124201, arXiv:1907.08311 [hep-ph].
  * [8] DUNE Collaboration, B. Abi et al. JINST 15 no. 08, (2020) T08008, arXiv:2002.02967 [physics.ins-det].
  * [9] DUNE Collaboration, B. Abi et al. arXiv:2002.03005 [hep-ex].
  * [10] V. De Romeri, E. Fernandez-Martinez, and M. Sorel JHEP 09 (2016) 030, arXiv:1607.00293 [hep-ph].
  * [11] DUNE Collaboration, V. Hewes et al. Instruments 5 no. 4, (2021) 31, arXiv:2103.13910 [physics.ins-det].
  * [12] DUNE Collaboration, B. Abi et al. Eur. Phys. J. C 81 no. 4, (2021) 322, arXiv:2008.12769 [hep-ex].
  * [13] C. Hasnip.  PhD thesis, Oxford U., 2023.
  * [14] C. M. Marshall, K. S. McFarland, and C. Wilkinson Phys. Rev. D 101 no. 3, (2020) 032002, arXiv:1910.10996 [hep-ex].
  * [15] D. Caratelli et al. arXiv:2203.00740 [physics.ins-det].
  * [16] J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly, and J. L. Raaf JHEP 02 (2020) 174, arXiv:1912.07622 [hep-ph].
  * [17] P. Ballett, T. Boschi, and S. Pascoli JHEP 03 (2020) 111, arXiv:1905.00284 [hep-ph].
  * [18] P. Coloma, E. Fernández-Martínez, M. González-López, J. Hernández-García, and Z. Pavlovic Eur. Phys. J. C 81 no. 1, (2021) 78, arXiv:2007.03701 [hep-ph].
  * [19] M. Breitbach, L. Buonocore, C. Frugiuele, J. Kopp, and L. Mittnacht JHEP 01 (2022) 048, arXiv:2102.03383 [hep-ph].
  * [20] V. De Romeri, K. J. Kelly, and P. A. N. Machado Phys. Rev. D 100 no. 9, (2019) 095010, arXiv:1903.10505 [hep-ph].
  * [21] B. Batell et al. in Snowmass 2021.  7, 2022.  arXiv:2207.06898 [hep-ph].
  * [22] V. Brdar, B. Dutta, W. Jang, D. Kim, I. M. Shoemaker, Z. Tabrizi, A. Thompson, and J. Yu Phys. Rev. D 107 no. 5, (2023) 055043, arXiv:2206.06380 [hep-ph].
  * [23] T. Kajita Rev. Mod. Phys. 88 (2016) 030501.
  * [24] A. B. McDonald Rev. Mod. Phys. 88 (2016) 030502.
  * [25] COHERENT Collaboration, D. Akimov et al. Science 357 (2017) 1123–1126, arXiv:1708.01294 [nucl-ex].
  * [26] COHERENT Collaboration, D. Akimov et al. Phys. Rev. Lett. 126 no. 1, (2021) 012002, arXiv:2003.10630 [nucl-ex].
  * [27] COHERENT Collaboration, D. Akimov et al. Phys. Rev. Lett. 129 (2022) 081801, arXiv:2110.07730 [hep-ex].
  * [28] J. Barranco, O. G. Miranda, and T. I. Rashba JHEP 12 (2005) 021, arXiv:hep-ph/0508299.
  * [29] M. Abdullah et al. arXiv:2203.07361 [hep-ph].
  * [30] D. K. Papoulias, T. S. Kosmas, and Y. Kuno Front. in Phys. 7 (2019) 191, arXiv:1911.00916 [hep-ph].
  * [31] V. De Romeri, O. G. Miranda, D. K. Papoulias, G. Sanchez Garcia, M. Tórtola, and J. W. F. Valle JHEP 04 (2023) 035, arXiv:2211.11905 [hep-ph].
  * [32] D. Akimov et al. in 2022 Snowmass Summer Study.  4, 2022.  arXiv:2204.04575 [hep-ex].
  * [33] CCM Collaboration, A. A. Aguilar-Arevalo et al. Phys. Rev. D 106 (2022) 012001, arXiv:2105.14020 [hep-ex].
  * [34] D. Baxter et al. JHEP 02 (2020) 123, arXiv:1911.00762 [physics.ins-det].
  * [35] S. S. Chatterjee, S. Lavignac, O. G. Miranda, and G. Sanchez Garcia Phys. Rev. D 107 no. 5, (2023) 055019, arXiv:2208.11771 [hep-ph].
  * [36] M. Lindner, W. Rodejohann, and X.-J. Xu JHEP 03 (2017) 097, arXiv:1612.04150 [hep-ph].
  * [37] W. Rodejohann, X.-J. Xu, and C. E. Yaguna JHEP 05 (2017) 024, arXiv:1702.05721 [hep-ph].
  * [38] D. Aristizabal Sierra, V. De Romeri, and N. Rojas Phys. Rev. D 98 (2018) 075018, arXiv:1806.07424 [hep-ph].
  * [39] T. S. Kosmas, D. K. Papoulias, M. Tortola, and J. W. F. Valle Phys. Rev. D 96 no. 6, (2017) 063013, arXiv:1703.00054 [hep-ph].
  * [40] C. Blanco, D. Hooper, and P. Machado Phys. Rev. D 101 no. 7, (2020) 075051, arXiv:1901.08094 [hep-ph].
  * [41] O. G. Miranda, D. K. Papoulias, O. Sanders, M. Tórtola, and J. W. F. Valle Phys. Rev. D 102 (2020) 113014, arXiv:2008.02759 [hep-ph].
  * [42] O. G. Miranda, D. K. Papoulias, O. Sanders, M. Tórtola, and J. W. F. Valle JHEP 12 (2021) 191, arXiv:2109.09545 [hep-ph].
  * [43] R. Calabrese, J. Gunn, G. Miele, S. Morisi, S. Roy, and P. Santorelli Phys. Rev. D 107 no. 5, (2023) 055039, arXiv:2212.11210 [hep-ph].
  * [44] V. De Romeri, V. M. Lozano, and G. Sanchez Garcia Phys. Rev. D 109 no. 5, (2024) 055014, arXiv:2307.13790 [hep-ph].
  * [45] D. McKeen and M. Pospelov Phys. Rev. D 82 (2010) 113018, arXiv:1011.3046 [hep-ph].
  * [46] P. D. Bolton, F. F. Deppisch, K. Fridell, J. Harz, C. Hati, and S. Kulkarni Phys. Rev. D 106 no. 3, (2022) 035036, arXiv:2110.02233 [hep-ph].
  * [47] V. Brdar, W. Rodejohann, and X.-J. Xu JHEP 12 (2018) 024, arXiv:1810.03626 [hep-ph].
  * [48] W.-F. Chang and J. Liao Phys. Rev. D 102 no. 7, (2020) 075004, arXiv:2002.10275 [hep-ph].
  * [49] W. Chao, T. Li, J. Liao, and M. Su Phys. Rev. D 104 no. 9, (2021) 095017, arXiv:2108.02341 [hep-ph].
  * [50] P. M. Candela, V. De Romeri, and D. K. Papoulias Phys. Rev. D 108 no. 5, (2023) 055001, arXiv:2305.03341 [hep-ph].
  * [51] D. Alonso-González, D. W. P. Amaral, A. Bariego-Quintana, D. Cerdeño, and M. de los Rios JHEP 12 (2023) 096, arXiv:2307.05176 [hep-ph].
  * [52] M. Schumann J. Phys. G 46 no. 10, (2019) 103003, arXiv:1903.03026 [astro-ph.CO].
  * [53] D. Baxter et al. Eur. Phys. J. C 81 no. 10, (2021) 907, arXiv:2105.00599 [hep-ex].
  * [54] J. Billard et al. Rept. Prog. Phys. 85 no. 5, (2022) 056201, arXiv:2104.07634 [hep-ex].
  * [55] J. Monroe and P. Fisher Phys. Rev. D 76 (2007) 033007, arXiv:0706.3019 [astro-ph].
  * [56] J. D. Vergados and H. Ejiri Nucl. Phys. B 804 (2008) 144–159, arXiv:0805.2583 [hep-ph].
  * [57] L. E. Strigari New J. Phys. 11 (2009) 105011, arXiv:0903.3630 [astro-ph.CO].
  * [58] J. Billard, L. Strigari, and E. Figueroa-Feliciano Phys. Rev. D 89 no. 2, (2014) 023524, arXiv:1307.5458 [hep-ph].
  * [59] C. A. J. O’Hare Phys. Rev. D 94 no. 6, (2016) 063527, arXiv:1604.03858 [astro-ph.CO].
  * [60] C. A. J. O’Hare Phys. Rev. Lett. 127 no. 25, (2021) 251802, arXiv:2109.03116 [hep-ph].
  * [61] D. G. Cerdeño, M. Fairbairn, T. Jubb, P. A. N. Machado, A. C. Vincent, and C. Bœhm JHEP 05 (2016) 118, arXiv:1604.01025 [hep-ph]. [Erratum: JHEP 09, 048 (2016)].
  * [62] B. Dutta, S. Liao, L. E. Strigari, and J. W. Walker Phys. Lett. B 773 (2017) 242–246, arXiv:1705.00661 [hep-ph].
  * [63] G. B. Gelmini, V. Takhistov, and S. J. Witte Phys. Rev. D 99 no. 9, (2019) 093009, arXiv:1812.05550 [hep-ph].
  * [64] R. Essig, M. Sholapurkar, and T.-T. Yu Phys. Rev. D 97 no. 9, (2018) 095029, arXiv:1801.10159 [hep-ph].
  * [65] D. W. P. d. Amaral, D. G. Cerdeno, P. Foldenauer, and E. Reid JHEP 12 (2020) 155, arXiv:2006.11225 [hep-ph].
  * [66] B. Dutta, R. F. Lang, S. Liao, S. Sinha, L. Strigari, and A. Thompson JHEP 09 (2020) 106, arXiv:2002.03066 [hep-ph].
  * [67] D. W. P. Amaral, D. G. Cerdeno, A. Cheek, and P. Foldenauer Eur. Phys. J. C 81 no. 10, (2021) 861, arXiv:2104.03297 [hep-ph].
  * [68] A. de Gouvêa, E. McGinness, I. Martinez-Soler, and Y. F. Perez-Gonzalez Phys. Rev. D 106 no. 9, (2022) 096017, arXiv:2111.02421 [hep-ph].
  * [69] M. Atzori Corona, W. M. Bonivento, M. Cadeddu, N. Cargioli, and F. Dordei Phys. Rev. D 107 no. 5, (2023) 053001, arXiv:2207.05036 [hep-ph].
  * [70] C. Giunti and C. A. Ternes Phys. Rev. D 108 no. 9, (2023) 095044, arXiv:2309.17380 [hep-ph].
  * [71] V. De Romeri, D. K. Papoulias, and C. A. Ternes arXiv:2402.05506 [hep-ph].
  * [72] XENON Collaboration, E. Aprile et al. Eur. Phys. J. C 77 no. 12, (2017) 881, arXiv:1708.07051 [astro-ph.IM].
  * [73] XENON Collaboration, E. Aprile et al. Phys. Rev. Lett. 121 no. 11, (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * [74] XENON Collaboration, E. Aprile et al. Phys. Rev. Lett. 129 no. 16, (2022) 161805, arXiv:2207.11330 [hep-ex].
  * [75] LZ Collaboration, D. S. Akerib et al. arXiv:1509.02910 [physics.ins-det].
  * [76] LZ Collaboration, D. S. Akerib et al. Nucl. Instrum. Meth. A 953 (2020) 163047, arXiv:1910.09124 [physics.ins-det].
  * [77] LZ Collaboration, J. Aalbers et al. Phys. Rev. Lett. 131 no. 4, (2023) 041002, arXiv:2207.03764 [hep-ex].
  * [78] LZ Collaboration, J. Aalbers et al. arXiv:2307.15753 [hep-ex].
  * [79] J. Wyenberg and I. M. Shoemaker Phys. Rev. D 97 no. 11, (2018) 115026, arXiv:1803.08146 [hep-ph].
  * [80] B. Carew, A. R. Caddell, T. N. Maity, and C. A. J. O’Hare arXiv:2312.04303 [hep-ph].
  * [81] D. Aristizabal Sierra, V. De Romeri, L. J. Flores, and D. K. Papoulias JCAP 01 no. 01, (2022) 055, arXiv:2109.03247 [hep-ph].
  * [82] A. Majumdar, D. K. Papoulias, and R. Srivastava Phys. Rev. D 106 no. 1, (2022) 013001, arXiv:2112.03309 [hep-ph].
  * [83] A. N. Khan Phys. Lett. B 837 (2023) 137650, arXiv:2208.02144 [hep-ph].
  * [84] D. W. P. Amaral, D. Cerdeno, A. Cheek, and P. Foldenauer JHEP 07 (2023) 071, arXiv:2302.12846 [hep-ph].
  * [85] T. D. Lee and C.-N. Yang Phys. Rev. 104 (1956) 254–258.
  * [86] P. Melas, D. K. Papoulias, and N. Saoulidou JHEP 07 (2023) 190, arXiv:2303.07094 [hep-ph].
  * [87] G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai Phys. Rev. D 98 no. 11, (2018) 115015, arXiv:1803.03262 [hep-ph].
  * [88] V. Brdar, A. Greljo, J. Kopp, and T. Opferkuch JCAP 01 (2021) 039, arXiv:2007.15563 [hep-ph].
  * [89] P. Vogel and J. Engel Phys. Rev. D 39 (1989) 3378.
  * [90] A. B. Balantekin and N. Vassh Phys. Rev. D 89 no. 7, (2014) 073013, arXiv:1312.6858 [hep-ph].
  * [91] T. Schwetz, A. Zhou, and J.-Y. Zhu JHEP 21 (2020) 200, arXiv:2105.09699 [hep-ph].
  * [92] M. Atkinson, P. Coloma, I. Martinez-Soler, N. Rocco, and I. M. Shoemaker JHEP 04 (2022) 174, arXiv:2105.09357 [hep-ph].
  * [93] M. Ovchynnikov, T. Schwetz, and J.-Y. Zhu Phys. Rev. D 107 no. 5, (2023) 055029, arXiv:2210.13141 [hep-ph].
  * [94] I. M. Shoemaker and J. Wyenberg Phys. Rev. D 99 no. 7, (2019) 075010, arXiv:1811.12435 [hep-ph].
  * [95] I. M. Shoemaker, Y.-D. Tsai, and J. Wyenberg Phys. Rev. D 104 no. 11, (2021) 115026, arXiv:2007.05513 [hep-ph].
  * [96] K. Chakraborty, A. Das, S. Goswami, and S. Roy JHEP 04 (2022) 008, arXiv:2111.08767 [hep-ph].
  * [97] S. Bilmis, I. Turan, T. M. Aliev, M. Deniz, L. Singh, and H. T. Wong Phys. Rev. D 92 no. 3, (2015) 033009, arXiv:1502.07763 [hep-ph].
  * [98] Z. Chen, T. Li, and J. Liao JHEP 05 (2021) 131, arXiv:2102.09784 [hep-ph].
  * [99] T. Schwemberger and T.-T. Yu Phys. Rev. D 106 no. 1, (2022) 015002, arXiv:2202.01254 [hep-ph].
  * [100] M. Demirci and M. F. Mustamin Phys. Rev. D 109 no. 1, (2024) 015021, arXiv:2312.17502 [hep-ph].
  * [101] M. Cirelli, E. Del Nobile, and P. Panci JCAP 10 (2013) 019, arXiv:1307.5955 [hep-ph].
  * [102] E. Del Nobile arXiv:2104.12785 [hep-ph].
  * [103] G. Belanger, F. Boudjema, A. Pukhov, and A. Semenov Comput. Phys. Commun. 185 (2014) 960–985, arXiv:1305.0237 [hep-ph].
  * [104] M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, F. Murgia, A. Prokudin, and S. Melis Nucl. Phys. B Proc. Suppl. 191 (2009) 98–107, arXiv:0812.4366 [hep-ph].
  * [105] T. W. Donnelly and R. D. Peccei Phys. Rept. 50 (1979) 1.
  * [106] M. Hoferichter, J. Menéndez, and A. Schwenk Phys. Rev. D 102 no. 7, (2020) 074018, arXiv:2007.08529 [hep-ph].
  * [107] T. S. Kosmas, V. K. B. Kota, D. K. Papoulias, and R. Sahu Phys. Rev. C 104 no. 6, (2021) 064618, arXiv:2111.08488 [nucl-th].
  * [108] V. Pandey Prog. Part. Nucl. Phys. 134 (2024) 104078, arXiv:2309.07840 [hep-ph].
  * [109] MINERvA Collaboration, J. Park et al. Phys. Rev. D 93 no. 11, (2016) 112007, arXiv:1512.07699 [physics.ins-det].
  * [110] L. Fields, 2017.  https://home.fnal.gov/~ljf26/DUNEFluxes/.
  * [111] https://glaucus.crc.nd.edu/DUNEFluxes/.
  * [112] DUNE Collaboration, R. Acciarri et al. arXiv:1512.06148 [physics.ins-det].
  * [113] P. Machado, H. Schulz, and J. Turner Phys. Rev. D 102 no. 5, (2020) 053010, arXiv:2007.00015 [hep-ph].
  * [114] A. de Gouvea, P. A. N. Machado, Y. F. Perez-Gonzalez, and Z. Tabrizi Phys. Rev. Lett. 125 no. 5, (2020) 051803, arXiv:1912.06658 [hep-ph].
  * [115] DUNE Collaboration, R. Acciarri et al. arXiv:1601.05471 [physics.ins-det].
  * [116] V. Mathur, I. M. Shoemaker, and Z. Tabrizi JHEP 10 (2022) 041, arXiv:2111.14884 [hep-ph].
  * [117] J. Bahcall. http://www.sns.ias.edu/~jnb/SNdata/sndata.html.
  * [118] J. N. Bahcall and R. K. Ulrich Rev. Mod. Phys. 60 (1988) 297–372.
  * [119] J. N. Bahcall Phys. Rev. D 49 (1994) 3923–3945, arXiv:astro-ph/9401024.
  * [120] J. N. Bahcall, E. Lisi, D. E. Alburger, L. De Braeckeleer, S. J. Freedman, and J. Napolitano Phys. Rev. C 54 (1996) 411–422, arXiv:nucl-th/9601044.
  * [121] D. Aristizabal Sierra, N. Rojas, and M. H. G. Tytgat JHEP 03 (2018) 197, arXiv:1712.09667 [hep-ph].
  * [122] XENON Collaboration, E. Aprile et al. Phys. Rev. D 102 no. 7, (2020) 072004, arXiv:2006.09721 [hep-ex].
  * [123] S. K. A., A. Majumdar, D. K. Papoulias, H. Prajapati, and R. Srivastava Phys. Lett. B 839 (2023) 137742, arXiv:2208.06415 [hep-ph].
  * [124] LZ Collaboration, J. Aalbers et al. Phys. Rev. D 108 no. 1, (2023) 012010, arXiv:2211.17120 [hep-ex].
  * [125] COHERENT Collaboration, D. Akimov et al. JINST 17 no. 10, (2022) P10034, arXiv:2111.02477 [physics.ins-det].
  * [126] E. Picciau.  PhD thesis, Cagliari U., 2022.
  * [127] G.-y. Huang, T. Ohlsson, and S. Zhou Phys. Rev. D 97 no. 7, (2018) 075009, arXiv:1712.04792 [hep-ph].
  * [128] A. Konovalov, 2018.  https://sites.duke.edu/coherent/files/2020/06/Konovalov_Experimental_program_2018.pdf.
  * [129] A. Thompson et al., 2009.  https://xdb.lbl.gov/.
  * [130] J.-W. Chen, H.-C. Chi, C. P. Liu, and C.-P. Wu Phys. Lett. B 774 (2017) 656–661, arXiv:1610.04177 [hep-ex].