# COHERENT production of a Dark Fermion

Pablo M. Candela<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<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
Dimitrios K. Papoulias<EMAIL_ADDRESS>Department of Physics,
National and Kapodistrian University of Athens, Zografou Campus GR-15772
Athens, Greece

###### Abstract

We consider the possible production of a new MeV-scale fermion at the COHERENT
experiment. The new fermion, belonging to a dark sector, can be produced
through the up-scattering process of neutrinos off the nuclei and the
electrons of the detector material, via the exchange of a light vector or
scalar mediator. We perform a detailed statistical analysis of the combined
COHERENT CsI and LAr data sets and obtain up-to-date constraints on the
couplings and masses of the dark fermion and mediators. We finally briefly
comment about the stability of the dark fermion.

neutrinos, CE$\nu$NS, COHERENT, light mediators, dark fermion

##  I Introduction

The neutrino physics picture, while not yet complete, is continuously
improving thanks to multiple experimental efforts. On the one hand, there are
oscillation experiments [1, 2] which aim at a precise determination of all
neutrino mixing parameters [3, 4]. On the other hand, experiments exploiting
neutral current neutrino interactions have been proven capable of providing
valuable and complementary information to this picture [5]. Among these, the
recent observation of coherent elastic neutrino-nucleus scattering
(CE$\nu$NS), about four decades after its first theoretical prediction [6],
has opened the window to a plethora of physics opportunities [7, 8]. The
CE$\nu$NS process has been observed by the COHERENT experiment using two
different detectors, one made out of cesium iodide (CsI) [9, 10] and the other
of liquid argon (LAr) [11]. The COHERENT experiment is one of the few that
currently use neutrinos from a $\pi$-DAR (pion decay-at-rest) source. In this
particular case, neutrinos are produced at the Spallation Neutron Source
(SNS). Other CE$\nu$NS experiments exploiting $\pi$-DAR neutrinos are the
COHERENT CAPTAIN-Mills (CCM) at the LANSCE Lujan Center [12] and the planned
European Spallation Source [13]. Several other experiments currently underway
also aim at observing the CE$\nu$NS process. They all rely on a different
neutrino source, namely nuclear reactors, and include facilities like CONNIE
[14], CONUS [15], $\nu$GEN [16], MINER [17], RICOCHET [18], NUCLEUS [19],
TEXONO [20], vIOLETA [21], RED-100 [22], NEON [23], NEWS-G [24] and the
Scintillating Bubble Chamber (SBC) [25]. Let us highlight a recent result
reported by the Dresden-II Collaboration, claiming a very strong preference
for the presence of a CE$\nu$NS component in their data using reactor
antineutrinos [26]. The low-energy tail of more energetic neutrino sources can
also be exploited to study CE$\nu$NS. This is the case, for instance, of the
$\nu$BDX-DRIFT Collaboration which aims at observing CE$\nu$NS induced by
decay-in-flight neutrinos produced at the Long Baseline Neutrino Facility
(LBNF), using a directional time projection chamber [27].

Measurements of CE$\nu$NS provide many opportunities for precision tests of
the Standard Model (SM) parameters [28, 29, 30, 31, 32, 33, 34, 35], but also
for probing the existence of new physics, for example in the form of new
interactions [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 43, 47, 48, 49, 50,
51, 52, 53, 54, 55, 33, 56], nontrivial electromagnetic properties [57, 58,
59, 60, 61, 62, 63, 64, 65] or new states [66, 67, 68, 69, 70, 71, 72, 73, 74,
33, 75]. Interestingly, new states with masses lying at the MeV scale
(dictated by kinematics arguments) can be produced at accelerator-based
facilities like COHERENT. One example is the possible production of a heavy
neutral lepton due to the presence of active-sterile neutrino transition
magnetic moments, via up-scattering processes [76, 71, 33]. In full
generality, not only sterile neutrinos but also other new particles belonging
to a dark sector can be within the reach of CE$\nu$NS experiments. In some
cases, they can constitute all or part of the dark matter (DM) of the
Universe. Indeed, sub-GeV DM particles, while too light to be observed in many
conventional DM detectors, may be produced at accelerators [77, 78, 79].
Various searches for sub-GeV DM particles have been performed both at COHERENT
and at CCM [80, 81, 82, 83, 84], leveraging the characteristic timing
composition of $\pi$-DAR sources which allow to reduce systematic
uncertainties [85, 86]. Among low-mass DM scenarios, a dark photon and a
leptophobic portal have been previously investigated [87, 88, 89, 84, 90].

In this paper, we go beyond the photon- or dark photon-mediated production of
new dark-sector particles, generalizing the up-scattering process to new
vector and scalar mediators. At this scope, we build upon previous analyses
[67, 91, 72, 73] and consider the possible production of a new MeV-scale dark
fermion (DF) at COHERENT. By means of a detailed statistical analysis [33], we
combine the most recent CsI and LAr data sets and infer constraints on the DF
parameter space. Reference [67] already considered the possible production of
a new fermion at CE$\nu$NS experiments, assuming that it couples to neutrinos
and quarks via a singlet scalar. In [91] loop effects in the up-scattering of
neutrinos with CE$\nu$NS, also mediated by a scalar particle, have been
studied. We improve upon these previous results by analyzing new COHERENT data
sets (CsI, 2021 [10] and LAr, 2020 [11]) and by including timing information
in the data analysis. We also consider the possible production of the DF from
the neutrino-electron scattering (ES) on an atomic nucleus, which in some
parts of the parameter space drastically enhances the expected signal at the
CsI detector. Going one step beyond Ref. [67] we further extend the analysis
addressing both (light) scalar and vector mediators. Notice that the possible
production of an exotic fermion from neutrino-electron scattering experiments
had been studied in [73], focusing on effective operators. Reference [72]
instead considered loop effects on vector and axial-vector contributions in
the up-scattering process at COHERENT, using the old CsI data. Finally, the
effective operators leading to the up-scattering interactions of interest here
have been investigated in the context of DM absorption at direct detection
experiments in [92, 93]. In this work, we recompute all the relevant cross
sections and provide their full expressions, to be able to explore the
explicit dependence of both CE$\nu$NS and ES cross sections on the mediator
and the DF masses. We conclude with a brief discussion about the stability of
the DF.

Our paper is organized as follows. In Sec. II we describe how the dark fermion
can be produced through up-scattering, from the interaction of an incoming
neutrino with the nuclei or the atomic electrons of the detector. In Sec. III
we discuss the expected signal and provide details of the statistical
analysis. In Sec. IV we present our results, in terms of exclusion regions in
the dark fermion parameter space. Finally, we summarize our results in Section
V.

##  II Production of a dark fermion at COHERENT

We consider the possible production of a DF $\chi$ at COHERENT, from the up-
scattering of neutrinos produced at the SNS on the nuclei of the detector:

$\nu_{\alpha}\mathcal{N}\to\chi\mathcal{N}\,,\qquad\nu_{\alpha}=\nu_{\mu},\,\overline{\nu}_{\mu},\,\nu_{e}\,.$
(1)

In the case of the COHERENT-CsI detector, ES events are also important as they
can mimic the CE$\nu$NS signal, so we also consider the DF production from
neutrinos scattering off atomic electrons:

$\nu_{\alpha}e^{-}\to\chi
e^{-}\,,\qquad\nu_{\alpha}=\nu_{\mu},\,\overline{\nu}_{\mu},\,\nu_{e}\,.$ (2)

The relevant Feynman diagrams are depicted in Fig. 1.

Figure 1: Feynman diagrams for the $\nu_{\alpha}f\to\chi f$ process
($f=u,\,d,\,e^{-}$), mediated by a vector (on the left) or a scalar (on the
right) particle.

In a simplified phenomenological scenario, we assume that at energies below
the electroweak scale the SM Lagrangian is extended by

$\displaystyle{\mathcal{L}}_{\mathrm{DF}}^{\mathrm{V}}$
$\displaystyle\supseteq$ $\displaystyle
V_{\mu}\,\overline{\chi}\gamma^{\mu}\left(g_{\chi_{L}}P_{L}+g_{\chi_{R}}P_{R}\right)\nu_{\alpha}+V_{\mu}\,\sum_{f}\overline{f}\gamma^{\mu}\left(g_{f_{L}}P_{L}+g_{f_{R}}P_{R}\right)f+\mathrm{H.c.}\,,$
$\displaystyle{\mathcal{L}}_{\mathrm{DF}}^{\mathrm{S}}$
$\displaystyle\supseteq$ $\displaystyle
S\,\overline{\chi}\left(g_{\chi_{L}}P_{L}+g_{\chi_{R}}P_{R}\right)\nu_{\alpha}+S\,\sum_{f}\overline{f}\left(g_{f_{L}}P_{L}+g_{f_{R}}P_{R}\right)f+\mathrm{H.c.}\,,$
(3)

thus encoding new neutrino-DF interactions mediated by either a real scalar
($S$) or a vector ($V_{\mu}$) mediator, also coupling to the charged fermions
of the first family, i.e., $f=u,\,d,\,e^{-}$. The strength of the interaction
is quantified by the couplings $g_{\chi_{L}}$ and $g_{\chi_{R}}$ for the
vertex of $\chi$ with the mediator and the neutrino (we assume this coupling
to be identical for all neutrino flavors), and $g_{f_{L}}$ and $g_{f_{R}}$ for
the vertex of $f$ and the mediator. The subscripts $L$ and $R$ in the
couplings denote their left- and right-handed components, while $P_{L}$ and
$P_{R}$ are the left- and right-handed chiral projectors, respectively. Given
the nature of neutrinos in the SM, we notice that $g_{\chi_{R}}=0$. While in
the present work we decide to remain agnostic of the origin of such
interactions, we refer the reader to [94, 95, 55] for a discussion on possible
UV-completed models leading to such low-scale interactions. Notice also that
we consider the most general case of light mediators, which can give rise to
interesting spectral features at CE$\nu$NS experiments like COHERENT. As a
consequence, a dependence on the mediator mass and the momentum transfer will
appear in the up-scattering cross sections, from the propagator. We should
mention that several constraints may apply to such light-mediator scenarios,
including laboratory limits from fixed-target and beam-dump experiments, rare
decays and accelerator data (see for instance [96, 97]). Most of these
constraints require the new mediator to couple to charged leptons and depend
on the Lorentz structure of the coupling. Scalar interactions are also
constrained by neutrino masses, to which they can contribute at the loop
level. Moreover, light new particles coupling to neutrinos and matter fields
can also affect the dynamics of stellar cooling and supernovae [94, 48] and
Big Bang Nucleosynthesis (BBN).

Given the Lagrangians in Eq. (II) and the diagrams shown in Fig. 1, we can
compute the relevant $\nu_{\alpha}\mathcal{N}\to\chi\mathcal{N}$ cross
sections for the vector (V) and scalar (S) mediators, in terms of the nuclear
recoil energy $E_{\mathrm{nr}}$. Note that the following expressions hold for
both neutrinos and antineutrinos 111We have checked that the squared amplitude
is the same in both cases, after applying the sum over spins. and are
therefore applicable to all components of the SNS neutrino flux. The up-
scattering cross section for the vector case reads

$\displaystyle\frac{\mathrm{d}\sigma_{\nu_{\alpha}\mathcal{N}}}{\mathrm{d}E_{\mathrm{nr}}}\Big{|}^{\mathrm{V}}_{\mathrm{CE\nu
NS}}(E_{\nu},E_{\mathrm{nr}})=$
$\displaystyle\,\dfrac{9\,m_{\mathcal{N}}}{4\pi\left(m_{V}^{2}+2m_{\mathcal{N}}E_{\mathrm{nr}}\right)^{2}}F_{W}^{2}(\left|\mathbf{q}\right|^{2})A^{2}g_{V}^{4}$
(4)
$\displaystyle\times\left[\left(2-\dfrac{m_{\mathcal{N}}E_{\mathrm{nr}}}{E_{\nu}^{2}}-\dfrac{2E_{\mathrm{nr}}}{E_{\nu}}+\dfrac{E_{\mathrm{nr}}^{2}}{E_{\nu}^{2}}\right)-\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\left(1+\dfrac{2E_{\nu}}{m_{\mathcal{N}}}-\dfrac{E_{\mathrm{nr}}}{m_{\mathcal{N}}}\right)\right]\,,\,$

where $E_{\nu}$ is the incident neutrino energy, $m_{\chi}$ is the mass of the
DF, $m_{V}$ is the vector mediator mass, $m_{\mathcal{N}}$ is the nuclear mass
and $A$ is the mass number. The vector coupling has been defined as
$g_{V}\equiv\sqrt{g_{\chi_{L}}g_{f}}$ (we refer to Appendix A for more
details). In the present work, nuclear-physics effects are incorporated
through the nuclear form factor $F_{W}(\left|\mathbf{q}\right|^{2})$, for
which we adopt the Klein-Nystrand parametrization, defined as

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

where $j_{1}(x)=\sin(x)/x^{2}-\cos(x)/x$ is the spherical Bessel function of
order one,
$\left|\mathbf{q}\right|\approx\sqrt{2m_{\mathcal{N}}E_{\mathrm{nr}}}$
represents the magnitude of the three-momentum transfer,
$R_{A}=1.23\,A^{1/3}~{}\mathrm{fm}$ is the nuclear radius and
$a_{k}=0.7~{}\mathrm{fm}$ is the Yukawa potential range.

The scalar mediator cross section, on the other hand, is given by

$\displaystyle\frac{\mathrm{d}\sigma_{\nu_{\alpha}\mathcal{N}}}{\mathrm{d}E_{\mathrm{nr}}}\Big{|}_{\mathrm{CE\nu
NS}}^{\mathrm{S}}(E_{\nu},E_{\mathrm{nr}})=$
$\displaystyle\,\dfrac{m_{\mathcal{N}}}{8\pi\left(m_{S}^{2}+2m_{\mathcal{N}}E_{\mathrm{nr}}\right)^{2}}C_{S}^{4}F_{W}^{2}(\left|\mathbf{q}\right|^{2})$
(6)
$\displaystyle\times\left(2+\dfrac{E_{\mathrm{nr}}}{m_{\mathcal{N}}}\right)\left(\dfrac{m_{\mathcal{N}}E_{\mathrm{nr}}}{E_{\nu}^{2}}+\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\right)\,.\,$

where $m_{S}$ is is the scalar mediator mass and $C_{S}$ is the scalar
coupling. It is defined as

$C_{S}^{2}\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),$
(7)

where $g_{S}\equiv\sqrt{g_{\chi_{L}}g_{f}}$, $Z$ is the number of protons of
the nucleus and $N=A-Z$ is the number of neutrons, $m_{p}$ and $m_{n}$ are the
proton and neutron masses, respectively, and $m_{q}$ is the quark mass. Here,
$f_{T_{q}}^{(p)}$ and $f_{T_{q}}^{(n)}$ indicate the quark-mass contributions
to the nucleon (proton and neutron) mass. While observing that the
determination of these values is still uncertain [98, 99], in the following we
fix [99]

$\displaystyle f_{T_{u}}^{(p)}$ $\displaystyle=0.026\,,\qquad
f_{T_{u}}^{(n)}=0.018\,,$ (8) $\displaystyle f_{T_{d}}^{(p)}$
$\displaystyle=0.038\,,\qquad f_{T_{d}}^{(n)}=0.056\,.$ (9)

We have checked that using a different value [100] would change the value of
$C_{S}^{2}$ by a factor of $\sim 30\%$, at most, however the sensitivity on
$g_{S}$ remains essentially unaffected.

Possible axial-vector and pseudoscalar contributions are in general suppressed
by the nuclear spin as detailed in Ref. [42]. Let us also note that in the
present analysis, the left- and right-handed DF couplings are taken to be
equal in Eq. (II), thus axial-vector and pseudoscalar contributions are
vanishing (see Appendix A).

Next, we provide expressions for the production of a DF via ES on an atomic
nucleus $\mathcal{A}$ containing $Z$ protons. Since the electrons are bound in
$\mathcal{A}$, for a given energy deposition $E_{\mathrm{er}}$ the relevant
cross sections are expressed as the product of the free ES cross section times
the effective charge $Z_{\text{eff}}^{\mathcal{A}}(E_{\mathrm{er}})$, as

$\displaystyle\frac{\mathrm{d}\sigma_{\nu_{\alpha}\mathcal{A}}}{\mathrm{d}E_{\mathrm{er}}}\Big{|}_{\mathrm{ES}}^{\mathrm{V}}(E_{\nu},E_{\mathrm{er}})=$
$\displaystyle\,Z_{\mathrm{eff}}^{\mathcal{A}}\left(E_{\mathrm{er}}\right)\dfrac{m_{e}g_{V}^{4}}{4\pi\left(m_{V}^{2}+2m_{e}E_{\mathrm{er}}\right)^{2}}$
(10)
$\displaystyle\times\left[\left(2-\dfrac{m_{e}E_{\mathrm{er}}}{E_{\nu}^{2}}-\dfrac{2E_{\mathrm{er}}}{E_{\nu}}+\dfrac{E_{\mathrm{er}}^{2}}{E_{\nu}^{2}}\right)-\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}\left(1+\dfrac{2E_{\nu}}{m_{e}}-\dfrac{E_{\mathrm{er}}}{m_{e}}\right)\right]\,,\,$

$\frac{\mathrm{d}\sigma_{\nu_{\alpha}\mathcal{A}}}{\mathrm{d}E_{\mathrm{er}}}\Big{|}_{\mathrm{ES}}^{\mathrm{S}}(E_{\nu},E_{\mathrm{er}})=\,Z_{\mathrm{eff}}^{\mathcal{A}}\left(E_{\mathrm{er}}\right)\dfrac{m_{e}g_{S}^{4}}{8\pi\left(m_{S}^{2}+2m_{e}E_{\mathrm{er}}\right)^{2}}\left(2+\dfrac{E_{\mathrm{er}}}{m_{e}}\right)\left(\dfrac{m_{\chi}^{2}}{2E_{\nu}^{2}}+\dfrac{m_{e}E_{\mathrm{er}}}{E_{\nu}^{2}}\right)\,,\,$
(11)

where $m_{e}$ is the electron mass. The effective charges
$Z_{\text{eff}}^{\mathcal{A}}(E_{\mathrm{er}})$ for Cs and I isotopes are
given in Appendix B.

A few comments are in order. First, let us note that in the limit $m_{\chi}\to
0$ one recovers exactly the expressions for CE$\nu$NS and ES with new light
mediators (see for instance [33, 101, 102]). Second, throughout the manuscript
we will assume that the new mediators couple equally to either $u,\,d$ or
$e^{-}$. The hypothesis of universal couplings reduces the degrees of freedom
simplifying the analysis, thus allowing us to express the results in terms of
a single coupling, namely $g_{\mathrm{S}}$ or $g_{\mathrm{V}}$, and the
corresponding mediator and DF masses. While this model is anomalous, it can be
made anomaly-free by introducing new particles, which may also belong to the
dark sector [55]. Moreover, while the scalar and vector DF cross sections can
in principle involve flavor-dependent terms, in the present work we restrict
to DF-production cross sections which are flavor-blind. It is finally
important to note that the DF mass is constrained from above:
$m_{\chi}\leq\sqrt{m\left(m_{\mu}+m\right)}-m$, where $m$ is the mass of the
fixed target ($m_{\mathcal{N}}$ or $m_{e}$) and $m_{\mu}$ is the muon mass.
This bound is dictated by kinematics and implies that the process shown in
Fig. 1 will not occur if $m_{\chi}$ is above that limit. For the sake of
completeness, let us conclude mentioning that the SM CE$\nu$NS and ES cross
sections are well-known results in the literature and can be found in, e.g.,
Refs. [7] and [103]. In the SM expressions for both CE$\nu$NS and ES, the
$V-A$ interference term appears with a different sign for neutrino or
antineutrino scattering. However, in the case of CE$\nu$NS, the axial
contributions are tiny for the nuclei of interest in the present work and
hence neglected in our calculations 222In the SM case, angular momentum
conservation implies that the axial-vector contribution to the CE$\nu$NS cross
section vanishes for even-even nuclei, while for the other cases the
contribution is tiny, especially for heavy nuclei [36].. Moreover, while in
the SM the CE$\nu$NS cross section is flavor-blind, the SM ES cross section is
different for $\nu_{e}-e^{-}$ compared to $\nu_{\mu,\tau}-e^{-}$ scattering.
The reason is that the former acquires contributions from both charged- and
neutral-currents, while for the latter only neutral-currents are relevant.
Before closing our discussion let us clarify that there is no interference
between the vector/scalar CE$\nu$NS cross sections and the SM one. The same
holds for the ES contribution.

##  III Data analysis

We now proceed to analyze the most recent COHERENT data sets, from the CsI
[10] and LAr [104] detectors. In each case, we take into account energy and
timing information as well as all relevant systematic effects, as done in Ref.
[33]. Moreover, in the CsI case, we also include possible contributions from
ES events [105, 106, 33], which turn out to be relevant especially in the
vector-mediated scenario. In the analysis of the LAr data set, the measurement
of the ratio of the integrated photomultiplier amplitude and the total
amplitude in the first 90 ns ($F_{90}$) allows to distinguish between
CE$\nu$NS-induced nuclear recoils and ES-induced electron recoils and
therefore we do not include the ES contribution.

For the COHERENT-CsI (COHERENT-LAr) analysis, we assume a detector mass
$m_{\rm det}=14.6~{}(24)$ kg located at a distance $L=19.3~{}(27.5)$ m from
the SNS source. For the calculation of CE$\nu$NS events,
$N_{ij}^{\mathrm{CE}\nu\mathrm{NS}}(\mathcal{N})$, we follow the procedure
described in Ref. [33]. The expected number of DF or SM events, on a nuclear
target $\mathcal{N}$, per neutrino flavor, $\nu_{\alpha}$, and in each nuclear
recoil energy bin $i$ is given by

$\displaystyle N_{i,\nu_{\alpha},\mathrm{\kappa}}^{\mathrm{CE\nu
NS}}(\mathcal{N})=$
$\displaystyle\,N_{\mathrm{target}}\int_{E_{\mathrm{nr}}^{i}}^{E_{\mathrm{nr}}^{i+1}}\mathrm{d}E_{\mathrm{nr}}\,\epsilon_{E}(E_{\mathrm{nr}})\int_{E^{\prime\mathrm{min}}_{\mathrm{nr}}}^{E^{\prime\text{max}}_{\text{nr}}}\mathrm{d}E^{\prime}_{\text{nr}}\,\mathcal{R}(E_{\text{nr}},E^{\prime}_{\text{nr}})$
$\displaystyle\times\int_{E_{\nu}^{\text{min}}(E^{\prime}_{\text{nr}})}^{E_{\nu}^{\text{max}}}\mathrm{d}E_{\nu}\,\frac{\mathrm{d}N_{\nu_{\alpha}}}{\mathrm{d}E_{\nu}}(E_{\nu})\frac{\mathrm{d}\sigma_{\nu_{\alpha}\mathcal{N}}}{\mathrm{d}E^{\prime}_{\mathrm{nr}}}\Big{|}_{\mathrm{CE\nu
NS}}^{\mathrm{\kappa}}(E_{\nu},E^{\prime}_{\mathrm{nr}}),$ (12)

where
$N_{\mathrm{target}}=N_{\mathrm{A}}m_{\mathrm{det}}/M_{\mathrm{\mathrm{target}}}$
is the number of target atoms in the detector, with $N_{\mathrm{A}}$ being the
Avogadro’s constant and $M_{\mathrm{\mathrm{target}}}$ the molar mass of the
detector material. Index $\mathrm{\kappa}$ accounts for the different
interactions, namely $\mathrm{\kappa=\\{SM,S,V\\}}$. The three components of
the differential $\pi$-DAR neutrino flux produced at the SNS read

$\displaystyle\frac{\mathrm{d}N_{\nu_{\mu}}}{\mathrm{d}E_{\nu}}(E_{\nu})$
$\displaystyle=\eta\,\delta\left(E_{\nu}-\frac{m_{\pi}^{2}-m_{\mu}^{2}}{2m_{\pi}}\right)\quad$
$\displaystyle(\text{prompt})\,,$ (13)
$\displaystyle\frac{\mathrm{d}N_{\bar{\nu}_{\mu}}}{\mathrm{d}E_{\nu}}(E_{\nu})$
$\displaystyle=\eta\frac{64E^{2}_{\nu}}{m_{\mu}^{3}}\left(\frac{3}{4}-\frac{E_{\nu}}{m_{\mu}}\right)\quad$
$\displaystyle(\text{delayed})\,,$
$\displaystyle\frac{\mathrm{d}N_{\nu_{e}}}{\mathrm{d}E_{\nu}}(E_{\nu})$
$\displaystyle=\eta\frac{192E^{2}_{\nu}}{m_{\mu}^{3}}\left(\frac{1}{2}-\frac{E_{\nu}}{m_{\mu}}\right)\quad$
$\displaystyle(\text{delayed})\,,$

($m_{\pi}$ is the pion mass) and are normalized to
$\eta=rN_{\mathrm{POT}}/4\pi L^{2}$, where $r$ denotes the number of neutrinos
per flavor produced for each proton on target (POT). For the CsI detector
$r=0.0848$ and $N_{\mathrm{POT}}=3.198\times 10^{23}$, while for the LAr
detector $r=0.009$ and $N_{\mathrm{POT}}=1.38\times 10^{23}$. The integration
limits read

$\displaystyle
E_{\nu}^{\mathrm{min}}(E^{\prime}_{\mathrm{nr}})=\dfrac{2m_{\mathcal{N}}\left(E^{\prime}_{\mathrm{nr}}\right)^{2}+m_{\chi}^{2}E^{\prime}_{\mathrm{nr}}}{4m_{\mathcal{N}}E^{\prime}_{\mathrm{nr}}}$
$\displaystyle\hskip
50.00008pt+\dfrac{\sqrt{\left(2m_{\mathcal{N}}\left(E^{\prime}_{\mathrm{nr}}\right)^{2}+m_{\chi}^{2}E^{\prime}_{\mathrm{nr}}\right)^{2}+2m_{\mathcal{N}}E^{\prime}_{\mathrm{nr}}\left(m_{\chi}^{4}+4m_{\mathcal{N}}^{2}\left(E^{\prime}_{\mathrm{nr}}\right)^{2}+4m_{\mathcal{N}}m_{\chi}^{2}E^{\prime}_{\mathrm{nr}}\right)}}{4m_{\mathcal{N}}E^{\prime}_{\mathrm{nr}}}$
$\displaystyle\hskip
50.00008pt\approx\sqrt{\dfrac{m_{\mathcal{N}}E^{\prime}_{\mathrm{nr}}}{2}}\left(1+\dfrac{m_{\chi}^{2}}{2m_{\mathcal{N}}E^{\prime}_{\mathrm{nr}}}\right),$
$\displaystyle E_{\nu}^{\mathrm{max}}=\dfrac{m_{\mu}}{2},$ $\displaystyle
E^{\prime\begin{subarray}{c}\scriptscriptstyle\mathrm{min}\\\
\scriptscriptstyle\mathrm{(max)}\end{subarray}}_{\mathrm{nr}}=\dfrac{2m_{\mathcal{N}}{(E_{\nu}^{\mathrm{max}})}^{2}-m_{\chi}^{2}\left(E_{\nu}^{\mathrm{max}}+m_{\mathcal{N}}\right){\begin{subarray}{c}-\\\
(+)\end{subarray}}E_{\nu}^{\mathrm{max}}\sqrt{4m_{\mathcal{N}}^{2}{(E_{\nu}^{\mathrm{max}})}^{2}-4m_{\mathcal{N}}m_{\chi}^{2}\left(E_{\nu}^{\mathrm{max}}+m_{\mathcal{N}}\right)+m_{\chi}^{4}}}{2m_{\mathcal{N}}\left(2E_{\nu}^{\mathrm{max}}+m_{\mathcal{N}}\right)}.$
(14)

Notice that when $m_{\chi}\to 0$ the integration limits relevant to the SM
case are recovered. The remaining detector-specific quantities, namely the
energy resolution function $\mathcal{R}(E_{\text{nr}},E^{\prime}_{\text{nr}})$
relating the true nuclear recoil energy ($E^{\prime}_{\text{nr}}$) with the
reconstructed one ($E_{\text{nr}}$) as well as the energy-dependent detector
efficiency $\epsilon_{E}(E_{\mathrm{nr}})$ are explained in detail in Appendix
C.

We further include timing information in the analysis, by distributing the
predicted $N_{i,\nu_{\alpha},\mathrm{\kappa}}^{\mathrm{CE\nu
NS}}(\mathcal{N})$ in each time bin $j$. We use the time distributions
$\mathcal{P}^{\nu_{\alpha}}_{T}(t_{\mathrm{rec}})$ of
$\nu_{\alpha}=\nu_{\mu},\,\overline{\nu}_{\mu},\,\nu_{e}$ given in [107, 10],
normalized to 6 $\mathrm{\mu s}$. Finally, the predicted event number, per
observed nuclear recoil energy and time bins $i,j$ is given by

$N_{ij,\mathrm{\kappa}}^{\mathrm{CE\nu
NS}}(\mathcal{N})=\sum_{\nu_{\alpha}=\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_{\alpha}}_{T}(t_{\mathrm{rec}},\alpha_{6})\epsilon_{T}(t_{\mathrm{rec}})N_{i,\nu_{\alpha},\mathrm{\kappa}}^{\mathrm{CE\nu
NS}}(\mathcal{N}),$ (15)

where $\epsilon_{T}(t_{\mathrm{rec}})$ is the time-dependent efficiency.
Notice that we include an additional nuisance parameter on the beam timing,
i.e. $\alpha_{6}$ [33] (for details, see the discussion below and Appendix C).
Then, the total CE$\nu$NS event rate is simply given by

$\displaystyle N_{ij}^{\mathrm{CE\nu NS}}=$
$\displaystyle\sum_{\mathrm{\kappa=SM,V}}\sum_{\mathcal{N}}N_{ij,\mathrm{\kappa}}^{\mathrm{CE\nu
NS}}(\mathcal{N})\qquad$ (vector mediator), (16) $\displaystyle
N_{ij}^{\mathrm{CE\nu NS}}=$
$\displaystyle\sum_{\mathrm{\kappa=SM,S}}\sum_{\mathcal{N}}N_{ij,\mathrm{\kappa}}^{\mathrm{CE\nu
NS}}(\mathcal{N})\qquad$ (scalar mediator),

where the inner sum runs over the Cs and I isotopes for the case of CsI
detector, while for the case of LAr the inner sum is trivially dropped.

Similarly to the CE$\nu$NS case described above, the expected number of DF (or
SM) ES events, on an atomic nucleus $\mathcal{A}$, per neutrino flavor,
$\nu_{\alpha}$, and in each electron recoil energy bin $i$ is given by

$\displaystyle N_{i,\nu_{\alpha},\mathrm{\kappa}}^{\mathrm{ES}}(\mathcal{A})=$
$\displaystyle\,N_{\mathrm{target}}\int_{E_{\mathrm{er}}^{i}}^{E_{\mathrm{er}}^{i+1}}\mathrm{d}E_{\mathrm{er}}\,\epsilon_{E}(E_{\mathrm{er}})\int_{E^{\prime\mathrm{min}}_{\mathrm{nr}}}^{E^{\prime\text{max}}_{\text{er}}}\mathrm{d}E^{\prime}_{\text{er}}\,\mathcal{R}(E_{\text{er}},E^{\prime}_{\text{er}})$
$\displaystyle\times\int_{E_{\nu}^{\text{min}}(E^{\prime}_{\text{er}})}^{E_{\nu}^{\text{max}}}\mathrm{d}E_{\nu}\,\frac{\mathrm{d}N_{\nu_{\alpha}}}{\mathrm{d}E_{\nu}}(E_{\nu})\frac{\mathrm{d}\sigma_{\nu_{\alpha}\mathcal{A}}}{\mathrm{d}E^{\prime}_{\mathrm{er}}}\Big{|}_{\mathrm{ES}}^{\mathrm{\kappa}}(E_{\nu},E^{\prime}_{\mathrm{er}})\,.$
(17)

The nuclear recoil energy and the corresponding electron-equivalent energy are
related via the quenching factor (QF) (see Appendix C for details), while the
integration limits are given by Eq. (III) via the substitutions
$m_{\mathcal{N}}\to m_{e}$ and $E^{\prime}_{\text{nr}}\to
E^{\prime}_{\text{er}}$ 333The approximated expression given for
$E_{\nu}^{\text{min}}(E^{\prime}_{\text{nr}})$ in Eq. (III) is only valid for
CE$\nu$NS, not for ES.. Then, the corresponding event rate per observed
electron recoil energy $i$ and time bin $j$ is given by

$N_{ij,\mathrm{\kappa}}^{\mathrm{ES}}(\mathcal{A})=\sum_{\nu_{\alpha}=\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_{\alpha}}_{T}(t_{\mathrm{rec}},\alpha_{6})\epsilon_{T}(t_{\mathrm{rec}})N_{i,\nu_{\alpha},\mathrm{\kappa}}^{\mathrm{ES}}(\mathcal{A})\,.$
(18)

Finally, the total ES event rate is simply given by

$\displaystyle N_{ij}^{\mathrm{ES}}=$
$\displaystyle\sum_{\mathrm{\kappa=SM,V}}\sum_{\mathcal{N}}N_{ij,\mathrm{\kappa}}^{\mathrm{ES}}(\mathcal{A})\qquad$
(vector mediator), (19) $\displaystyle N_{ij}^{\mathrm{ES}}=$
$\displaystyle\sum_{\mathrm{\kappa=SM,S}}\sum_{\mathcal{N}}N_{ij,\mathrm{\kappa}}^{\mathrm{ES}}(\mathcal{A})\qquad$
(scalar mediator),

see the discussion below Eq. (16).

For the statistical analysis of COHERENT data we rely on the following
Poissonian least-squares function [33]

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

The predicted number of events, including CE$\nu$NS, DF and background events,
is defined as

$\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{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}}\,.$
(21)

These expressions involve several nuisances ($\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 the backgrounds, Beam
Related Neutrons (BRN), Neutrino Induced Neutrons (NIN) and Steady State
Background (SSB), respectively. Finally, $\sigma_{5}=3.8\%$ is associated to
the QF. Notice that the number of events $N_{ij}^{\mathrm{CE\nu NS}}$,
$N_{ij}^{\mathrm{DF}}$ and $N_{ij}^{\mathrm{ES}}$ also include nuisance
parameters: $\alpha_{4}$, which enters the nuclear form factor and thus
affects only the CE$\nu$NS number of events 444This is done by introducing the
nuisance parameter $\alpha_{4}$ in the nuclear radius entering Eq. (5) such
that $R_{A}=1.23\,A^{1/3}(1+\alpha_{4})$., with $\sigma_{4}=5\%$; $\alpha_{6}$
accommodates the uncertainty in beam timing with no prior assigned, while
$\alpha_{7}$ allows for deviations of the uncertainty in the CE$\nu$NS
efficiency. We refer the reader to Ref. [33] for further details about the
statistical analysis.

For the COHERENT-LAr $\chi^{2}$ analysis we instead consider the following
Gaussian least-squares function, based on [55, 33]

$\chi^{2}_{\mathrm{LAr}}=\sum_{i=1}^{12}\sum_{j=1}^{10}\left(\frac{N^{\mathrm{th}}_{ij}-N_{ij}^{\text{exp}}}{\sigma_{ij}}\right)^{2}+\sum_{k=0,3,4,8}\left(\dfrac{\beta_{k}}{\sigma_{k}}\right)^{2}+\sum_{k=1,2,5,6,7}\left(\beta_{k}\right)^{2}\,,$
(22)

with the experimental uncertainty being
$\sigma_{ij}^{2}=N_{ij}^{\text{exp}}+N_{ij}^{\mathrm{SSB}}/5$, while the
theoretical signal is calculated as

$\displaystyle N^{\mathrm{th}}_{ij}=$
$\displaystyle(1+\beta_{0}+\beta_{1}\Delta_{\mathrm{CE\nu
NS}}^{F_{90+}}+\beta_{1}\Delta_{\mathrm{CE\nu
NS}}^{F_{90-}}+\beta_{2}\Delta_{\mathrm{CE\nu
NS}}^{\mathrm{t_{trig}}})N_{ij}^{\mathrm{CE\nu
NS}}+(1+\beta_{3})N_{ij}^{\mathrm{SSB}}$ (23)
$\displaystyle+(1+\beta_{4}+\beta_{5}\Delta_{\mathrm{pBRN}}^{E_{+}}+\beta_{5}\Delta_{\mathrm{pBRN}}^{E_{-}}+\beta_{6}\Delta_{\mathrm{pBRN}}^{t_{\text{trig}}^{+}}+\beta_{6}\Delta_{\mathrm{pBRN}}^{t_{\text{trig}}^{-}}+\beta_{7}\Delta_{\mathrm{pBRN}}^{t_{\text{trig}}^{\text{w}}})N_{ij}^{\mathrm{pBRN}}$
$\displaystyle+(1+\beta_{8})N_{ij}^{\mathrm{dBRN}}\,.$

Here, the nuisance parameters $\beta_{0},~{}\beta_{3},~{}\beta_{4}$ and
$\beta_{8}$ account for the normalization uncertainties of CE$\nu$NS, SS,
prompt BRN (pBRN) and delayed BRN (dBRN) rates, with
$\\{\sigma_{0},~{}\sigma_{3},~{}\sigma_{4},~{}\sigma_{8}\\}$ =
$\\{0.13,~{}0.0079,~{}0.32,~{}1.0\\}$ [11]. Relevant systematic effects
affecting the shape uncertainties of the CE$\nu$NS and pBRN rates are also
taken into account through the nuisance parameters
$\beta_{1},\beta_{2},\beta_{5},\beta_{6}$ and $\beta_{7}$. In particular,
$\beta_{1}$ and $\beta_{2}$ account for the uncertainty on the CE$\nu$NS shape
due to existing systematic uncertainties on the $\pm 1\sigma$ energy
distributions of the $F_{90}$ parameter ($\Delta_{\mathrm{CE\nu
NS}}^{F_{90\pm}}$) and due to the mean time to trigger distribution
($\Delta_{\mathrm{CE\nu NS}}^{\mathrm{t_{trig}}}$), respectively. On the other
hand $\beta_{5}$, $\beta_{6}$ and $\beta_{7}$ are introduced to quantify the
pBRN shape uncertainty due to the corresponding uncertainty on the $\pm
1\sigma$ energy distributions ($\Delta_{\mathrm{pBRN}}^{E_{\pm}}$), the $\pm
1\sigma$ mean time to trigger distributions
($\Delta_{\mathrm{pBRN}}^{t_{\text{trig}}^{\pm}}$) as well as the trigger
width distribution ($\Delta_{\mathrm{pBRN}}^{t_{\text{trig}}^{\text{w}}}$).
The latter five distributions are defined as

$\Delta_{\lambda}^{\xi_{\lambda}}=\frac{N_{ij}^{\lambda,\xi_{\lambda}}-N_{ij}^{\lambda,\mathrm{CV}}}{N_{ij}^{\lambda,\mathrm{CV}}}\,,$
(24)

with $\lambda=$ {CE$\nu$NS, pBRN} and $\xi_{\lambda}$ corresponding to the
different source uncertainties affecting the CE$\nu$NS or pBRN shapes. Ιn the
definition of Eq. (24) the superscript “CV” denotes the central values of the
CE$\nu$NS or pBRN distributions, available in the COHERENT-LAr data release
[104]. Before closing this discussion, let us note that the $\beta_{0}$
component includes several uncertainties accounting for the flux (10%),
efficiency (3.6%), energy calibration (0.8%), the calibration of the pulse-
shape discrimination parameter $F_{90}$ (7.8%), QF (1%), and nuclear form
factor (2%) [11].

##  IV Results

In this section, we present the obtained constraints on the DF parameter space
from the combined analysis of COHERENT data from the CsI [10] and LAr [11]
detectors. We start discussing the results for the vector-mediator case and
then we proceed to the scalar one. We further compare our results with
existing limits in the literature (when applicable). We finally briefly
comment about the stability of the DF and a possible connection to the dark
matter.

###  IV.1 Vector mediator

In the following, we refer to the Lagrangian given in Eq. (II) and to the
cross sections given in Eqs. (4) and (10). Thanks to the redefinition of the
coupling $g_{V}\equiv\sqrt{g_{\chi_{L}}g_{f}}$, with $f=u,\,d$ for CE$\nu$NS
and $f=e^{-}$ for ES, we can express our results in terms of the relevant
parameters, namely the DF mass ($m_{\chi}$), the mediator mass ($m_{V}$) and
the coupling ($g_{V}$). This simplification is based on the assumption that
the vector mediator has universal couplings to both quarks and electrons.

Figure 2: $90\%$ C.L. (2 d.o.f.) exclusion regions in the $m_{V}-g_{V}$ plane,
for the vector-mediator case. Different fixed values of $m_{\chi}$ are
considered. The top row shows the analysis done with the CsI (left) and LAr
(right) COHERENT data, while the bottom row shows the combined analysis (left)
and a comparison between the different analyses for the most constraining case
(right).

Figure 3: $90\%$ C.L. (2 d.o.f.) exclusion regions in the $m_{\chi}-g_{V}$
plane, for the vector-mediator case. Different fixed values of $m_{V}$ are
considered.

Figure 4: $90\%$ C.L. (2 d.o.f.) exclusion regions in the $m_{V}-g_{V}$ plane,
for the vector-mediator case. Different values of $m_{\chi}$ proportional to
$m_{V}$ are considered.

Figure 5: $90\%$ C.L. (2 d.o.f.) exclusion regions in the $m_{\chi}-g_{V}$
plane, for the vector-mediator case. Different values of $m_{\chi}$
proportional to $m_{V}$ are considered.

We start our analysis by fixing $m_{\chi}$ and letting the mediator mass and
the coupling vary. In Fig. 2 we choose several fixed values of $m_{\chi}$
($m_{\chi}=0.1,10$ and $50$ MeV, shown with darker to lighter shades) and plot
the $90\%$ C.L. exclusion regions (assuming 2 d.o.f.) for the CsI (top-left)
and LAr (top-right) COHERENT data. We also perform a combined CsI+LAr data
analysis (bottom-left) and we finally compare the three results (CsI, LAr and
the combined CsI+LAr) for the most constraining case, corresponding to
$m_{\chi}=0.1~{}\mathrm{MeV}$ (bottom-right). The most notable difference
among the two data sets is that CsI data lead to a significantly improved
sensitivity compared to the LAr data in the low mass region, due to the
inclusion of the ES events in the analysis. However, only the
$m_{\chi}=0.1~{}\mathrm{MeV}$ case seems to be sensitive to this effect. This
is due to the fact that the mass of the produced DF is bounded from above.
Indeed, as explained in Sec. II, there is an upper limit on $m_{\chi}$
dictated by the kinematics of the up-scattering process. In particular, for
CE$\nu$NS on Cs or I nuclei, $m_{\chi}\lesssim 53~{}\mathrm{MeV}$, while for
ES the bound is lower, $m_{\chi}\lesssim 7~{}\mathrm{MeV}$.

In the intermediate region, LAr is more sensitive than CsI, whereas in the
high mass region, CsI data lead to more stringent constraints. In the limit of
heavy mediator masses 555This limit corresponds to the case of effective
couplings $\epsilon^{V}\equiv\sqrt{2}g_{V}^{2}/\left(G_{F}m_{V}^{2}\right)$
explored in Ref. [73]., i.e. $m_{V}\gg|\mathbf{q}|$, or equivalently
$m_{V}\gtrsim 100~{}\mathrm{MeV}$, the data start to loose sensitivity on
$g_{V}$, since the cross section gets suppressed due to the presence of the
mediator mass in the denominator. On the contrary, when $m_{V}\ll|\mathbf{q}|$
the cross section becomes independent of $m_{V}$ and the exclusion contour
reaches a plateau. We further find that the $m_{\chi}=50~{}\mathrm{MeV}$
benchmark leads to a less stringent constraint on $g_{V}$ in the low-mass
region, compared to the $m_{\chi}=10~{}\mathrm{MeV}$ one. This is due to the
fact that the term proportional to $m_{\chi}$ in the kinetic factor of the
cross section (see Eq. (4)) is subtracting the term on the left. Thus, a
larger coupling $g_{V}$ is required in order to produce a sizeable number of
events, comparable to the experimental one. On the contrary, if $m_{\chi}$
decreases, then the coupling can also decrease and the corresponding
constraint is stronger. Finally, in the limit $m_{\chi}\to 0$, the light
mediator case from [33] is recovered, barring a small factor due to the
absence of interference terms between the SM and the new physics interactions
in the DF scenario. We have furthermore checked that the limits saturate for
$m_{\chi}\lesssim 0.1$ MeV, at low $m_{V}$.

Next, we perform the same analysis, but fixing the mediator mass, $m_{V}$.
Figure 3 shows the corresponding result at 90% C.L. in the $m_{\chi}-g_{V}$
plane for $m_{V}=0.1,\,10,\,100$ and $500$ MeV, depicted with darker to
lighter shades. In the CsI (top-left) panel, the impact of the ES events below
the bound of $m_{\chi}\sim 7~{}\mathrm{MeV}$ is clearly visible. Notice
however that the effect of the ES contribution only appears for
$m_{V}=0.1~{}\mathrm{MeV}$ 666The case of $m_{V}=1.0~{}\mathrm{MeV}$ was also
computed, but it overlaps the $m_{V}=0.1~{}\mathrm{MeV}$ contour and it is not
shown in the plot.. For higher $m_{V}$ values the ES cross section given in
Eq. (10) gets suppressed thus making its contribution to the total number of
events negligible. Let us emphasize that, without ES, LAr data (top-right)
provide better sensitivity on the coupling than CsI. For this reason, the
combination of both data sets (bottom row) produces a stronger constraint.
Notice also the kinematic limit which is again clearly visible in the figure,
$m_{\chi}\sim 53~{}\mathrm{MeV}$ for every possible value of $m_{V}$.

In the last set of analyses we set $m_{\chi}$ proportional to $m_{V}$ and we
let both the coupling and the mediator mass vary freely. In Fig. 4 we plot the
$90\%$ C.L. bounds on $g_{V}$ versus the mediator mass, for
$m_{\chi}=0.1,\,0.5,\,1,\,5,\,10~{}m_{V}$, from dark to light contours. The
main difference with respect to the previous figures is that the different
contours are not shifted vertically but horizontally. Take, for example, the
case with $m_{\chi}=10~{}m_{V}$. When $m_{V}=5~{}\mathrm{MeV}$, $m_{\chi}$ is
close to its upper bound and for values higher than that, the up-scattering
production of the DF is kinematically forbidden, hence only the SM CE$\nu$NS
remains. If we focus now on the contour for $m_{\chi}=0.5~{}\mathrm{MeV}$, the
upper bound on $m_{\chi}$ is achieved at higher values of $m_{V}$, hence it is
more constraining.

Finally, Fig. 5 depicts the same data as Fig. 4, except that it is plotted
with respect to the DF mass. The upper bound of $m_{\chi}\sim
53~{}\mathrm{MeV}$ is the same for every contour line. In the CsI panel (top-
left) the ES contribution is shifted to the left when $m_{\chi}$ gets smaller
with respect to $m_{V}$. The reason is that when the ES contribution is
allowed ($m_{\chi}\lesssim 7~{}\mathrm{MeV}$), the mediator mass is large and
the cross section gets suppressed.

###  IV.2 Scalar mediator

Similarly to the vector-mediator discussion, in this subsection we refer to
the Lagrangian given in Eq. (II) and to the cross sections given in Eqs. (6)
and (11). Thanks to the redefinition of the coupling
$g_{S}\equiv\sqrt{g_{\chi_{L}}g_{f}}$, with $f=u,\,d$ for CE$\nu$NS and
$f=e^{-}$ for ES, we can express our results in terms of the relevant
parameters, namely the DF mass ($m_{\chi}$), the mediator mass ($m_{S}$) and
one coupling ($g_{S}$). In analogy to the vector-mediator case described
above, this simplification is based on the assumption that the scalar mediator
has universal couplings to both quarks and electrons.

In Figs. 6, 7, 8 and 9 we present the $90\%$ C.L. (2 d.o.f.) exclusion regions
for the same benchmarks described in the previous subsection. Because of the
similarities in the general interpretation of the two results, we refer the
reader to the vector mediator discussion in subsection IV.1, while here we
will only comment about the main differences between scalar and vector
scenarios.

The first relevant variation is that the ES contribution — affecting the CsI
data analysis — is negligible in the scalar-mediator case. Notice, for
instance, that in the corresponding constraints obtained in Figs. 2 and 3, the
ES contribution led to a “bump” in the exclusion regions at $m_{\chi}\lesssim
7$ MeV, that is now absent. This is due to the different expressions for the
ES cross sections given in Eqs. (10) and (11). The scalar cross section is
suppressed by a factor $m_{e}E_{\mathrm{er}}/E_{\nu}^{2}$ compared to the
leading terms of the vector case. Indeed, in the ES cross section for the
vector mediator there is a constant term that is always present regardless of
the values of the masses or the energies appearing in the kinematic terms.
Without the ES contribution, the inferred bounds on the coupling $g_{S}$
saturate already at $m_{\chi}\lesssim 10~{}\mathrm{MeV}$ in Fig. 6 and at
$m_{S}\lesssim 10~{}\mathrm{MeV}$ in Fig. 7.

At last, let us comment about how our results compare to those of Ref. [67].
The bounds obtained in that reference rely on some assumptions, for instance
$F_{W}^{2}=1$. Taking into account differences in the definition of the
relevant coupling ($\bar{y}_{s}$ vs $g_{S}$), we notice that our results are
in general comparable: while we include a more recent data set with more
statistics, we also perform a more sophisticated statistical analysis.

Figure 6: $90\%$ C.L. (2 d.o.f.) exclusion regions for the scalar-mediator
case, in the $m_{S}-g_{S}$ plane. Different fixed values of $m_{\chi}$ are
considered.

Figure 7: $90\%$ C.L. (2 d.o.f.) exclusion regions for the scalar mediator
case, in the $m_{\chi}-g_{S}$ plane. Different fixed values of $m_{S}$ are
considered.

Figure 8: $90\%$ C.L. (2 d.o.f.) exclusion regions for the scalar mediator
case, in the $m_{S}-g_{S}$ plane. Different values of $m_{\chi}$ proportional
to $m_{S}$ are considered.

Figure 9: $90\%$ C.L. (2 d.o.f.) exclusion regions for the scalar mediator
case, in the $m_{\chi}-g_{S}$ plane. Different values of $m_{\chi}$
proportional to $m_{S}$ are considered. (With respect to Fig. 8 only the
x-axis is changed.)

###  IV.3 Comparison to other constraints

Figure 10: $90\%$ C.L. (2 d.o.f.) exclusion regions for the vector- (left) and
scalar-mediator (right) case. The orange contours are obtained from the
combined analysis of CsI+LAr data, and assuming a massless DF. We refer to the
main text for more details regarding this approximation.

In this subsection we discuss the limits obtained in this work in comparison
with other existing bounds in the literature. In principle, limits from other
experiments should be carefully recast in terms of the up-scattering scenario
that we considered. However we recall here that our exclusion contours
saturate for $m_{\chi}\lesssim 0.1$ MeV, and eventually coincide with the
limits expected for $m_{\chi}\to 0$. Given these considerations, we can
compare to existing constraints on CE$\nu$NS or neutrino ES with light vector
mediators. In Fig. 10 we hence show, as for comparison, our result from the
combined analysis together with limits from other experimental searches on
light mediators (vector on the left and scalar on the right). We display
limits from other CE$\nu$NS experiments, in particular CONNIE [108], CONUS
[15] and Dresden-II [32], from collider experiments [55], like LHCb [109] and
BaBar [110], from rare meson decays at NA48 [111], from Borexino solar
neutrino data [112] and from the analyses [113] of multi-ton DM experiments
(XENONnT [114] and LZ [115]). Let us finally comment that limits from deep
inelastic scattering (CHARM-II), from reactor neutrinos (TEXONO) and from
solar neutrinos (BOREXINO) have been set on effective interactions leading to
the up-scattering production of a DF [73]. While a direct comparison with
these results is not possible as we are considering light mediators, we can
notice that COHERENT can explore DF masses above the reach of TEXONO and
BOREXINO, and a bit below CHARM-II. To conclude, we see that the limits here
obtained using COHERENT data improve upon existing ones for mediators with
masses in the range $70~{}\mathrm{MeV}\lesssim m_{V}\lesssim
110~{}\mathrm{MeV}$ and $m_{S}\gtrsim 40$ MeV.

###  IV.4 The dark fermion as the dark matter

The DF produced through neutrino up-scattering at COHERENT, belonging to a
dark sector, could also constitute part or the totality of the DM observed in
the Universe. The evaluation of the cosmic abundance of $\chi$ will depend on
all possible production mechanisms. In this respect, in addition to those
induced by the interactions in Eq. (II), new ones could arise in UV-
completions of our phenomenological model. While a rigorous discussion about
the viability of $\chi$ as the DM would require a dedicated analysis and it is
therefore beyond the scope of this work, in the following, we will make few
considerations about its stability, focusing on the vector-mediator scenario.
Similar considerations can also be made for the scalar-mediator case. We refer
the reader to Refs. [67, 116, 73, 91] for further discussions.

First of all we can note that, given the interactions allowed by the
Lagrangian in Eq. (II), there are only two possibilities for $\chi$ to decay,
at tree-level. If $m_{\chi}>m_{V}$, $\chi$ can decay into $\chi\to
V\nu_{\alpha}$. If instead $m_{\chi}<m_{V}$, the three-body decay
$\chi\to\nu_{\alpha}f\bar{f}$ is allowed, via an off-shell mediator. Given the
interactions in Eq. (II) and the $\chi$ mass range accessible at COHERENT,
eventually only the channel $\chi\to\nu_{\alpha}e^{+}e^{-}$ is kinematically
allowed. In principle, the tree-level decay
$\chi\to\nu_{\alpha}\nu_{\beta}\bar{\nu}_{\beta}$ could also be allowed,
depending on the nature of the interactions given in Eq. (II). However, if the
new mediator $V$ couples only to charged fermions (e.g. in a dark-photon
scenario, see for example Ref. [117]), the tree-level decay into three
neutrinos is forbidden. We assume this possibility and we hence estimate the
tree-level decay widths as following

$\displaystyle\Gamma(\chi\to V\nu_{\alpha})$
$\displaystyle=\dfrac{g_{\chi_{L}}^{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_{\alpha}e^{-}e^{+})$
$\displaystyle=\dfrac{g_{V}^{4}}{768\pi^{3}}\,m_{\chi}^{5}\,I(m_{\chi},m_{V},\mu_{2})\,,$
(25)

where $I(m_{\chi},m_{V},\mu_{2})$ reads

$\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{\\}}.$
(26)

Here $\mu_{2}\equiv m_{e}^{2}/m_{\chi}^{2}$ and
$\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$ is the Källén function.
Equation (25) is valid in both the heavy- and light-mediator regimes. Let us
comment that to compare with other results in the literature, for instance
[118, 117, 119], where the heavy-mediator approximation is assumed ($m_{V}\gg
m_{\chi}$), one may use
$2I(m_{\chi},m_{V},\mu_{2})\approx[I_{2}(0,\sqrt{\mu_{2}},\sqrt{\mu_{2}})+2I_{1}(0,\sqrt{\mu_{2}},\sqrt{\mu_{2}})]/m_{V}^{4}$.
The kinematical functions $I_{1}$ and $I_{2}$ are defined, e.g., in Ref.
[120]. From these expressions, we can notice that for $m_{\chi}<2m_{e}$ no
tree-level decays are kinematically allowed. Radiative decays such as
$\chi\to\nu\gamma$, $\chi\to\nu\gamma\gamma\gamma$ or $\chi\to\nu\nu\nu$ could
be, on the other hand, open. These decays have been computed assuming
effective operators in [93].
Based on the previous decay rates, we estimate the lifetime of $\chi$ and
compare it to the age of the Universe. We find that when assuming a very light
mediator ($m_{V}\sim 1$ keV) and values of couplings testable at COHERENT, the
tree-level decays given in Eq. (25) are efficient enough to make $\chi$ decay
before $t\sim 1$ sec, independently of its mass. In this case, $\chi$ could
not be the DM, but also would not affect the BBN as it would decay at earlier
epochs. When the mediator is heavy ($m_{V}\gtrsim 1$ MeV), $\chi$ would still
decay faster than the age of the Universe, in the region of parameters where
the decay channels in Eq. (25) are kinematically allowed. Finally, should the
mediator be heavy ($m_{V}\gtrsim 1$ MeV) and the DF light enough ($m_{\chi}<1$
MeV) then the tree-level decays would not be kinematically accessible.
Estimations of radiative decays [93, 92, 91] also indicate that they would not
be very efficient thus making $\chi$ a long-lived particle and possibly a
viable DM candidate. Of course in all cases, very small couplings would always
allow $\chi$ stability, but would definitely be outside the reach of COHERENT.
As already mentioned, in a UV-completed model the possible presence of
additional interactions beyond those considered in Eq. (II) could be relevant
for the production of $\chi$ as the DM in the early Universe. We notice that
small couplings of the order of those probed by COHERENT already indicate a
possible nonthermal production, e.g. via the freeze-in mechanism. We conclude
mentioning that, even if stable, light (sub-)MeV DM could affect the
predictions for the BBN, thus requiring a careful evaluation of all
cosmological constraints.

Finally, going back to the region of parameter space where the DF is unstable,
one can investigate whether it decays fast enough to give rise to novel
signatures inside the COHERENT detector. Based on the COHERENT sensitivity
reach demonstrated e.g. in Fig. 4, we find that $\chi$ can decay inside the
detector only when $m_{\chi}>m_{V}$ i.e., via the $\chi\to V\nu_{\alpha}$
mode. However, if $m_{V}$ is also larger than twice the electron mass, $V$ can
subsequently decay into a $e^{-}e^{+}$ pair via

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

and possibly lead to a detectable electronic recoil signal within the COHERENT
detector.

##  V Conclusions

In this paper we have investigated the possible production of a dark fermion
$\chi$ at the COHERENT experiment, through neutrinos up-scattering off both
the nuclei and the atomic electrons of the detector. We have performed a
detailed statistical analysis of two COHERENT data-sets: the most recent one,
obtained with the CsI detector (2021) and the LAr one (2020). With the idea of
going beyond up-scattering through the dipole operator previously studied in
the literature, we have focused on two possible mediators: a light vector and
a light scalar. We have provided the relevant cross sections for the processes
$\nu\mathcal{N}\to\chi\mathcal{N}$ and $\nu e\to\chi e$, for both vector and
scalar mediators.

We have obtained $90\%$ C.L. exclusion regions on the relevant coupling and
mediator/dark fermion mass parameter space. Our results show that including
the ES in the CsI analysis improves the bounds on the couplings only in the
vector-mediator case. In the scalar-mediator scenario the ES events turn out
not to be relevant, nevertheless the most recent CsI data still allow to
improve the sensitivity on the couplings, compared to LAr data. Furthermore,
in comparison with other experiments and under the assumption of a very light
dark fermion, our combined analysis improves the existing bounds for mediators
with masses in the range $70~{}\mathrm{MeV}\lesssim m_{V}\lesssim
110~{}\mathrm{MeV}$ and $m_{S}\gtrsim 40$. Future data from COHERENT will
allow to further probe the possible up-scattering production of a dark
fermion.

We have finally briefly discussed the possibility for $\chi$ to comprise the
dark matter of the Universe. We have found that, when assuming a very light
mediator, given the masses and couplings accessible at COHERENT, the dark
fermion would not be stable enough to be the dark matter. In other regions of
the parameter space, the dark fermion could be stable however a careful
discussion of its role as the cosmological dark matter would require a
dedicated analysis.

## Acknowledgments

We thank Sergio Pastor for reading this manuscript and for discussions. This
work has been partially supported by the Spanish grants PID2020-113775GB-I00
(MCIN/AEI/10.13039/ 501100011033) and CIPROM/2021/054 (Generalitat
Valenciana). VDR acknowledges financial support by the CIDEGX/2022/20 grant
(project “D’AMAGAT”) funded by Generalitat Valenciana. PMC is supported by the
grant CIACIF/2021/281 also funded by Generalitat Valenciana. DKP 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 Details of the cross section calculation

In this appendix we provide the necessary steps for obtaining the up-
scattering cross section formulae given in Eqs. (4) and (6). To this purpose,
starting from the quark-level Lagrangians in Eq. (II) we need to follow two
steps: first we translate the quark-level cross sections to the nucleon level
and then from nucleons to the nucleus (see below for details).

###  A.1 Vector mediator

For the process $\nu_{\alpha}(p_{1})f(p_{2})\to\chi(p_{3})f(p_{4})$, the
corresponding tree-level amplitude at the quark level is

$\displaystyle i\mathcal{M}_{V}=$
$\displaystyle\,\dfrac{i}{t-m_{V}^{2}}\,\left[\overline{u}_{\chi}(p_{3},s_{3})\gamma_{\mu}g_{\chi_{L}}P_{L}u_{\nu}(p_{1},s_{1})\right]\left\\{\overline{u}_{f}(p_{4},s_{4})\gamma^{\mu}\left(g_{f_{L}}P_{L}+g_{f_{R}}P_{R}\right)u_{f}(p_{2},s_{2})\vphantom{\dfrac{m_{\chi}m_{f}}{m_{V}^{2}}}\right.$
$\displaystyle+\left.\dfrac{m_{\chi}m_{f}}{m_{V}^{2}}\,\overline{u}_{f}(p_{4},s_{4})\left[\left(g_{f_{L}}-g_{f_{R}}\right)P_{L}+\left(g_{f_{R}}-g_{f_{L}}\right)P_{R}\right]u_{f}(p_{2},s_{2})\right\\},$
(28)

where the Dirac equation $\not pu(p,s)=mu(p,s)$ has been used together with
the fact that only left-handed neutrinos and right-handed antineutrinos exist
within the SM, thus $P_{R}\,u_{\nu}(p_{1},s_{1})=0$. The amplitude has been
rewritten in this way to highlight the explicit cancellation of axial terms,
under the assumption of chirality-blind couplings. Indeed, in order to
simplify the calculations, in what follows the left- and right-handed
components of the charged-fermion couplings are assumed to be equal, i.e.,
$g_{f_{L}}=g_{f_{R}}\equiv g_{f}$. In addition, quark universal couplings are
assumed, $g_{u}=g_{d}=g_{f}$. Finally we also work under the approximation of
massless neutrinos, i.e., $m_{\nu}\to 0$. With this in mind, the amplitude
given in Eq. (A.1) is simplified considerably and eventually reads

$i\mathcal{M}_{V}=\,\dfrac{i}{t-m_{V}^{2}}\,g_{\chi_{L}}g_{f}\left[\overline{u}_{\chi}(p_{3},s_{3})\gamma_{\mu}P_{L}u_{\nu}(p_{1},s_{1})\right]\left[\overline{u}_{f}(p_{4},s_{4})\gamma^{\mu}u_{f}(p_{2},s_{2})\right]\,.$
(29)

Notice that only the vector component survives, while the axial and
pseudoscalar components vanish.

Next, we go from quarks to nucleons. This is achieved by computing the quark-
operator matrix elements between nucleon states following the procedure
usually adopted for DM direct detection searches and described in [99, 98].
For a vector-type interaction, we make the following change

$g_{f}\overline{u}_{f}\gamma^{\mu}u_{f}\rightarrow\sum_{f=u,\,d}g_{f}\left(c_{f}^{(p)}\overline{u}_{p}\gamma^{\mu}u_{p}+c_{f}^{(n)}\overline{u}_{n}\gamma^{\mu}u_{n}\right),$
(30)

where the coefficients for the proton and neutron are
$c_{u}^{(p)}=c_{d}^{(n)}=2$, $c_{d}^{(p)}=c_{u}^{(n)}=1$.

In order to compute the scattering amplitude for the whole nucleus, an
additional step is required. Let $Z$ be the number of protons inside the
nucleus and $N=A-Z$ the number of neutrons, with $A$ being the mass number,
then one should proceed by changing

$g_{f}\left(\overline{u}_{p}\gamma^{\mu}u_{p}+\overline{u}_{n}\gamma^{\mu}u_{n}\right)\rightarrow
g_{f}\left(ZF_{W_{p}}+NF_{W_{n}}\right)\overline{u}_{\mathcal{N}}\gamma^{\mu}u_{\mathcal{N}},$
(31)

where $F_{W_{p}}$ and $F_{W_{n}}$ are the nuclear form factors for protons and
neutrons, respectively. Finally, we assume that both form factors are equal
777This is an accurate assumption for the typical momentum transfer involved
at COHERENT., $F_{W_{p}}=F_{W_{n}}\equiv F_{W}$, we apply the transformations
described above and perform the quark summation to obtain

$i\mathcal{M}_{V}=\,\dfrac{i}{t-m_{V}^{2}}\,3g_{\chi_{L}}g_{f}F_{W}(\left|\mathbf{q}\right|^{2})A\left[\overline{u}_{\chi}(p_{3},s_{3})\gamma_{\mu}P_{L}u_{\nu}(p_{1},s_{1})\right]\left[\overline{u}_{\mathcal{N}}(p_{4},s_{4})\gamma^{\mu}u_{\mathcal{N}}(p_{2},s_{2})\right].$
(32)

We define the vector coupling $g_{V}\equiv\sqrt{g_{\chi_{L}}g_{f}}$ to reduce
the number of parameters by one. At this point, one has just to perform the
spin summation, rewrite the four-momenta in the lab frame and compute the
cross section as a function of the nuclear recoil energy to obtain Eq. (4).

Once we have computed the CE$\nu$NS cross section, obtaining the ES one is
fairly simple. Since quarks and electrons are fermions, their quantum fields
have the same structure, hence we can start with Eq. (29) and take $f=e^{-}$.
Quark universality of $g_{f}$ is extended to fermion universality, which now
includes electrons, i.e., $g_{f}\equiv g_{e}=g_{u}=g_{d}$. The cross section
obtained with this amplitude corresponds to a neutrino up-scattering off a
single and isolated electron. However, we need to take into account that
electrons are bound inside the atomic nucleus. To obtain Eq. (10), it is
required to weigh the free ES cross section with the effective charge,
$Z_{\text{eff}}^{\mathcal{A}}$, corresponding to an energy deposition
$E_{\text{er}}$. The effective charges for the Cs and I isotopes are given in
Appendix B.

###  A.2 Scalar mediator

The procedure followed to obtain the DF-production cross section through a
scalar interaction is very similar to the one described above for the case of
a vector mediator. For the process
$\nu_{\alpha}(p_{1})f(p_{2})\to\chi(p_{3})f(p_{4})$, the corresponding tree-
level amplitude is

$i\mathcal{M}_{S}=\,\dfrac{i}{t-m_{S}^{2}}\left[\overline{u}_{\chi}(p_{3},s_{3})g_{\chi_{L}}P_{L}u_{\nu}(p_{1},s_{1})\right]\left[\overline{u}_{f}(p_{4},s_{4})\left(g_{f_{L}}P_{L}+g_{f_{R}}P_{R}\right)u_{f}(p_{2},s_{2})\right].$
(33)

Like before, we assume that the left- and right-handed components of the
couplings are the same, and also that fermion universality is satisfied. Next,
we go from quarks to nucleons by computing the relevant quark current in
nucleons. For a scalar-type interaction, it proves convenient to perform the
following modification [99, 98]

$g_{f}\overline{u}_{f}u_{f}\rightarrow\sum_{f=u,\,d}g_{f}\left(\dfrac{m_{p}}{m_{f}}f_{T_{f}}^{(p)}\,\overline{u}_{p}u_{p}+\dfrac{m_{n}}{m_{f}}f_{T_{f}}^{(n)}\,\overline{u}_{n}u_{n}\right),$
(34)

where $m_{p}$ is the proton mass, $m_{n}$ is the neutron mass, $m_{f}$ is the
quark mass and $f_{T_{f}}^{(p)}$ and $f_{T_{f}}^{(n)}$ express the quark-mass
contributions to the nucleon mass. The chosen values are given in Eq. (8).

In the last step, we translate from the nucleon to the nuclear level in the
same way as done for the vector interaction (see Eq. (31)). After performing
the quark summation, we obtain

$i\mathcal{M}_{S}=\,\dfrac{i}{t-m_{S}^{2}}C_{S}^{2}F_{W}(\left|\mathbf{q}\right|^{2})\left[\overline{u}_{\chi}(p_{3},s_{3})P_{L}u_{\nu}(p_{1},s_{1})\right]\left[\overline{u}_{\mathcal{N}}(p_{4},s_{4})u_{\mathcal{N}}(p_{2},s_{2})\right],$
(35)

where $C_{S}^{2}$ is the scalar coupling found in Eq. (7). We have defined
$g_{S}\equiv\sqrt{g_{\chi_{L}}g_{f}}$ to reduce the number of free parameters.
Given this amplitude one can compute the differential cross section in Eq.
(6). To obtain the ES one, we start from the amplitude in Eq. (33) and repeat
the procedure explained previously for the case of a vector mediator. The
result is given in Eq. (11).

##  Appendix B Effective electron charge for Cs and I

In this Appendix we provide the effective electron charge for Cs and I (see
table 1).

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

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

Table 1: Effective electron charge for Cs and I as a function of the energy
deposition $E_{\mathrm{er}}$ [121].

##  Appendix C Detector-specific quantities for CsI and LAr

We convert the true nuclear recoil energy into electron equivalent energy
$E_{\mathrm{er}}$ through the QF. For the CsI detector this is done in terms
of the light yield $\text{LY}=13.35~{}\mathrm{PE/keV_{ee}}$, with
PE$=\mathrm{LY}\times E_{\mathrm{er}}$, where PE denotes the number of
photoelectrons while the electron-equivalent energy $E_{\mathrm{er}}$ is
obtained from the scintillation curve, as
$E_{\mathrm{er}}=x_{1}E_{\mathrm{nr}}^{\prime}+x_{2}E_{\mathrm{nr}}^{\prime
2}+x_{3}E_{\mathrm{nr}}^{\prime 3}+x_{4}E_{\mathrm{nr}}^{\prime 4}$
($x_{1}=0.0554628,~{}x_{2}=4.30681,~{}x_{3}=-111.707,~{}x_{4}=840.384$) [10].
For the LAr detector the QF is given by
$\mathrm{QF}(E_{\mathrm{nr}}^{\prime})=0.246+7.8\times
10^{-4}E_{\mathrm{nr}}^{\prime}(\mathrm{keV_{nr}})$ [11].

The energy-dependent efficiency for the CsI detector is given by

$\epsilon_{E}(x)=\frac{a}{1+e^{-b(x-c)}}+d\,,$ (36)

where $x=\mathrm{PE}+\alpha_{7}$ and $a=1.32045$, $b=0.285979$, $c=10.8646$,
$d=-0.333322$ [10]. Following Ref. [33] we account for the $1\sigma$
uncertainty on the efficiency curve through the parameter $\alpha_{7}$, see
Eq. (21), which is allowed to float freely between $[-1,+1]\times$PE. For the
LAr detector there is no analytical efficiency function and we utilize the
data provided in Ref. [104] without introducing any nuisance parameter.

The time-dependent efficiency at the CsI detector is given by [10]

$\epsilon_{T}(t_{\mathrm{rec}})=\begin{cases}1,&\text{if}\
t_{\mathrm{rec}}<t^{\prime}\\\
e^{-k(t_{\mathrm{rec}}-t^{\prime})},&\text{if}~{}t_{\mathrm{rec}}\geq
t^{\prime}\end{cases}$ (37)

with $t^{\prime}=0.52~{}\mathrm{\mu s}$ and $k=0.0494\mathrm{/\mu s}$. In this
case we allow for a $\pm 250$ ns variation of the timing distribution by
introducing the nuisance parameter $\alpha_{6}$, with no prior assigned. For
the LAr detector it holds $\epsilon_{T}(t_{\mathrm{rec}})=1$.

Finally, the resolution function
$\mathcal{R}(E_{\text{nr}},E^{\prime}_{\text{nr}})$ for the CsI detector is
given by

$\mathcal{R}(E_{\text{nr}},E^{\prime}_{\text{nr}})=\frac{(\mathfrak{a}(1+\mathfrak{b}))^{1+\mathfrak{b}}}{\Gamma(1+\mathfrak{b})}\cdot
x^{\mathfrak{b}}\cdot e^{-\mathfrak{a}(1+\mathfrak{b})x},$ (38)

where $x$ is the reconstructed recoil energy expressed in PE units i.e.
PE$(E_{\mathrm{nr}})$, while $\mathfrak{a}$ and $\mathfrak{b}$ instead depend
on the true quenched energy deposition:
$\mathfrak{a}=0.0749/E_{\mathrm{er}}(E^{\prime}_{\text{nr}})$,
$\mathfrak{b}=9.56\times E_{\mathrm{er}}(E^{\prime}_{\text{nr}})$ [10]. For
the LAr detector instead, the resolution function is approximated by a
normalized Gaussian with resolution power [104]

$\dfrac{\sigma_{E_{\mathrm{er}}}}{E_{\mathrm{er}}}=\dfrac{0.58}{\sqrt{E_{\mathrm{er}}(\mathrm{keV_{ee}})}}.$
(39)

## References

  * [1] T. Kajita Rev. Mod. Phys. 88 (2016) 030501.
  * [2] A. B. McDonald Rev. Mod. Phys. 88 (2016) 030502.
  * [3] 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].
  * [4] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou JHEP 09 (2020) 178, arXiv:2007.14792 [hep-ph].
  * [5] M. Sajjad Athar, A. Fatima, and S. K. Singh Prog. Part. Nucl. Phys. 129 (2023) 104019, arXiv:2206.13792 [hep-ph].
  * [6] D. Z. Freedman Phys. Rev. D 9 (1974) 1389–1392.
  * [7] M. Abdullah et al. arXiv:2203.07361 [hep-ph].
  * [8] D. K. Papoulias, T. S. Kosmas, and Y. Kuno Front. in Phys. 7 (2019) 191, arXiv:1911.00916 [hep-ph].
  * [9] COHERENT Collaboration, D. Akimov et al. Science 357 (2017) 1123–1126, arXiv:1708.01294 [nucl-ex].
  * [10] COHERENT Collaboration, D. Akimov et al. Phys. Rev. Lett. 129 (2022) 081801, arXiv:2110.07730 [hep-ex].
  * [11] COHERENT Collaboration, D. Akimov et al. Phys. Rev. Lett. 126 (2021) 012002, arXiv:2003.10630 [nucl-ex].
  * [12] CCM Collaboration, A. A. Aguilar-Arevalo et al. Phys. Rev. D 106 (2022) 012001, arXiv:2105.14020 [hep-ex].
  * [13] D. Baxter et al. JHEP 02 (2020) 123, arXiv:1911.00762 [physics.ins-det].
  * [14] CONNIE Collaboration, A. Aguilar-Arevalo et al. JHEP 05 (2022) 017, arXiv:2110.13033 [hep-ex].
  * [15] CONUS Collaboration, H. Bonet et al. JHEP 05 (2022) 085, arXiv:2110.02174 [hep-ph].
  * [16] $\nu$GeN Collaboration, I. Alekseev et al. Phys.Rev.D 106 (2022) L051101, arXiv:2205.04305 [nucl-ex].
  * [17] MINER Collaboration, G. Agnolet et al. Nucl. Instrum. Meth. A 853 (2017) 53–60, arXiv:1609.02066 [physics.ins-det].
  * [18] J. Billard et al. J. Phys. G 44 (2017) 105101, arXiv:1612.09035 [physics.ins-det].
  * [19] R. Strauss et al. Eur. Phys. J. C 77 (2017) 506, arXiv:1704.04320 [physics.ins-det].
  * [20] H. T.-K. Wong The Universe 3 (2015) 22–37, arXiv:1608.00306 [hep-ex].
  * [21] G. Fernandez-Moroni, P. A. N. Machado, I. Martinez-Soler, Y. F. Perez-Gonzalez, D. Rodrigues, and S. Rosauro-Alcaraz JHEP 03 (2021) 186, arXiv:2009.10741 [hep-ph].
  * [22] D. Y. Akimov et al. JINST 17 (2022) T11011, arXiv:2209.15516 [physics.ins-det].
  * [23] NEON Collaboration, J. J. Choi et al. Eur. Phys. J. C 83 (2023) 226, arXiv:2204.06318 [hep-ex].
  * [24] NEWS-G Collaboration, L. Balogh et al. Phys. Rev. D 105 (2022) 052004, arXiv:2109.01055 [physics.ins-det].
  * [25] SBC, CE$\nu$NS Theory Group at IF-UNAM Collaboration, L. J. Flores et al. Phys. Rev. D 103 (2021) L091301, arXiv:2101.08785 [hep-ex].
  * [26] J. Colaresi, J. I. Collar, T. W. Hossbach, C. M. Lewis, and K. M. Yocum Phys. Rev. Lett. 129 (2022) 211802, arXiv:2202.09672 [hep-ex].
  * [27] $\nu$BDX-DRIFT Collaboration, D. Aristizabal Sierra, J. L. Barrow, B. Dutta, D. Kim, D. Snowden-Ifft, L. Strigari, and M. H. Wood arXiv:2210.08612 [hep-ex].
  * [28] D. K. Papoulias, T. S. Kosmas, R. Sahu, V. K. B. Kota, and M. Hota Phys. Lett. B 800 (2020) 135133, arXiv:1903.03722 [hep-ph].
  * [29] P. Coloma, I. Esteban, M. C. Gonzalez-Garcia, and J. Menendez JHEP 08 (2020) 030, arXiv:2006.08624 [hep-ph].
  * [30] M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, C. A. Ternes, and Y. Y. Zhang Phys. Rev. C 104 (2021) 065502, arXiv:2102.06153 [hep-ph].
  * [31] A. Majumdar, D. K. Papoulias, R. Srivastava, and J. W. F. Valle arXiv:2208.13262 [hep-ph].
  * [32] D. Aristizabal Sierra, V. De Romeri, and D. K. Papoulias JHEP 09 (2022) 076, arXiv:2203.02414 [hep-ph].
  * [33] 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].
  * [34] D. A. Sierra arXiv:2301.13249 [hep-ph].
  * [35] M. Atzori Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, and G. Masia arXiv:2303.09360 [nucl-ex].
  * [36] J. Barranco, O. G. Miranda, and T. I. Rashba JHEP 12 (2005) 021, arXiv:hep-ph/0508299.
  * [37] T. Ohlsson Rept. Prog. Phys. 76 (2013) 044201, arXiv:1209.2710 [hep-ph].
  * [38] O. G. Miranda and H. Nunokawa New J. Phys. 17 (2015) 095002, arXiv:1505.06254 [hep-ph].
  * [39] J. B. Dent, B. Dutta, S. Liao, J. L. Newstead, L. E. Strigari, and J. W. Walker Phys. Rev. D 96 (2017) 095007, arXiv:1612.06350 [hep-ph].
  * [40] Y. Farzan and M. Tortola Front. in Phys. 6 (2018) 10, arXiv:1710.09360 [hep-ph].
  * [41] J. Liao and D. Marfatia Phys. Lett. B 775 (2017) 54–57, arXiv:1708.04255 [hep-ph].
  * [42] D. Aristizabal Sierra, V. De Romeri, and N. Rojas Phys. Rev. D 98 (2018) 075018, arXiv:1806.07424 [hep-ph].
  * [43] M. Abdullah, J. B. Dent, B. Dutta, G. L. Kane, S. Liao, and L. E. Strigari Phys. Rev. D 98 (2018) 015005, arXiv:1803.01224 [hep-ph].
  * [44] J. Billard, J. Johnston, and B. J. Kavanagh JCAP 11 (2018) 016, arXiv:1805.01798 [hep-ph].
  * [45] P. B. Denton, Y. Farzan, and I. M. Shoemaker JHEP 07 (2018) 037, arXiv:1804.03660 [hep-ph].
  * [46] T. Han, J. Liao, H. Liu, and D. Marfatia JHEP 11 (2019) 028, arXiv:1910.03272 [hep-ph].
  * [47] C. Giunti Phys. Rev. D 101 (2020) 035039, arXiv:1909.00466 [hep-ph].
  * [48] D. Aristizabal Sierra, B. Dutta, S. Liao, and L. E. Strigari JHEP 12 (2019) 124, arXiv:1910.12437 [hep-ph].
  * [49] D. Aristizabal Sierra, V. De Romeri, and N. Rojas JHEP 09 (2019) 069, arXiv:1906.01156 [hep-ph].
  * [50] P. B. Denton and J. Gehrlein JHEP 04 (2021) 266, arXiv:2008.06062 [hep-ph].
  * [51] L. J. Flores, N. Nath, and E. Peinado JHEP 06 (2020) 045, arXiv:2002.12342 [hep-ph].
  * [52] M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, and Y. Y. Zhang JHEP 01 (2021) 116, arXiv:2008.05022 [hep-ph].
  * [53] D. W. P. Amaral, D. G. Cerdeno, A. Cheek, and P. Foldenauer Eur. Phys. J. C 81 (2021) 861, arXiv:2104.03297 [hep-ph].
  * [54] L. J. Flores, N. Nath, and E. Peinado Phys. Rev. D 105 (2022) 055010, arXiv:2112.05103 [hep-ph].
  * [55] M. Atzori Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, E. Picciau, C. A. Ternes, and Y. Y. Zhang JHEP 05 (2022) 109, arXiv:2202.11002 [hep-ph].
  * [56] A. Elpe, E. Akyumuk, T. M. Aliev, L. Selbuz, and I. Turan arXiv:2212.06861 [hep-ph].
  * [57] J. Schechter and J. W. F. Valle Phys.Rev.D 24 (1981) 1883–1889.
  * [58] P. B. Pal and L. Wolfenstein Phys. Rev. D 25 (1982) 766.
  * [59] B. Kayser Phys. Rev. D 26 (1982) 1662.
  * [60] J. F. Nieves Phys. Rev. D 26 (1982) 3152.
  * [61] R. E. Shrock Nucl. Phys. B 206 (1982) 359–379.
  * [62] T. S. Kosmas, O. G. Miranda, D. K. Papoulias, M. Tortola, and J. W. F. Valle Phys. Rev. D92 (2015) 013011, arXiv:1505.03202 [hep-ph].
  * [63] B. C. Canas, O. G. Miranda, A. Parada, M. Tortola, and J. W. F. Valle Phys. Lett. B 753 (2016) 191–198, arXiv:1510.01684 [hep-ph]. [Addendum: Phys.Lett.B 757, 568–568 (2016)].
  * [64] O. G. Miranda, D. K. Papoulias, M. Tórtola, and J. W. F. Valle JHEP 07 (2019) 103, arXiv:1905.03750 [hep-ph].
  * [65] O. Miranda et al. Phys.Lett.B 808 (2020) 135685, arXiv:2007.01765 [hep-ph].
  * [66] T. S. Kosmas, D. K. Papoulias, M. Tortola, and J. W. F. Valle Phys. Rev. D 96 (2017) 063013, arXiv:1703.00054 [hep-ph].
  * [67] V. Brdar, W. Rodejohann, and X.-J. Xu JHEP 12 (2018) 024, arXiv:1810.03626 [hep-ph].
  * [68] C. Blanco, D. Hooper, and P. Machado Phys. Rev. D 101 (2020) 075051, arXiv:1901.08094 [hep-ph].
  * [69] 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].
  * [70] 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].
  * [71] P. D. Bolton, F. F. Deppisch, K. Fridell, J. Harz, C. Hati, and S. Kulkarni Phys. Rev. D 106 (2022) 035036, arXiv:2110.02233 [hep-ph].
  * [72] W. Chao, T. Li, J. Liao, and M. Su Phys. Rev. D 104 (2021) 095017, arXiv:2108.02341 [hep-ph].
  * [73] Z. Chen, T. Li, and J. Liao JHEP 05 (2021) 131, arXiv:2102.09784 [hep-ph].
  * [74] D. Aristizabal Sierra, O. G. Miranda, D. K. Papoulias, and G. S. Garcia Phys. Rev. D 105 (2022) 035027, arXiv:2112.12817 [hep-ph].
  * [75] R. Calabrese, J. Gunn, G. Miele, S. Morisi, S. Roy, and P. Santorelli Phys. Rev. D 107 (2023) 055039, arXiv:2212.11210 [hep-ph].
  * [76] D. McKeen and M. Pospelov Phys. Rev. D 82 (2010) 113018, arXiv:1011.3046 [hep-ph].
  * [77] P. deNiverville, M. Pospelov, and A. Ritz Phys. Rev. D 84 (2011) 075020, arXiv:1107.4580 [hep-ph].
  * [78] A. Berlin, P. deNiverville, A. Ritz, P. Schuster, and N. Toro Phys. Rev. D 102 (2020) 095011, arXiv:2003.03379 [hep-ph].
  * [79] B. Batell et al. in 2022 Snowmass Summer Study.  7, 2022.  arXiv:2207.06898 [hep-ph].
  * [80] S.-F. Ge and I. M. Shoemaker JHEP 11 (2018) 066, arXiv:1710.10889 [hep-ph].
  * [81] CCM Collaboration Collaboration, A. A. Aguilar-Arevalo et al. Phys. Rev. D 106 (Jul, 2022) 012001. https://link.aps.org/doi/10.1103/PhysRevD.106.012001.
  * [82] COHERENT Collaboration, D. Akimov et al. Phys. Rev. D 102 (2020) 052007, arXiv:1911.06422 [hep-ex].
  * [83] COHERENT Collaboration, D. Akimov et al. Phys. Rev. Lett. 130 (2023) 051803, arXiv:2110.11453 [hep-ex].
  * [84] CCM Collaboration, A. A. Aguilar-Arevalo et al. Phys. Rev. Lett. 129 (2022) 021801, arXiv:2109.14146 [hep-ex].
  * [85] B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, and L. E. Strigari Phys. Rev. Lett. 124 (2020) 121802, arXiv:1906.10745 [hep-ph].
  * [86] B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, L. E. Strigari, and A. Thompson JHEP 01 (2022) 144, arXiv:2006.09386 [hep-ph].
  * [87] P. deNiverville, M. Pospelov, and A. Ritz Phys. Rev. D 92 (2015) 095005, arXiv:1505.07805 [hep-ph].
  * [88] COHERENT Collaboration, D. Y. Akimov et al. Phys. Rev. D 106 (2022) 052004, arXiv:2205.12414 [hep-ex].
  * [89] B. Dutta, D. Kim, S. Liao, J.-C. Park, S. Shin, and L. E. Strigari Phys. Rev. Lett. 124 (Mar, 2020) 121802. https://link.aps.org/doi/10.1103/PhysRevLett.124.121802.
  * [90] B. Dutta, W.-C. Huang, and J. L. Newstead arXiv:2302.10250 [hep-ph].
  * [91] T. Li and J. Liao JHEP 02 (2021) 099, arXiv:2008.00743 [hep-ph].
  * [92] J. A. Dror, G. Elor, and R. McGehee Phys. Rev. Lett. 124 (May, 2020) 181301. https://link.aps.org/doi/10.1103/PhysRevLett.124.181301.
  * [93] J. A. Dror, G. Elor, R. McGehee, and T.-T. Yu Phys. Rev. D 103 (2021) 035001, arXiv:2011.01940 [hep-ph]. [Erratum: Phys.Rev.D 105, 119903 (2022)].
  * [94] Y. Farzan, M. Lindner, W. Rodejohann, and X.-J. Xu JHEP 05 (2018) 066, arXiv:1802.05171 [hep-ph].
  * [95] E. Bertuzzo, G. Grilli di Cortona, and L. M. D. Ramos JHEP 06 (2022) 075, arXiv:2112.04020 [hep-ph].
  * [96] M. Bauer, P. Foldenauer, and J. Jaeckel JHEP 07 (2018) 094, arXiv:1803.05466 [hep-ph].
  * [97] P. Ilten, Y. Soreq, M. Williams, and W. Xue JHEP 06 (2018) 004, arXiv:1801.04847 [hep-ph].
  * [98] M. Cirelli, E. Del Nobile, and P. Panci JCAP 10 (2013) 019, arXiv:1307.5955 [hep-ph].
  * [99] E. Del Nobile arXiv:2104.12785 [hep-ph].
  * [100] J. Ellis, N. Nagata, and K. A. Olive Eur. Phys. J. C 78 (2018) 569, arXiv:1805.09795 [hep-ph].
  * [101] 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)].
  * [102] A. Majumdar, D. K. Papoulias, and R. Srivastava Phys. Rev. D 106 (2022) 013001, arXiv:2112.03309 [hep-ph].
  * [103] C. Giunti and C. W. Kim.  Oxford University Press, 03, 2007.  https://doi.org/10.1093/acprof:oso/9780198508717.001.0001.
  * [104] COHERENT Collaboration, D. Akimov et al. arXiv:2006.12659 [nucl-ex].
  * [105] P. Coloma, I. Esteban, M. C. Gonzalez-Garcia, L. Larizgoitia, F. Monrabal, and S. Palomares-Ruiz JHEP 05 (2022) 037, arXiv:2202.10829 [hep-ph].
  * [106] M. Atzori Corona, M. Cadeddu, N. Cargioli, F. Dordei, C. Giunti, Y. F. Li, C. A. Ternes, and Y. Y. Zhang JHEP 09 (2022) 164, arXiv:2205.09484 [hep-ph].
  * [107] E. Picciau.  PhD thesis, Cagliari U., 2022.
  * [108] CONNIE Collaboration, A. Aguilar-Arevalo et al. JHEP 04 (2020) 054, arXiv:1910.04951 [hep-ex].
  * [109] LHCb Collaboration, R. Aaij et al. Phys. Rev. Lett. 120 (2018) 061801, arXiv:1710.02867 [hep-ex].
  * [110] BaBar Collaboration, J. P. Lees et al. Phys. Rev. Lett. 113 (2014) 201801, arXiv:1406.2980 [hep-ex].
  * [111] NA48/2 Collaboration, J. R. Batley et al. Phys. Lett. B 746 (2015) 178–185, arXiv:1504.00607 [hep-ex].
  * [112] P. Coloma, M. C. Gonzalez-Garcia, M. Maltoni, J. a. P. Pinheiro, and S. Urrea JHEP 07 (2022) 138, arXiv:2204.03011 [hep-ph].
  * [113] S. K. A., A. Majumdar, D. K. Papoulias, H. Prajapati, and R. Srivastava Phys. Lett. B 839 (2023) 137742, arXiv:2208.06415 [hep-ph].
  * [114] (XENON Collaboration)††, XENON Collaboration, E. Aprile et al. Phys. Rev. Lett. 129 (2022) 161805, arXiv:2207.11330 [hep-ex].
  * [115] LZ Collaboration, J. Aalbers et al. arXiv:2207.03764 [hep-ex].
  * [116] N. Hurtado, H. Mir, I. M. Shoemaker, E. Welch, and J. Wyenberg Phys. Rev. D 102 (2020) 015006, arXiv:2005.13384 [hep-ph].
  * [117] Y. Jho, J. Kim, P. Ko, and S. C. Park arXiv:2008.12598 [hep-ph].
  * [118] E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal Phys. Rev. Lett. 121 (2018) 241801, arXiv:1807.09877 [hep-ph].
  * [119] K. Jodłowski and S. Trojanowski JHEP 05 (2021) 191, arXiv:2011.04751 [hep-ph].
  * [120] J. C. Helo, S. Kovalenko, and I. Schmidt Nucl. Phys. B 853 (2011) 80–104, arXiv:1005.1607 [hep-ph].
  * [121] A. Thompson et al., 2009.  https://xdb.lbl.gov/.