Neutrino non-standard interactions and dark matter searches with multi-ton scale detectors
D. Aristizabal Sierra N. Rojas M. H. G. Tytgat
hep-ph









  
  
  
  
  
  
  

  


  Neutrino non-standard interactions 
    and dark matter
    

    searches with multi-ton scale detectors
  

  D. Aristizabal Sierra^a,b,[e-mail address: daristizabal@ulg.ac.be],
  N. Rojas^a,[e-mail address: nicolas.rojasro@usm.cl],
  M.H.G. Tytgat^c,[email address: mtytgat@ulb.ac.be]
  


  ^a

  ^b

  ^c





  
  
  Future dark matter (DM) direct detection searches will be subject to
  irreducible neutrino backgrounds that will challenge the
  identification of an actual WIMP signal in experiments without
  directionality sensitivity. We study the impact of neutrino-quark
  non-standard interactions (NSI) on this background, assuming the
  constraints from neutrino oscillations and the recent COHERENT
  experiment data, which are relevant for NSI mediated by light
  mediators, m_ med≲ O(GeV). We calculate the
  expected number of neutrino-nucleus elastic scattering events in a
  Xe-based ton-size dark matter detector, including solar neutrino
  fluxes from the pp chain and CNO cycle as well as sub-GeV
  atmospheric fluxes and taking into account NSI effects in both
  propagation and detection. We find that sizable deviations from the
  standard model expectation are possible, but are more pronounced for
  flavor-diagonal couplings, in particular for electron neutrinos. We
  show that neutrino NSI can enhance or deplete the neutrino-nucleus
  event rate, which may impact DM searches in multi-ton detectors.







§ INTRODUCTION


A wealth of astrophysical and cosmological data support the idea that
the dominant form of matter in the universe is made of dark matter
(DM).  Several DM candidates have been proposed ranging from
astrophysical observable objects <cit.> up to
fundamental particles <cit.>. Among the latter, the most
studied candidates (both theoretically and experimentally) are massive
particles whose abundance is determined by thermal freeze-out (weakly
interacting massive particles or WIMPs).

WIMPs are being searched for using direct detection experiments since
the late 80's, following a proposal by Goodman and Witten
<cit.>. These experiments rely on energy deposition
measurements from WIMP-induced nuclear recoil processes, which are
expected to result as the Earth (detector) moves in the galactic DM
halo. The non-observation of such signatures are translated into upper
limits for the WIMP-nucleon (spin-independent or spin-dependent) cross
section as a function of the DM mass. At present the most stringent
bound is given by the XENON1T collaboration, which in the absence of a
nuclear recoil signal has derived the limit
σ_χ-n≲ 7.7× 10^-47 cm^2/c^2 for
m_χ=35 GeV/c^2 at 90% CL <cit.>. Searches
beyond these values and down to, say,
σ_χ-n∼ 10^-49 cm^2/c^2 require multi-ton
target masses, as planned for the DARWIN detector
<cit.>.

The reach of these experiments, however, will be limited by an
irreducible neutrino background. Indeed, coherent neutrino nucleus
elastic scattering (CνNES) induced by solar and, eventually,
atmospheric neutrino fluxes will mimic a possible WIMP-nucleus signal
<cit.>. While the precise extent of this background is
affected by various factors, including neutrino flux uncertainties and
astrophysical parameters <cit.>, in the
SM, the CνNES cross section is well understood and is not a source
of uncertainty. This can potentially change in the presence of new
physics. The most obvious way through which this can occur is via
neutrino-quark non-standard interactions (NSI), which historically
were introduced by Wolfenstein as a possible way to account for the
solar neutrino deficit <cit.>.

Neutrino NSI are constrained by several observations, including
neutrino oscillation data <cit.>, neutrino
deep inelastic scattering experimental data (CHARM
<cit.> and NuTeV <cit.>) and,
potentially, the measurement of the CνNES process itself, recently
observed by the COHERENT experiment <cit.> (see
<cit.> for a review and <cit.> for the
detailed analysis). CHARM and NuTeV set strong constraints on
ν_μ-related NSI parameters, but their relevance for constraining
CνNES depends on the range of application of the NSI effective
operators. In particular, these constraints do not apply if the NSI
are mediated by a particle with mass m_ med <  O(1) GeV
<cit.>.  In what follows we will consider the
possibility that NSI are mediated by a particle in the mass range
10 MeV≲ m_ med≲ 1 GeV, range in
which CHARM and NuTeV bounds do not apply, but one is still allowed to
use the constraints on NSI as extracted from a combined analysis of
COHERENT and oscillation data in ref.<cit.>. It is
worth pointing out that in this mass range bounds are less severe for
all flavor-diagonal NSI couplings, and so compared to the case where
deep inelastic scattering data apply one expects larger deviations on
the neutrino-nucleus event rates.

Specifically, using current constraints on neutrino NSI derived from
neutrino oscillation and COHERENT data <cit.>, in this
paper we quantify the deviations on the CνNES expected number of
events in a Xe-based DM detector using an exposure
ℰ=1 (ton·year). For that aim we
consider solar neutrino fluxes for which we include neutrinos from the
pp chain (pp, hep, ^8B, ^7Be, pep) and CNO cycle
(^13N, ^15O, ^17F). We consider as well
sub-GeV atmospheric neutrino fluxes. For the solar sector analysis we
calculate NSI propagation effects in the mass dominance limit
(Δ m_31^2→∞) and for the calculation of the electron
neutrino survival probability we take into account the distribution of
neutrino production in the Sun according to the standard solar model
(SSM) BS05 <cit.>.  We quantify matter-induced effects
in the atmospheric sector and find they are negligible
(H_mat/H_vac≲ 0.1) <cit.>,
however we include neutrino flavor oscillations in vacuum in the
calculation of the number of neutrino-nucleus events.

The study of beyond-the-standard model physics in neutrino-quark
interactions in the context of multi-ton DM detectors, has been
recently the subject of various analyses. Ref. <cit.>
studied the capability of such detectors to constraint light mediators
scenarios, while ref. <cit.> instead considered
scenarios where the new physics couples both neutrinos and
DM. Ref. <cit.> considered the implications of neutrino
NSI on neutrino backgrounds. The analysis was done for solar neutrinos
and in particular for the ^8B neutrino flux, for which solar
neutrino fluxes can mimic a WIMP-nucleus scattering event. We extend
upon this study by reanalyzing NSI effects in the case of lighter
mediators, and so in the light of the recently released COHERENT
data. We furthermore consider NSI effects all over solar neutrino
energies and include sub-GeV atmospheric neutrino fluxes. The former
is particularly relevant if one aims at setting e.g. background-free
sensitivities for light WIMPS.

The rest of this paper is organized as follows. In sec. <ref>
we introduce neutrino NSI, first assuming they are generated by heavy
mediators (q^2≪ m_X^2) [They might as well proceed
  through the exchange of light mediators (q^2≫ m_X^2). The
  results for this case will be presented elsewhere
  <cit.>.] and review current experimental allowed
ranges. In sec. <ref> we discuss CνNES in the
presence of NSI, calculate electron neutrino survival probabilities in
the Sun for various NSI data sets and study atmospheric vacuum
conversion probabilities. In sec. <ref> we
present our main results, namely the maximum and minimum number of
CνNES events in ton-size DM detectors (Xe based). We perform a
single-parameter analysis, c'est-à-dire we assume all NSI
parameters are zero but one and vary that coupling within its
experimental allowed range. Although constrained, this analysis enable
us to capture the main features of what should be expected from
neutrino NSI in near-future DM direct detection experiments. Finally,
in sec. <ref> we summarize and present our conclusions.



§ NEUTRINO-QUARK NON-STANDARD INTERACTIONS


NSI can be parametrized in terms of the following effective Lagrangian
(see e.g. <cit.>):

    ℒ_Eff=-√(2)G_F∑_q=u,dνγ_μ P_L(ϵ^qV 
        + ϵ^qAγ_5 )ν  qγ^μ q ,

where ν^T=(ν_e, ν_μ, ν_τ) and
ϵ^qV and ϵ^qA refer to
3×3 matrices in flavor space whose elements correspond to
vector and vector-axial couplings. Hereafter we will denote
flavor-related matrices in bold-face. This parametrization assumes
that the neutrino NSI are related to new physics at a scale
m_med≫ q, with q the characteristic energy momentum
transferred in the relevant process.  It can be as well that the new
physics is related with a scale way below the typical momentum
exchange (light mediator scenario), in which case a phenomenological
analysis requires to take into account the dynamics of the
mediator. Such models with light mediators have been discussed in, for
instance <cit.>. Here, we
focus on the intermediate possibility where CHARM and NuTeV
constraints do not apply, while an effective approach is still valid
so that the constraints on NSI parameters from the COHERENT data
analyzed in <cit.> hold. In this regime, the
constraints on the flavor-diagonal NSI couplings are substantially
weaker. For instance (see e.g. <cit.>)
  
    |ϵ_μμ^dV| < 0.042    ()
        ⟶ -0.075 < ϵ_μμ^dV  < 0.33
           ()

Accelerator constraints are irrelevant for mediator masses below
O(1 GeV). The precise value is not important here, as long as
it is larger than about 100 MeV, the maximal momentum exchange scale
we will consider.  Since <cit.> assumes
m_ med≳ 10, our analysis will be
therefore limited to
10  < m_ med < 1. While a
  priori this seems a narrow window in parameter space, it is however
for that window where the largest possible deviations—compared to SM
expectations—on the neutrino-nucleus scattering rate are expected.

The interactions in (<ref>) lead to forward coherent
scattering (order G_F interaction) and scattering (order G_F^2
interaction) processes such as CνNES. The former are responsible
for matter potentials in the Sun and in the Earth and are entirely
controlled by ϵ^qV. For the latter, the vector
current determines spin-independent processes while the axial-vector
current accounts for spin-dependent ones.  Reactions rates driven by
vector-axial currents are relatively smaller (in particular for large
nuclei) and so we will not consider them. From now on then we drop the
index V.

Constraints on NSI are derived from neutrino oscillation
<cit.>
and COHERENT data <cit.>. Let us discuss this in more
detail. The presence of the vector current in (<ref>)
induce additional diagonal and non-diagonal matter potential terms
which affect neutrino flavor evolution in the Sun and in the Earth
(depending on the neutrino energy). They produce distortions in the
neutrino oscillation probabilities that can be constrained by using
solar, KamLAND, atmospheric and long-baseline neutrino data
<cit.>. Limits on all flavor non-diagonal NSI
couplings and two combinations of diagonal couplings
e.g. ϵ_ii^q-ϵ_μμ^q (i=e,τ) can be derived
this way [Since the evolution of the neutrino flavor
  eigenstates is invariant under a constant shift of the Hamiltonian
  H → H - 𝒞𝕀 (with 𝒞 an
  arbitrary constant and 𝕀 the identity matrix), neutrino
  oscillation experiments can only constraint diagonal elements
  differences.]. For comparison we show the 90% CL limits reported in
<cit.> for the large mixing angle (LMA) case
[The combined analysis of oscillation and COHERENT data
  disfavor the so-called LMA-dark solution at the 3.1σ
  (3.6σ) level for up (down) quarks <cit.>.]  in
table <ref>.


The vector current modifies as well the CνNES cross section, and
so data from the recent measurement of the CνNES process by the
COHERENT experiment <cit.> can be used to set bounds
on the vector NSI parameters. The low-energy neutrinos
(E_ν≲ 50 MeV) used by COHERENT result from
accelerator-driven protons that strike a fixed mercury target thus
producing pions. Their decay produce a monochromatic muon neutrino
flux and muons that afterward decay, thus producing anti-muon and
electron continuous neutrino fluxes. These data therefore enable
placing bounds on ϵ^q_ei and ϵ^q_μ i
(i=e, μ, τ) <cit.>. Overall, bounds
can be placed by two independent experimental sources (order G_F and
order G_F^2 processes), but can be as well derived by combining
both, resulting in more competitive limits. Such an analysis has been
done in ref.<cit.>, and these are the limits we use in
the study of the impact of neutrino NSI in near-future DM detectors in
sec. <ref>. For reference they are displayed in table
<ref> (table to the right).



§ NEUTRINO-NUCLEUS ELASTIC SCATTERING EVENTS




 §.§ Coherent neutrino-nucleus elastic cross section and NSI


CνNES is a neutral-current process in which the momentum exchange
is of order q≲ 100 MeV. In the SM it is therefore
well-described by the effective Lagrangian

    ℒ_Eff^SM=-√(2)G_F
      ∑_q=u,d(ν_iγ_μ P_Lν_i)
      qγ^μ(g^q_V + g^q_Aγ_5)q ,

where g^q_V and g^q_A are the SM vector and axial-vector
couplings. The cross section for CνNES with spinless nuclei is
given by <cit.>

    dσ(E_ν,E_r)/dE_r=
      G_F^2/4π Q^2_SM m_N
      (1-E_r/E_r^max(E_ν))F^2(E_r) .

Here m_N is the nucleus mass and E_r the nuclear recoil energy,
whose maximum value, E_r^max, is given by

    E_r^max=2E_ν^2/m_N+2E_ν  .

Note that since E_ν≪ m_N (E_ν≲ 100 MeV for the
neutrino-nucleus scattering process to be coherent),
E_r/E_r^max≃ E_rm_N/2E_ν^2. The atomic number and
number of neutrons are encoded in
Q_SM=N-(1- 4sin^2θ_W)Z, from which one can see the
N^2 enhancement of the CνNES cross section. Finally, we use the
Helm <cit.> form factor

    F^2(q r_n)=3 j_1(q r_n)/q r_n e^-q^2 s^2/2 ,

where q can be traded to recoil energy through
q=6.92× 10^-3√(A E_r) fm^-1, j_1(x) is the
order-one spherical Bessel function of the first kind and for the
different parameters in (<ref>) we adopt the conventions
from <cit.>. Nuclear radius according to

    r_n^2=c^2 + 7/3π^2a^2-5s^2 ,

with skin thickness given by s=0.9fm, a=0.52fm and
c=(1.23 A^1/3 - 0.6)fm.

CνNES can take place with non-zero spin nuclei as well.  For
example, of the six most abundant Xe isotopes (^129Xe, ^130Xe,
^131Xe, ^132Xe, ^134Xe, ^136Xe) ^129Xe is a
J=1/2 fermion state whereas ^131Xe is a J=3/2. To properly
account for scattering with non-zero spin nuclei one has to include
contributions from the axial-vector current
in (<ref>). These contributions, being related with the spin
operator, are suppressed compared with those from the vector
current. Thus, to a fairly good approximation, the CνNES is
well-described by (<ref>) regardless of the spin of
the target nuclei.

Once NSI are introduced, the CνNES cross section is no longer
flavor diagonal and so different neutrino flavors have different cross
sections. For ν_i-N it reads <cit.>:

    dσ_ν_i/dE_r=  G_F^2/4π Q^2_ν_i m_N
      (1-E_r/E_r^max(E_ν))F^2(E_r) .

where the only difference with respect to (<ref>)
arises from Q_SM→ Q_ν_i, with

    Q_ν_i^2=4
      [
        -Q_SM/2
        +
        N(ϵ_ii^u + 2ϵ_ii^d)
        +
        Z(2ϵ_ii^u + ϵ_ii^d)
      ]^2
      +
      ∑_j≠ i
      4
      [
        N(ϵ_ij^u + 2ϵ_ij^d)
        +
        Z(2ϵ_ij^u + ϵ_ij^d)
      ]^2 .

Note that for light mediators (m_X^2≪ q^2) this is no longer true,
and a different recoil energy dependence is expected. If this is the
case, the ν-N and WIMP-nucleon recoil spectra might differ, thus
alleviating the effects of the solar and atmospheric neutrino
backgrounds <cit.>.


 §.§ Solar neutrino flavor conversion and 
  neutrino-nucleus scattering rates


Solar neutrinos (ν_e) are produced in regions near the solar core
(r≲ 0.2 R_⊙), with the exact distribution of neutrino
production determined by the SSM.  Once produced, electron neutrinos
are subject to flavor conversion, governed by the vacuum and matter
Hamiltonians according to

    id/dr|ν⟩
      =[
        1/2E_νU H_vac U^†
        +
        H_mat]|ν⟩ .

Here |ν⟩^T= |ν_e, ν_μ, ν_τ⟩^T
refers to the neutrino flavor eigenstate basis, r to the neutrino
propagation path,
U=U(θ_23)U(θ_13)U(θ_12)
is the leptonic mixing matrix (U(θ_ij) is a 3× 3
rotation matrix),
H_vac=diag(0,Δ m_21^2,Δ
m_31^2)
and in the absence of NSI
H_mat=√(2)G_F n_e(r)diag(1,0,0),
with n_e(r) the solar electron number density. Neutrino NSI induce
additional matter potentials which change the matter Hamiltonian and
thus affect neutrino flavor evolution. In full generality the
evolution equation in (<ref>) becomes

    id/dr|ν⟩
      =
      [
        1/2E_νU H_vac U^†
        +
        √(2)G_Fn_e(r)
        ∑_f=e,u,dε^f]|ν⟩ ,

where the NSI coupling matrices ε^f read

    ε^f=
       [ 1 + ε_ee^f     ε_eμ^f     ε_eτ^f;     ε_eμ^f     ϵ_μμ^f     ε_μτ^f;     ε_eτ^f     ε_μτ^f     ε_ττ^f;            ] ,

with ε_ij^f(r)=Y_f(r)ϵ_ij^f (f=e,u,d) and
Y_f(r)=n_f(r)/n_e(r). The up- and down-quark relative abundances can
be written in terms of the neutron relative abundance according to

    Y_u = 2 + Y_n    and     Y_d = 1 + 2 Y_n ,

with the neutron number density in turn calculated from the ^4He and
^1H mass fractions (metallicity in the Sun amounts to less than
1.5% so it can be safely neglected), namely

    n_n(r)≃X(^4He)/2X(^1H)+X(^4He) .

Since we are interested in CνNES, from now on we set
ϵ^e_ij=0 and following
refs. <cit.> we do not consider
simultaneous contributions from up- and down-quark couplings. For all
the SSM related quantities we use those from the BS05 SSM
<cit.>.

To a fairly good approximation, flavor conversion probabilities can be
calculated in the mass dominance limit, Δ m_31^2→∞.
In this limit, and with vanishing CP-violating phases, neutrino
propagation takes place in the basis
|ν̃⟩=U_13^TU_23^T|ν⟩,
which neglecting terms proportional to sinθ_13 corresponds
to the following neutrino eigenstates:

    ν̃_e≃ν_e ,    ν̃_μ= cosθ_23ν_μ - sinθ_23ν_τ ,    ν̃_τ≃sinθ_23ν_μ + cosθ_23ν_τ .

In this basis, which we dub as propagation basis,
ν̃_e-ν̃_τ and ν̃_μ-ν̃_τ
mixing is of order
ξ_ij=G_Fn_q(r)ϵ_ij^q/(Δ m_31^2/E_ν).  Taking
n_q≃ 10^26 cm^-3 (up-quark number density at about
r=0.1 R_⊙), E_ν≲ 20MeV and
Δ m_31^2=2.55× 10^-3 eV^2 (BFPV for normal
ordering) <cit.>, one can see that even for order-one
NSI couplings ξ_ij≪ 1.  This means that in the mass dominance
limit, ν̃_τ is decoupled and so neutrino flavor evolution
can be studied in the two-flavor approximation with electron neutrino
flavor conversion determined by
1-⟨𝒫_ee(E_ν)⟩. Thus, in what follows we
discuss the calculation of the averaged survival probability
⟨𝒫_ee(E_ν)⟩. First of all, due to the
quark and electron densities 𝒫_ee depends not only on
E_ν but also on the neutrino propagation path
r. 𝒫_ee(E_ν,r) can be written as <cit.>

    𝒫_ee(E_ν,r)=cos^4θ_13 𝒫_eff(E_ν,r)
      +
      sin^4θ_13 ,

where the r dependence is introduced by the effective probability
given by <cit.> [We have checked that
  γ^-1≪ 1, which guarantees neutrino adiabatic propagation,
  and so have taken the level-crossing probability P_c→ 0.]

    𝒫_eff(E_ν,r)=
      1+cos2θ_M(r)cos2θ_12/2 .

The mixing angle in matter is calculated from the diagonalization of
the 2× 2 Hamiltonian

    H=1/4E_ν[ -Δ m_21^2 cos2θ_12 + A  Δ m_21^2 sin2θ_12 + B;  Δ m_21^2 sin2θ_12 + B  Δ m_21^2 cos2θ_12 - A ] ,

where

    A=4√(2)E_ν G_F n_e(r) 
      [cos^2θ_13/2 - Y_q(r)ε_D] ,    
      B=4√(2)E_ν G_F n_e(r) Y_q(r)ε_N  .

Note that in the limit ϵ_ij^q=0 and cosθ_13=0,
A reduces to the SM term and B vanishes. The parameters
ε_D and ε_N result from the rotation from the
flavor to the propagation basis and read
<cit.>:

    ε_D   =
      -c_13^2/2ϵ_ee^q
      +[c_13^2 - (s_23^2 - s_13^2 c_23^2)]/2ϵ_μμ^q
      +(s_23^2 - c_23^2s_13^2)/2ϵ_ττ^q + s_13c_13s_23ϵ_eμ^q
      
          
      + s_13c_13c_23ϵ_eμ^q
      - c_23s_23ϵ_μτ^q ,
      
    ε_N   = 
      -s_13c_23s_23ϵ_μμ^q
      + s_13c_23s_23ϵ_ττ^q
      + c_13c_23ϵ_eμ^q
      - c_13s_23ϵ_eτ^q
      + s_13(s_23^2 - c_23^2)ϵ_μτ^q .


Here we have adopted the notation cosθ_ij=c_ij and
sinθ_ij=s_ij. From eq. (<ref>) the mixing
angle in matter can be straightforwardly calculated:

    cos2θ_M(r) = Δ m_12^2cos2θ_12-A/√((Δ m_12^2cos2θ_12-A)^2
          +
      (Δ m_12^2sin2θ_12+B)^2) .

Eqs. (<ref>) and (<ref>) combined with
eqs. (<ref>)-(<ref>) allow the determination of
𝒫_ee(E_ν,r) in terms of neutrino oscillation
parameters, electron and quark number densities and NSI
parameters. The averaged survival probability is then obtained by
integrating over r taking into account the distribution of neutrino
production in the Sun <cit.>:

    ⟨𝒫_ee(E_ν)  ⟩=
      ∑_αΦ_α(E_ν)∫_0^1
         drρ_α(r) 𝒫_ee(E_ν,r)/∑_αΦ_α(E_ν) ,

where Φ_α(E_ν) stands for neutrino fluxes and
ρ_α(r) for the distribution of neutrino production,
α=pp,^7Be (ground and excited state), pep, hep,
^8B, ^13N, ^15O, ^17F.
In practice, to calculate ⟨𝒫_ee(E_ν) ⟩ we
have followed a single-parameter analysis, varying only one NSI
coupling at a time in the ranges allowed by the constraints derived
from neutrino oscillation and COHERENT data
<cit.>. These ranges have been split according to
ϵ_ij^q|_n+1=ϵ_ij^q|_n+δϵ_ij^q
with
δϵ_ij^q=
(ϵ_ij^q|_max-ϵ_ij^q|_min)/10.
And the neutrino energy, E_ν, has been varied from 0.145 MeV to
14.5 MeV. The neutrino oscillation parameters have been fixed to their
BFPVs according to ref.<cit.>. For each NSI parameter
we thus have generated 11 data sets (11 averaged survival
probabilities), including the SM case (ϵ_ij^q=0).

Focusing on up-quark NSI couplings (results for down-quark parameters
resemble those of the up-quark case), the averaged survival
probability deviates the most from the SM expectation for
flavor-diagonal NSI and ϵ_eτ^u, being more pronounced
for the former than for the latter. Fig. <ref> shows the
result for ⟨𝒫_ee⟩ for the cases
ϵ_ee^u≠ 0 and ϵ_μμ^u≠ 0 (results for
ϵ_ττ^u≠ 0 are pretty close to those found for
ϵ_μμ^u≠ 0).

It can be seen that the range of variation of
⟨𝒫_ee⟩ has about the same amplitude in both
cases. However its behavior with increasing values of the NSI
couplings is rather different. For ϵ_ee^u≠ 0, the
averaged survival probability decreases with increasing values of
ϵ_ee^u and its largest value (as a function of E_ν) is
indeed found when ϵ_ee=0. For ϵ_μμ^u≠ 0,
it is the other way around. Increasing values of the parameter leads
to larger ⟨𝒫_ee⟩. This behavior can be
readily understood as follows. In the SM case limit,
ϵ_ij^u=0, eq. (<ref>) reduces to

    cos2θ_M=Δ m_21^2cos2θ_12
        - 2√(2)E_ν G_Fn_e c_13^2/√((Δ m_21^2cos2θ_12
        - 2√(2)E_ν G_Fn_e c_13^2)^2
      +Δ m_21^4sin^22θ_12) .

The presence of the ϵ_ee^u and ϵ_μμ^u
couplings changes this equation by shifting either the first term in
the denominator (numerator) or both. As can be seen from
eqs. (<ref>)-(<ref>) the shift in the numerator
of eq. (<ref>) due to ϵ_ee^u is
always negative. Thus, as ϵ_ee^u increases cos2θ_M
decreases and accordingly ⟨𝒫_ee⟩ decreases
(see eqs. (<ref>) and (<ref>)). The shift
due to ϵ_μμ^u is instead positive thus leading to
larger values of cos2θ_M as ϵ_μμ^u increases.

As we have pointed out, in the mass dominance limit solar neutrino
flavor conversion can be reduced to a two-flavor oscillation problem,
in which a mainly |ν_e⟩ state oscillates to a
|ν̃_μ⟩ state which is a superposition of
|ν_μ⟩ and |ν_τ⟩ flavor eigenstates. Thus, in
an Earth-based detector, neutrino-nucleus scattering will take place
with either |ν_e⟩ or |ν̃_μ⟩ and so the
calculation of neutrino-nucleus scattering rates has to be done in the
propagation basis as well. Rotating the effective Lagrangian
in (<ref>) to this basis, after dropping the
axial-vector coupling matrix, results in

    ℒ_Eff=-√(2)G_F∑_q=u,dν̃γ_μ P_Lϵ^q ν̃  qγ^μ q ,

with the NSI couplings in the propagation basis related with the NSI
couplings in the flavor basis through the following transformation

    ϵ^q=
      U(θ_13)^T U(θ_23)^T ϵ^q U(θ_23) U(θ_13) .

The differential cross sections for ν̃_e-N≃ν_e-N and
ν̃_μ-N can then be obtained from (<ref>)
and (<ref>) by ϵ_ij^q→ϵ̃_ij^q.
Results for the cases ϵ_ee^u≠ 0 and
ϵ_μμ^u≠ 0 are shown for both neutrino flavors in
fig. <ref>. As can be seen, deviations on
either dσ_ν_e/dE_r or dσ_ν̃_μ/dE_r are
possible. The presence of NSI can enhance but also can deplete the
flavored cross sections, thus potentially leading to substantial
deviations of the expected neutrino-nucleus scattering event rate.


Assuming a 100% efficiency and with the oscillation probabilities and
the ν-N differential cross sections calculated in the propagation
basis, the ν-N event rate as a function of the recoil energy E_r
can be written according to <cit.>:

    dR_ν^⊙/dE_r=N_T ∑_α∫_E_min^∞ Φ_α(E_ν)
      [
       ⟨𝒫_ee⟩dσ_ν_e/dE_r
        +
         (1- ⟨𝒫_ee⟩)
        dσ_ν̃_μ/dE_r] dE_ν .

Here N_T=N_A/A (N_A=6.022× 10^29 ton^-1),
E_r^min=√(m_N E_r/2) and we include neutrinos emitted in
the pp and CNO chains, α=pp,^7Be (ground
  and excited state), pep, hep, ^8B,
^13N,
^15O,
^17F,
as in (<ref>). The neutrino fluxes are in
units of cm^-2 year^-1 keV^-1,
so the differential rates are in units of
ton^-1 year^-1 keV^-1
and the total event rates, R_ν^⊙,
are then measured in ton^-1 year^-1.
The latter are calculated from (<ref>) by
integration over E_r
from E_r^th
to 100 keV, with E_r^th⊂
[10^-3,10^2] keV and varied in steps
δ E_r=(E_r^n+1-E_r^n)/50.


 §.§ Atmospheric neutrino-nucleus scattering rates


In contrast to solar neutrinos, matter effects for sub-GeV atmospheric
neutrinos are negligible and flavor conversion can be well described
by vacuum oscillations <cit.>. This can be readily
seen by considering the PREM model <cit.>, for which
⟨
n_e^⊕⟩=3.25× 10^24 cm^-3 and ⟨
Y_n^⊕⟩=1.137 both in the Earth core. With these values,
using the BFPV for Δ
m_31^2 <cit.> and taking
E_ν=100 MeV
so to minimize the vacuum term contribution one gets (see
eq. (<ref>))

    H_mat/H_vac≃
      0.2×(2.55× 10^-3/Δ m_31^2)
      (E_ν/100 MeV)
      (n^⊕_e/3.25× 10^24 cm^-3)
      (Y_u/3.137) ϵ_ij^u .

Note that this is the largest contribution from H_mat
to the neutrino flavor evolution equation, thus implying that matter
effects for sub-GeV neutrinos can be fairly ignored.


Since only vacuum oscillations are relevant, the full three-flavor treatment
is straightforward. In this case the oscillation probabilities can be
further simplified by considering the limit Δ
m_21^2→ 0. They read <cit.>

    𝒫_ij(E_ν,L)=4|U_i3|^2|U_j3|^2
      sin^2(Δ m_31^2/4E_ν L) ,

where L
refers to the distance covered by a neutrino from production at a height
h (perpendicular to the surface) to detection and is given by:

    L=√((R_⊕ + h)^2-R_⊕^2sin^2θ_z)
      -
      R_⊕cosθ_z ,

with θ_z
the zenith angle of the neutrino trajectory. The oscillation
probability in (<ref>) is therefore a
function of production height h
and zenith angle θ_z.
The distribution of vertical production height for sub-GeV neutrinos
was calculated in <cit.> and its mean was found to
be ⟨
h⟩≃
18 km. Thus, rather than using the distribution we fix
h
to this value and integrate over the zenith angle. The oscillation
probabilities relevant for the determination of the ν-N
event rate can then be written as

    𝒫_ij(E_ν)=1/π∫_-1^1 dcosθ 𝒫_ij(E_ν,L) .

As expected, neutrino flavor oscillations in this case are only
relevant for ν_μ-ν_τ
conversion. The result is shown in fig. <ref>.

Since we set CP-violating phases to zero, the differential event rates
associated with ν_e+ν̅_e and ν_μ+ν̅_μ fluxes can
be written according to

    dR_ν_e+ν̅_e^Atm/dE_r   =   N_T∫_E_min^∞Φ_ν_e+ν̅_e(E_ν)
      [
      (1-𝒫_eμ-𝒫_eτ)dσ_ν_e/dE_r
      +
      𝒫_eμdσ_ν_μ/dE_r
      +
      𝒫_eτdσ_ν_τ/dE_r] ,
      
    dR_ν_μ+ν̅_μ^Atm/dE_r   =   N_T∫_E_min^∞Φ_ν_μ+ν̅_μ(E_ν)
      [
      𝒫_eμdσ_ν_e/dE_r
      +
      (1-𝒫_eμ-𝒫_μτ)
      dσ_ν_μ/dE_r
      +
      𝒫_μτdσ_ν_τ/dE_r] .

In (<ref>) we assume time reversal
invariance. We use the atmospheric neutrino fluxes from
ref. <cit.>, derived with the FLUKA Monte Carlo
simulation package <cit.> and which provides fluxes for
ν_e, ν_μ and the corresponding anti-neutrinos. The full
differential ν-N event rate is obtained from
(<ref>) and (<ref>)
and the total event rate as in the solar case, integrating over E_r
from E_r^th to 100 keV, with
E_r^th⊂ [10^-3,10^2] keV and varied in steps
δ E_r=(E_r^n+1-E_r^n)/50. The flavored cross sections are
then given by (<ref>) for which one finds sizable
deviations from the SM only for diagonal NSI couplings, as expected
given the tight constraints on the non-diagonal, particularly
ϵ_eμ^u and ϵ_μτ^u. The result for
diagonal couplings is shown in fig. <ref>, in which all
couplings but one have been assumed to be zero.



§ NSI AT TON-SIZE DM XE DETECTORS


We now turn to the discussion of the impact of neutrino NSI on the
expected number of events in a Xe-based DM detector. For concreteness
we assume an exposure of
ℰ=1 ton×year. For practical reasons,
and using the same approach we used in the calculation of the averaged
survival probability in the Sun (see sec. <ref>),
we perform a single-parameter analysis, considering only up-type
couplings (results for down-type couplings resemble those found for
the up-quark case). We first calculate the neutrino-nucleus
differential event rate for ten equally-spaced values of the
corresponding NSI parameter:
ϵ_ij^u|_a+1=ϵ_ij^u|_a+δϵ, with
δϵ=
(ϵ_ij^u|_max-ϵ_ij^u|_min)/10.
From these results we then identify the parameters that
maximize/minimize the recoil spectra. And for these parameters we
calculate the neutrino-nucleus scattering events. For solar neutrinos
fluxes we employ the results derived in sec. <ref>,
eq.  (<ref>) combined with the averaged survival
probabilities. For atmospheric neutrino fluxes we use the results from
sec. <ref>. The results of our calculation are
displayed in fig. <ref>, which we regard
as the main result of this paper.

One can see that even those parameters that are forced to be small by
phenomenological constraints can have a sizable impact on the number
of events. Departures from the SM value (horizontal line at 1 in
fig. <ref>) are energy dependent and are
particularly pronounced in the solar-atmospheric transition, at
E_r≃ 5.8 keV when hep neutrinos reach their energy tail.
Note that in the absence of NSI propagation effects, the SM and NSI
neutrino-nucleus scattering rates differ only by a global numerical
factor, Q_SM→ Q_ν_i (see eq. (<ref>)).  Thus,
this recoil energy dependence is a consequence of NSI propagation
effects.  Neutrino NSI not only lead to enhancements of the total
number of events but can as well produce depletions (as firstly noted
in ref. <cit.>), with both being relevant for DM direct
detection experiments. The former places stronger limits on the
sensitivity one can reach in a multi-ton DM detector, while the latter
enables measurements of WIMP-nucleon cross sections that otherwise
would be challenging reaching.


Overall, one observes that no matter the coupling a deviation from the
SM expectation is always possible. However, diagonal couplings lead to
the largest deviations which can either be enhancements or depletions
of the event rate. Particularly important for solar neutrinos, are
deviations in the energy range where the background is dominated by
^8B, since it is in that range where solar neutrinos can mimic a
WIMP-nucleus signals <cit.>. In there, our findings
show that ϵ_μμ^u leads to enhancements of
R_ν^NSI which can exceed the SM value by about 40%. For
ϵ_ττ^u, enhancements are slightly smaller, but still
substantial. For ϵ_ee^u they reach values of order 2 or so,
being for this coupling the largest possible. This is somehow
expected, as is for this coupling that bounds are less severe. Note
that compared to what was found in ref. <cit.> this
enhancement is larger, with the difference just due to the range
employed. In ref. <cit.>
ϵ_ee^u|_max=0.45, value derived from from
oscillation and deep inelastic scattering data <cit.>,
while in our case ϵ_ee^u|_max=0.6. This difference,
although small, leads to a sufficiently large enhancement of the
neutrino-nucleus cross section (see
fig. <ref>) that accounts for the larger
deviation we found.

Since the largest deviations are found for diagonal NSI we focus in
these cases, discussing in a bit more detail the results for
ϵ_ee^u and ϵ_μμ^u. Maximization or
minimization of R_ν is determined by the way in which the NSI
parameters affect detection and propagation. The latter relevant only
in the solar sector. For ϵ_ee^u we found
R_ν^NSI|_max at ϵ_ee^u|_max=0.6
for both, the region where recoils are controlled by solar neutrinos
fluxes (E_r≲ 5.8 keV) and the region governed by atmospheric
neutrinos (E_r≳ 5.8 keV).  This can be understood as
follows. In the solar region dσ_ν̃_μ/dE_r does not
depend on ϵ_ee^u and so it keeps always its SM value (see
fig. <ref>). The ν̃_e-N cross
section instead reaches its maximum exactly at
ϵ_ee^u|_max. The survival probability behaves in the
opposite way, as ϵ_ee^u increases 𝒫_ee
decreases globally, reaching minimum values exactly at
ϵ_ee^u|_max (see fig. <ref>, left
graph). This effect tends to deplete the contribution of the
ν̃_e-N scattering process to the total number of
events. However, since for ϵ_ee^u=0.6 the survival
probability decreases at most 20% from its SM value, the overall
effect is not sufficiently strong resulting in
R_ν^NSI|_max at ϵ_ee^u|_max.
R_ν^NSI|_min is found at about the point where
dσ_ν̃_e/dE_r reaches its minimum,
ϵ_ee^u≃ 0.2. In the atmospheric sector, due
𝒫_eμ, 𝒫_eτ≪ 1 the event rate is
mainly controlled by the ν_e-N cross section, which can be largely
modified by NSI effects (see figs. <ref> and
<ref>). This in turn is reflected in the non-depletion of
R_ν^NSI when passing from the solar to the atmospheric
regions.

The parameter dependence in the case of ϵ_μμ^u is
rather different.  For solar neutrinos, R_ν^NSI|_max
and R_ν^NSI|_min are found where the
ν̃_μ-N cross section is maximized/minimized (see
fig. <ref>, right graph). When solar fluxes
are no longer relevant and atmospheric kick in, the maximum expected
number of events decreases. From (<ref>) one
can see that contributions from ν_e+ν̅_e fluxes are mainly
SM, sizable deviations are due only to ν_μ+ν̅_μ
fluxes. They are however not as large as in the case of
ϵ_ee^u due to the sizable value of 𝒫_μτ
which tends to deplete the value of the second term in
eq. (<ref>), thus leading to small
R_ν^NSI|_max.

The parameter dependence in the case ϵ_ττ^u is rather
similar to the one found in the ϵ_μμ^u scenario in the
solar sector, but in the atmospheric region the behavior is
different. From eqs. (<ref>) and
(<ref>) this is somehow
expected. Contributions from electron neutrino fluxes barely deviate
from SM values, substantial deviations can only arise from muon
neutrino fluxes, with the ν_τ-N cross section weighted by the
ν_u-ν_τ oscillation probability which is of order 0.35 (see
fig. <ref>). This results in slightly larger deviations
that those found in the ϵ_μμ^u case. Non-diagonal NSI
couplings are subject to tighter constraints, thus resulting in
smaller effects. The exception being the case of the
ϵ_eτ^u parameter, for which we found no depletion and
enhancements of order 40%, inline with what was reported in
<cit.>.

The exact value (range) of the mediator mass for which our analysis is
valid depends—of course—on the recoil energy and hence on the
neutrino sector, solar or atmospheric. For the
ϵ_μμ^u≠ 0 case
(fig. <ref> left top panel), we show the
minimum mass that the mediator should have. In the solar sector masses
below ∼ 120 MeV suffice, with the exact value depending on the
specific recoil energy. For recoil energies associated with
atmospheric neutrinos fluxes, m_med≳ 120MeV. Thus,
while for solar neutrinos most of the mediator mass range is
available, for atmospheric neutrinos the analysis is limited to a
narrow window. However, as we have already pointed out is for this
window where NSI effects may be expected to be maximized.


 §.§ Implications for background-free sensitivities


CνNES leads to recoil energy spectra which are closely degenerate
with that induced by WIMP-nucleus scattering. As a result, once a
direct detection experiment becomes sensitive to the neutrino
background there is a limit on the WIMP cross section
(σ_DM-n) that such experiment can explore, due to the
inability to distinguish a neutrino- from a WIMP-induced signal. This
limit defines the so-called neutrino “floor”.

With the CνNES process enhanced (diminished) by the
neutrino-nucleus NSI, one expects the neutrino floor to be more (less)
severe. Thus, in order to quantify the impact of NSI on the reach of
next-generation DM detectors, we have quantified the 90% CL
sensitivities achievable. To do so, we have followed the method
introduced in ref. <cit.> which enables the
determination of the smallest WIMP cross section that can be measured
for any WIMP mass in the presence of one-neutrino background event. It
is worth pointing out that the limit on the WIMP cross section derived
this way is just an estimate of what one should really expect.  The
precise determination of the saturation point, where the WIMP cross
section one can access does not reduce with increasing (feasible)
exposure, requires abandoning the background-free paradigm and using
instead a likelihood analysis <cit.>. Here, however,
we aim just at a preliminary understanding of the possible impact that
NSI might have on direct detection sensitivities, and therefore stick
to a simple background-free analysis.


We calculate the WIMP cross section as a function of the WIMP mass
by fixing the WIMP-nucleus scattering event rate, R_DM,
to 2.3 (90% CL sensitivity limits), namely

    σ_DM-n=2.3/ℰ(E_r^th)
      ∫_E^th_r^E_r^max . dR_DM/dE_r|_σ_DM-n=1 .

Here the WIMP event rate, which depends upon astrophysical, nuclear
physics and particle physics parameters <cit.> reads

    dR_DM/dE_r=ℰ ρ_DM σ_DM-N/2 m_DM μ_N^2
      F^2(E_r)
      k_0/k1/(π v_0^2)^3/2 ∫^v_esc_v_min 
      e^-|v+v_lab|^2/v_0^2 d^3v/|v| ,

where ℰ is the exposure, ρ_DM is the DM
density in the local halo, σ_DM-N is the DM-nucleus
cross section, μ_N is the DM-nucleus reduced mass, v_0 is the
most likely speed of the DM in the halo, F^2(E_r) is the form factor
given in (<ref>) and v is the velocity of
the DM particles. The velocity distribution is Maxwellian and
truncated at the DM escape velocity:
|v+v_lab|=v_esc,
with v_lab the laboratory (Earth) velocity
relative to the DM local halo. The factor k_0/k is a normalization
factor obtained by integrating the Maxwellian DM velocity distribution
from 0 to v_esc and is given by

    k_0/k=
      [
        erf(v_esc/v_0)
        -
        2/√(π)v_esc/v_0e^-v_esc^2/v_0^2]^-1 .

Rather than integrating numerically the velocity distribution in
(<ref>) we have used its exact analytic form,
namely <cit.>

    ∫_v_min^v_esc 
      e^-|v+v_lab|^2/v_0^2 d^3v/|v|
      =
      2π v_0^2 {√(π)/4v_0/v_lab[
        erf(v_min+v_lab/v_0)
        -
        erf(v_min-v_lab/v_0)
        ]
        -
        e^-v_esc^2/v_0^2} ,

with v_min=√(2E_r/m_DM r) × c and
r=4m_DMm_N/(m_DM+m_N)^2. We have written as well the
DM-nucleus cross section in terms of the DM-nucleon cross section

    σ_DM-N=μ_N^2/μ_n^2 A^2 σ_DM-n ,

with μ_n the DM-nucleon reduced mass. For the different
astrophysical quantities we have used: ρ_DM=0.3
GeV/c^2/cm^3, v_0=220 km/s, v_lab=232 km/s and
v_esc=544 km/s. For the integration we varied E_r^th
from 10^-3 to 10^2 keV (with 10^3 partitions) and took

    E_r^max=m_DM/2(v_esc+v_lab/c)^2
      r/2 ,

with r given after
eq. (<ref>). The exposure has been
tuned so to assure that each E_r^th isocontour in the
σ_DM-n-m_DM plane contains one-neutrino
event. This has been done by calculating the exposure as follows:

    ℰ(E^th_r)=1 neutrino event/∫_E^th_RdR_ν/dE_r .


The result is shown in fig. <ref>,
where we have plotted the SM sensitivities as well and included
current (future) exclusion limits from XENON1T and PandaX-II
<cit.> (XENONnT and PandaX-4T
<cit.>). The sensitivities shown in the plot refer to those
obtained by calculating the exposure in (<ref>) for the
values of ϵ_ee^u that maximize/minimize dR_ν/dE_r
(lower and upper curves in fig. <ref>,
left-top graph). The result in
fig. <ref> show that the effects of
neutrino-nucleus NSI, if present, work either way. They can worsen
sensitivities or improve them. Note, that despite the effect being
small it might have consequences. For example, if the NSI-induced
depletion is effective (lower purple curve) the background-free
sensitivities at XENONnT will be less affected. If instead the
NSI-induced enhancement is present, XENONnT or PandaX-4T could start
measuring neutrino scattering events sooner than expected.


§ CONCLUSIONS


The identification of a DM signal in next-to-next generation
(multi-ton scale) detectors will be seriously challenged by
irreducible solar and atmospheric neutrino backgrounds. Identifying
all the features of this background is therefore mandatory to
determine the implications it will have and the challenges it will
pose. With this in mind, in this paper—following
<cit.>—we have entertained the possibility that
neutrinos are subject to NSI and have studied the impact they would
have on the expected number of neutrino-nucleus scattering events in a
Xe-based detector assuming an exposure of
1 ton·year. We have as well determine the
background-free sensitivities achievable if these interactions are
present, showing that they can either worsen or improve them, although
only slightly in both cases.

Despite the recent observation of the CνNES process by the
COHERENT experiment and its consistency with the SM expectation, NSI
couplings can still be relatively large, even when the COHERENT
scattering data is combined with data from neutrino oscillation
experiments. Thus, by using current constraints applicable in
scenarios where neutrino NSI arise from mediators with masses in the
range [10, 10^3] MeV, we have calculated the expected number
of neutrino-nucleus scattering events in the presence of NSI including
propagation and detection effects and extending the analysis over the
solar and atmospheric neutrino sectors.  In the former case, we have
included all pp and CNO solar neutrino fluxes, while in the latter
case we have included ν_e, ν_μ, ν̅_e and ν̅_μ
fluxes for neutrino energies up to ∼ 100 MeV. To simplify the
computational burden, we have carried out a single-parameter analysis
in which all couplings were put to zero but one and repeated the
analysis for the remaining NSI couplings. Since we used constraints
from <cit.>, we analyzed up- and down-quark couplings
independently and presented only the results for the former, as the
results for the latter are rather similar.

For the solar sector, we have studied propagation in the mass
dominance limit and have taken into account the distribution of
neutrino production within the Sun, for the calculation of the
averaged neutrino survival probability. In the atmospheric region,
Earth matter effects are negligible due to the largest relevant
neutrino energies being below ∼ 100 MeV. In this sector then, NSI
affect the expected number of neutrino-nucleus scattering events only
through detection. We have anyway included vacuum oscillation
probabilities (properly integrated over the zenith angle) in the
determination of the number of scattering events.

Our main results are given by fig. <ref>
and fig. <ref> and can be summarized
as follows. Largest deviations from the SM expectation are more
pronounced for flavor-diagonal NSI and can operate either enhancing
the number of events or depleting them. For the most relevant recoil
energy interval in the solar sector (that in which the neutrino flux
is determined by ^8B neutrinos in which the background can mimic
WIMP-nucleon scattering events) the largest enhancement is found for
ϵ_ee^u and can exceed the SM expectation up to a factor
∼ 2. For ϵ_μμ^u and ϵ_ττ^u,
enhancements are smaller but still sizable. Depletions in that energy
interval can amount to about 30-40% for all diagonal NSI
parameters. Moving to higher neutrino energies we found that when
atmospheric neutrino fluxes dominate the background, enhancements as
large as twice the SM expectation are possible when ϵ_ee^u
is present, in all other cases they are below this value but still are
sizable. Depletions in the atmospheric sector can amount to about
50% if ϵ_ττ^u is at work. For the off-diagonal
couplings ϵ_eμ^u and ϵ_eτ^u we did not find
any depletion, but we did for ϵ_μτ^u although rather
small (order 5%) and so probably not experimentally
measurable. Thus, if the neutrino-nucleus event rate for energies
dominated by atmospheric neutrino fluxes is found to be below SM
expectations—arguably—that will support the presence of new
physics featuring flavor-diagonal NSI, assuming one can disentangle the
new interactions above any other effect.

Deviations, though small, might have an effect on the sensitivities
achievable. To estimate the extent of such effects we have calculated
the background-free sensitivities. We have shown that neutrino-nucleus
NSI, if present, may worsen sensitivities or improve them. The effect
is small, but might have observable consequences. We have stressed
that if neutrino NSI are present either XENONnT or PandaX-4T, among
others, will start observing neutrino-induced nuclear recoil events
sooner than expected.


§ ACKNOWLEDGMENTS

D.A.S. would like to thank Atri Bhattacharya, Jean-Rene Cudell, Thomas
Hambye, Julian Heeck, Concha Gonzalez-Garcia, Florica Stancu and
Mariam Tortola. In particular, Concha Gonzalez-Garcia for discussions
regarding refs. <cit.>, Mariam
Tortola for discussions on neutrino flavor oscillations in the Sun and
Jean-Rene Cudell and Florica Stancu for useful conversations on
nuclear physics. He would like also to thank the “Service de Physique
Théorique” of the “Université Libre de Bruxelles” and the STAR
research unit of the “Université de Liège” for the warm
hospitality during the early stages of this work and partial financial
support. This research has been supported by F.R.S.-FNRS through the
program “Bourses de Séjour Scientifique” and by the Chilean grant
“Unraveling new physics in the high-intensity and high-energy
frontiers”, Fondecyt No 1171136. NR was funded by proyecto FONDECYT
Postdoctorado Nacional (2017) num. 3170135.  The work of M.T.  is
supported by the FNRS-FRS, the Belgian Federal Science Policy Office
through the Interuniversity Attraction Pole P7/37 and the Belgian
Institute Interuniversitaire des Sciences Nucléaires (IISN).


