A Direct Detection View of the Neutrino NSI Landscape

TL;DR

This paper investigates how direct detection experiments, especially xenon-based detectors, can explore non-standard neutrino interactions by analyzing solar neutrino scattering, offering a new approach complementary to existing methods.

Contribution

It introduces a novel NSI parametrisation that includes electron contributions and visualises experimental complementarity, applying it to recent and future xenon detector data.

Findings

First limits on NSI from XENONnT and LUX-ZEPLIN data

Projections show future experiments can significantly constrain NSI

New software package SNuDD developed for analysis

Abstract

In this article, we study the potential of direct detection experiments to explore the parameter space of general non-standard neutrino interactions (NSI) via solar neutrino scattering. Due to their sensitivity to neutrino-electron and neutrino-nucleus scattering, direct detection provides a complementary view of the NSI landscape to that of spallation sources and neutrino oscillation experiments. In particular, the large admixture of tau neutrinos in the solar flux makes direct detection experiments well-suited to probe the full flavour space of NSI. To study this, we develop a re-parametrisation of the NSI framework that explicitly includes a variable electron contribution and allows for a clear visualisation of the complementarity of the different experimental sources. Using this new parametrisation, we explore how previous bounds from spallation source and neutrino oscillation experiments are impacted. For the first time, we compute limits on NSI from the first results of the XENONnT and LUX-ZEPLIN experiments, and we obtain projections for future xenon-based experiments. These computations have been performed with our newly developed software package, SNuDD. Our results demonstrate the importance of using a more general NSI parametrisation and indicate that next generation direct detection experiments will become powerful probes of neutrino NSI.

Tables

Table 1. Table 1: Solar neutrino rates relevant for our Borexino analysis. Shown are the measured rates from the Phase-II run of Borexino Borexino:2017rsf , our calculated theoretical rates, and the assumed fractional uncertainties in our calculation.

Population	Phase-II Rate [${\left(100\mathrm{ton}\mathrm{day}\right)}^{-1}$]	Theoretical Rate [${\left(100\mathrm{ton}\mathrm{day}\right)}^{-1}$]	Fractional Uncertainty
$pp$	$134\pm 10$	133	$1.1\%$
${}_{}{}^7\mathrm{Be}$	$48.3\pm 1.1$	48.5	$5.8\%$
$pep$	$2.43\pm 0.36$	2.78	$1.5\%$

Table 2. Table 2: 90 % percent 90 90\% CL allowed intervals for NSI in the up quark, down quark, proton, and neutron directions. Shown are the results from our analysis of the LZ 2022 data LZ:2022ufs and those of the global fit study of Ref.  Coloma:2019mbs . Note that the latter do not quote NSI in the neutron direction.

	LZ 2022 (this work)	Global Fits Coloma:2019mbs
${\epsilon }_{ee}^u$	$[-0.545,\mathrm{\hspace{0.17em}1.222}]$	$[-0.031,\mathrm{\hspace{0.17em}0.476}]$
${\epsilon }_{\mu \mu }^u$	$[-0.971,\mathrm{\hspace{0.17em}1.397}]$	$[-0.029,\mathrm{\hspace{0.17em}0.068}]\oplus [0.309,\mathrm{\hspace{0.17em}0.415}]$
${\epsilon }_{\tau \tau }^u$	$[-0.645,\mathrm{\hspace{0.17em}1.598}]$	$[-0.029,\mathrm{\hspace{0.17em}0.068}]\oplus [0.309,\mathrm{\hspace{0.17em}0.414}]$
${\epsilon }_{e\mu }^u$	$[-0.630,\mathrm{\hspace{0.17em}0.679}]$	$[-0.048,\mathrm{\hspace{0.17em}0.020}]$
${\epsilon }_{e\tau }^u$	$[-0.721,\mathrm{\hspace{0.17em}0.558}]$	$[-0.077,\mathrm{\hspace{0.17em}0.095}]$
${\epsilon }_{\mu \tau }^u$	$[-1.120,\mathrm{\hspace{0.17em}0.518}]$	$[-0.006,\mathrm{\hspace{0.17em}0.007}]$
${\epsilon }_{ee}^d$	$[-0.540,\mathrm{\hspace{0.17em}1.084}]$	$[-0.034,\mathrm{\hspace{0.17em}0.426}]$
${\epsilon }_{\mu \mu }^d$	$[-0.863,\mathrm{\hspace{0.17em}1.233}]$	$[-0.027,\mathrm{\hspace{0.17em}0.063}]\oplus [0.275,\mathrm{\hspace{0.17em}0.371}]$
${\epsilon }_{\tau \tau }^d$	$[-0.576,\mathrm{\hspace{0.17em}1.241}]$	$[-0.027,\mathrm{\hspace{0.17em}0.067}]\oplus [0.274,\mathrm{\hspace{0.17em}0.372}]$
${\epsilon }_{e\mu }^d$	$[-0.542,\mathrm{\hspace{0.17em}0.635}]$	$[-0.050,\mathrm{\hspace{0.17em}0.020}]$
${\epsilon }_{e\tau }^d$	$[-0.655,\mathrm{\hspace{0.17em}0.455}]$	$[-0.076,\mathrm{\hspace{0.17em}0.097}]$
${\epsilon }_{\mu \tau }^d$	$[-0.982,\mathrm{\hspace{0.17em}0.461}]$	$[-0.006,\mathrm{\hspace{0.17em}0.007}]$
${\epsilon }_{ee}^p$	$[-1.805,4.195]$	$[-0.086,\mathrm{\hspace{0.17em}0.884}]\oplus [1.083,\mathrm{\hspace{0.17em}1.605}]$
${\epsilon }_{\mu \mu }^p$	$[-3.330,4.791]$	$[-0.097,\mathrm{\hspace{0.17em}0.220}]\oplus [1.063,\mathrm{\hspace{0.17em}1.410}]$
${\epsilon }_{\tau \tau }^p$	$[-2.209,5.710]$	$[-0.098,\mathrm{\hspace{0.17em}0.221}]\oplus [1.063,\mathrm{\hspace{0.17em}1.408}]$
${\epsilon }_{e\mu }^p$	$[-2.209,2.249]$	$[-0.124,\mathrm{\hspace{0.17em}0.058}]$
${\epsilon }_{e\tau }^p$	$[-2.434,2.006]$	$[-0.239,\mathrm{\hspace{0.17em}0.244}]$
${\epsilon }_{\mu \tau }^p$	$[-3.849,1.772]$	$[-0.013,\mathrm{\hspace{0.17em}0.021}]$
${\epsilon }_{ee}^n$	$[-1.714,2.915]$	—
${\epsilon }_{\mu \mu }^n$	$[-2.331,3.282]$	—
${\epsilon }_{\tau \tau }^n$	$[-1.564,2.705]$	—
${\epsilon }_{e\mu }^n$	$[-1.426,1.846]$	—
${\epsilon }_{e\tau }^n$	$[-1.829,1.147]$	—
${\epsilon }_{\mu \tau }^n$	$[-2.275,1.250]$	—

Table 3. Table 3: Limits on electron NSI, most extracted from Ref.  Farzan:2017xzy . The limits derived from our analysis of the recent XENONnT ER results are highlighted in orange. The limits derived by Ref.  Coloma:2022umy using the Borexino Phase-II data, which we compare to in Fig.   5 , are highlighted in green.

	L	R	V	Ref.
${\epsilon }_{ee}^e$	$[-0.021,0.052]$	$[-0.18,0.51]$	–	SK & KamLAND Bolanos:2008km
	$[-0.046,0.053]$	$[-0.21,0.16]$	–	Borexino Phase-I Agarwalla:2012wf
	–	–	$[-0.56,0.24]$	Borexino & COHERENT Dutta:2020che
	$[-1.37,-1.29]\oplus$ $[0.03,0.06]$	$[-0.23,0.07]$	$[-0.09,0.14]$	Borexino Phase-II Coloma:2022umy
	–	–	$[-2.65,0.78]$	XENONnT 2022 (this work)
${\epsilon }_{\mu \mu }^e$	–	–		SK atm. Gonzalez-Garcia:2011vlg
	$[-0.03,0.03]$	$[-0.03,0.03]$	–	react. + acc. Davidson:2003ha ; Barranco:2007ej
	–	–	$[-0.58,0.72]$	Borexino & COHERENT Dutta:2020che
	$[-0.20,0.13]\oplus$ $[0.58,0.81]$	$[-0.36,0.37]$	$[-0.51,0.35]$	Borexino Phase-II Coloma:2022umy
	–	–	$[-2.19,2.34]$	XENONnT 2022 (this work)
${\epsilon }_{\tau \tau }^e$	$[-0.12,0.060]$	$[-0.99,0.23]$	–	SK & KamLAND Bolanos:2008km
	$[-0.23,0.87]$	$[-0.98,0.73]$	–	Borexino Phase-I Agarwalla:2012wf
	–	–		SK atm. Gonzalez-Garcia:2011vlg
	–	–	$[-0.60,0.72]$	Borexino & COHERENT Dutta:2020che
	$[-0.26,0.26]\oplus$ $[0.45,0.86]$	$[-0.58,0.47]$	$[-0.66,0.52]$	Borexino Phase-II Coloma:2022umy
	–	–	$[-2.09,2.20]$	XENONnT 2022 (this work)
${\epsilon }_{e\mu }^e$	$[-0.13,0.13]$	$[-0.13,0.13]$	–	react. + acc. Barranco:2007ej
	–	–	$[-0.58,0.60]$	Borexino & COHERENT Dutta:2020che
	$[-0.17,0.29]$	$[-0.21,0.41]$	$[-0.34,0.61]$	Borexino Phase-II Coloma:2022umy
	–	–	$[-1.03,1.41]$	XENONnT 2022 (this work)
${\epsilon }_{e\tau }^e$	$[-0.33,0.33]$	$[-0.28,-0.05]\oplus$ $[0.05,0.28]$	–	react. + acc. Barranco:2007ej
	–	$[-0.19,0.19]$	–	TEXONO TEXONO:2010tnr
	–	–	$[-0.60,0.62]$	Borexino & COHERENT Dutta:2020che
	$[-0.26,0.23]$	$[-0.35,0.31]$	$[-0.48,0.47]$	Borexino Phase-II Coloma:2022umy
	–	–	$[-1.26,1.11]$	XENONnT 2022 (this work)
${\epsilon }_{\mu \tau }^e$	–	–	$[-0.035,0.018]$	SK atm. Gonzalez-Garcia:2011vlg
	–	–	$[-0.20,0.07]$	MINOS MINOS:2013hmj
	–	–	$[-0.018,0.016]$	IceCube Salvado:2016uqu
	–	–	$[-0.67,0.62]$	Borexino & COHERENT Dutta:2020che
	$[-0.62,-0.52]\oplus$ $[-0.09,0.14]$	$[-0.26,0.23]$	$[-0.25,0.36]$	Borexino Phase-II Coloma:2022umy
	–	–	$[-1.57,1.50]$	XENONnT 2022 (this work)

Equations

$$\text{diff}R{E}_R=N_T{\int }_{E_{\nu }^{\min }}\text{diff}{\varphi }_{\nu }{E}_{\nu }\mathrm{Tr}\left[\mathbf{\rho }\text{diff}\mathbf{\zeta }{E}_R\right]\text{dl}{E}_{\nu },$$

$${\mathcal{L}}_{\text{NSI}}=-2\sqrt{2}{G}_F\sum _{\begin{array}{c}f=e,u,d\\ \alpha ,\beta =e,\mu ,\tau \end{array}}{\epsilon }_{\alpha \beta }^{fP}\left[{\bar{\nu }}_{\alpha }{\gamma }_{\rho }{P}_L{\nu }_{\beta }\right]\left[\bar{f}{\gamma }^{\rho }Pf\right],$$

$${\epsilon }_{\alpha \beta }^{fP}={\epsilon }_{\alpha \beta }^{\eta ,\phi }{\xi }^{fP},$$

$${\epsilon }_{\alpha \beta }^f$$

$${\tilde{\epsilon }}_{\alpha \beta }^f$$

$${\epsilon }_{\alpha \beta }^p=2{\epsilon }_{\alpha \beta }^u+{\epsilon }_{\alpha \beta }^d,$$

$${\epsilon }_{\alpha \beta }^n={\epsilon }_{\alpha \beta }^u+2{\epsilon }_{\alpha \beta }^d.$$

$${\xi }^e$$

$${\xi }^p$$

$${\xi }^n$$

$$H^{\nu }$$

$$H^{\bar{\nu }}$$

$$H_{\mathrm{vac}}=U_{\mathrm{PMNS}}\frac{1}{2E_{\nu }}\left(\begin{array}{ccc}0& 0& 0\\ 0& \Delta m_{21}^2& 0\\ 0& 0& \Delta m_{31}^2\end{array}\right)U_{\mathrm{PMNS}}^{†},$$

$$U_{\mathrm{PMNS}}=\underset{\equiv R_{23}}{\underbrace{\left(\begin{array}{ccccc}1& & 0& & 0\\ 0& & c_{23}& & s_{23}\\ 0& & -s_{23}& & c_{23}\end{array}\right)}}\underset{\equiv R_{13}}{\underbrace{\left(\begin{array}{ccccc}{c}_{13}& & 0& & s_{13}\\ 0& & 1& & 0\\ -s_{13}& & 0& & c_{13}\end{array}\right)}}\underset{\equiv U_{12}}{\underbrace{\left(\begin{array}{ccccc}{c}_{12}& & s_{12}{e}^{i{\delta }_{\mathrm{CP}}}& & 0\\ -s_{12}{e}^{-i{\delta }_{\mathrm{CP}}}& & c_{12}& & 0\\ 0& & 0& & 1\end{array}\right)}}.$$

$$H_{\mathrm{mat}}=\sqrt{2}{G}_F{N}_e(x)\left(\begin{array}{ccc}1+{\mathcal{E}}_{ee}(x)& {\mathcal{E}}_{e\mu }(x)& {\mathcal{E}}_{e\tau }(x)\\ {\mathcal{E}}_{e\mu }^{\ast }(x)& {\mathcal{E}}_{\mu \mu }(x)& {\mathcal{E}}_{\mu \tau }(x)\\ {\mathcal{E}}_{e\tau }^{\ast }(x)& {\mathcal{E}}_{\mu \tau }^{\ast }(x)& {\mathcal{E}}_{\tau \tau }(x)\end{array}\right),$$

$${\mathcal{E}}_{\alpha \beta }=\sum _f\frac{N_f(x)}{N_e(x)}{\epsilon }_{\alpha \beta }^f,$$

$${\mathcal{E}}_{\alpha \beta }={\epsilon }_{\alpha \beta }^e+{\epsilon }_{\alpha \beta }^p+Y_n(x){\epsilon }_{\alpha \beta }^n=\left[{\xi }^e+{\xi }^p+Y_n(x){\xi }^n\right]{\epsilon }_{\alpha \beta }^{\eta ,\phi },$$

$$H^{\nu }=\frac{1}{2E_{\nu }}\left(\begin{array}{ccc}{c}_{13}^2{A}_{\text{cc}}+s_{12}^2\Delta m_{21}^2& s_{12}{c}_{12}{e}^{i{\delta }_{\mathrm{CP}}}\Delta m_{21}^2& s_{13}{c}_{13}{A}_{\text{cc}}\\ s_{12}{c}_{12}{e}^{-i{\delta }_{\mathrm{CP}}}\Delta m_{21}^2& c_{12}^2\Delta m_{21}^2& 0\\ s_{13}{c}_{13}{A}_{\text{cc}}& 0& s_{13}^2{A}_{\text{cc}}+\Delta m_{31}^2\end{array}\right)+\sqrt{2}{G}_F{N}_e(x)O^{†}\mathcal{E}O,$$

$${\epsilon }_{\alpha \beta }^{\eta ,\phi }\ll \frac{\sqrt{2}\Delta m_{31}^2}{A_{\mathrm{cc}}(x)\left[{\xi }^e+{\xi }^p+Y_n(x){\xi }^n\right]}.$$

$$H_{\mathrm{vac}}^{\mathrm{eff}}\equiv \frac{\Delta m_{21}^2}{4E_{\nu }}\left(\begin{array}{cc}-\cos 2{\theta }_{12}& \sin 2{\theta }_{12}\\ \sin 2{\theta }_{12}& \cos 2{\theta }_{12}\end{array}\right),$$

$$H_{\mathrm{mat}}^{\mathrm{eff}}\equiv \sqrt{2}{G}_F{N}_e(x)\left[\left(\begin{array}{cc}{c}_{13}^2& 0\\ 0& 0\end{array}\right)+\left[{\xi }^e+{\xi }^p+Y_n(x){\xi }^n\right]\left(\begin{array}{cc}-{\epsilon }_D^{\eta ,\phi }& {\epsilon }_N^{\eta ,\phi }\\ {\epsilon }_N^{\eta ,\phi }& {\epsilon }_D^{\eta ,\phi }\end{array}\right)\right].$$

$${\epsilon }_D^{\eta ,\phi }\equiv$$

$$-\frac{c_{13}^2}{2}\left({\epsilon }_{ee}^{\eta ,\phi }-{\epsilon }_{\mu \mu }^{\eta ,\phi }\right)+\frac{s_{23}^2-s_{13}^2{c}_{23}^2}{2}\left({\epsilon }_{\tau \tau }^{\eta ,\phi }-{\epsilon }_{\mu \mu }^{\eta ,\phi }\right),$$

$${\epsilon }_N^{\eta ,\phi }\equiv c_{13}\left(c_{23}{\epsilon }_{e\mu }^{\eta ,\phi }-s_{23}{\epsilon }_{e\tau }^{\eta ,\phi }\right)+s_{13}\left[s_{23}^2{\epsilon }_{\mu \tau }^{\eta ,\phi }-c_{23}^2{\epsilon }_{\mu \tau }^{\eta ,\phi }+c_{23}{s}_{23}\left({\epsilon }_{\tau \tau }^{\eta ,\phi }-{\epsilon }_{\mu \mu }^{\eta ,\phi }\right)\right].$$

$$U_{12}^m=\left(\begin{array}{cc}\cos {\theta }_{12}^m& \sin {\theta }_{12}^m\\ -\sin {\theta }_{12}^m& \cos {\theta }_{12}^m\end{array}\right),$$

$$E_{1,2}^m$$

$$p\equiv$$

$$q\equiv$$

$$\Delta E_{21}^m\equiv E_2^m-E_1^m=\frac{\Delta m_{21}^2}{2E_{\nu }}\sqrt{p^2+q^2}.$$

$$\sin 2{\theta }_{12}^m$$

$$\cos 2{\theta }_{12}^m$$

$$\tan 2{\theta }_{12}^m$$

$$i\text{diff}x\left(\begin{array}{c}{\hat{\nu }}_e\\ {\hat{\nu }}_{\mu }\\ {\hat{\nu }}_{\tau }\end{array}\right)=\underset{\equiv \tilde{S}}{\underbrace{\left(\begin{array}{ccc}\mathrm{Evol}[H^{\mathrm{eff}}]& & 0\\ 0& & \exp [-i\frac{\Delta m_{31}^2}{2E_{\nu }}L]\end{array}\right)}}\left(\begin{array}{c}{\hat{\nu }}_e\\ {\hat{\nu }}_{\mu }\\ {\hat{\nu }}_{\tau }\end{array}\right).$$

$$\mathrm{Evol}[H^{\mathrm{eff}}]=\underset{\Delta x\to 0}{\lim }\prod _{n=0}^N{U}_{\mathrm{PMNS}}^m(x_n)\exp \left[-i\left(U_{12}^{m†}{H}^{\mathrm{eff}}{U}_{12}^m-i{U}_{12}^{m†}{\dot{U}}_{12}^m\right)\Delta x\right]U_{\mathrm{PMNS}}^m(x_n{)}^{†},$$

$$i\text{diff}x\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\\ {\nu }_{\tau }\end{array}\right)=\underset{S}{\underbrace{O\tilde{S}{O}^{†}}}\left(\begin{array}{c}{\nu }_e\\ {\nu }_{\mu }\\ {\nu }_{\tau }\end{array}\right),$$

Full text

††preprint: IPPP/23/08††preprint: IFT-UAM/CSIC-23-19††preprint: FT-UAM-23-1

1. I Introduction

2. II Solar Neutrino Physics and Non-Standard Interactions

3. II.1 NSI parametrisation

4. II.2 Three-flavour neutrino oscillations in the presence of NSI

5. II.3 Solar neutrino density matrix

6. II.4 Generalised neutrino cross sections

7. II.4.1 CE$\nu$NS cross section

8. II.4.2 E$\nu$ES cross section

9. III Extending current constraints to the full NSI parameter space

10. III.1 The CENNS-10 LAr Experiment

11. III.2 The Borexino Experiment

12. IV Direct Detection Experiments

13. IV.1 Expected Number of Events and Statistical Procedure

14. IV.2 Sensitivities in the nucleon NSI plane

15. IV.3 Sensitivities in the charged NSI plane

16. IV.4 Final Remarks

17. V Conclusions

18. A Solar neutrino transition rate

19. B $\Delta\chi^{2}$ Plots for CENNS-10 LAr and Borexino

I Introduction

Neutrinos are among the most mysterious particles of the Standard Model (SM). The discovery of their flavour oscillations remains one of the strongest pieces of evidence for new physics, since it requires neutrinos to be massive Davis:1968cp ; SNO:2001kpb . In the SM, however, neutrinos are described by a purely left-handed spinorial field forming part of an $SU(2)_{L}$ doublet Weinberg:1967tq , which, by the principles of gauge invariance, disallows a mass term for neutrinos at the renormalisable level. Hence, the neutrino sector has inspired a myriad of extensions beyond the SM (BSM) (see, for example, Ref. deGouvea:2016qpx for a review).

A convenient parametrisation of new physics in neutrino interactions has been established in terms of the low-energy effective field theory (EFT) of non-standard interactions (NSI) Wolfenstein:1977ue ; Guzzo:1991hi ; Guzzo:1991cp ; Gonzalez-Garcia:1998ryc ; Bergmann:2000gp ; Guzzo:2000kx ; Guzzo:2001mi ; GonzalezGarcia:2004wg . This formalism contemplates modifications to neutrino interactions with SM particles while respecting the SM vector current structure. Over the last decades, a variety of experimental bounds have been derived for the NSI couplings Fornengo:2001pm ; Huber:2002bi ; Davidson:2003ha ; Maltoni:2004ei ; Kopp:2007mi ; Kopp:2007ne ; Kopp:2008ds ; Antusch:2008tz ; Bolanos:2008km ; Biggio:2009kv ; Biggio:2009nt ; Escrihuela:2009up ; Kopp:2010qt ; Escrihuela:2011cf ; Agarwalla:2014bsa ; Gonzalez-Garcia:2015qrr ; Liao:2016orc ; Liao:2016reh ; Farzan:2017xzy ; Coloma:2017ncl ; Liao:2017awz ; Coloma:2017egw ; Liao:2017uzy ; Esteban:2018ppq ; Denton:2018xmq ; AristizabalSierra:2018eqm ; Coloma:2019mbs ; Han:2019zkz ; Esteban:2019lfo ; Proceedings:2019qno ; Miranda:2020tif ; Denton:2020hop ; Martinez-Mirave:2021cvh ; Masud:2021ves ; Chatterjee:2021wac ; Chaves:2021pey ; deSalas:2021aeh ; IceCubeCollaboration:2021euf ; Denton:2022nol ; Denton:2022pxt ; DeRomeri:2022twg ; Brahma:2023wlf , most importantly from neutrino oscillation and spallation source experiments. The latter have recently succeeded in observing the coherent elastic scattering of neutrinos with nuclei (CE$\nu$NS) COHERENT:2015mry ; COHERENT:2017ipa with a rate consistent with the SM prediction Freedman:1973yd ; Drukier:1984vhf , providing stronger constraints on new physics contributions (see, for example, Refs. Kosmas:2017zbh ; Miranda:2020syh ; Abdullah:2018ykz ; Bauer:2018onh ; Miranda:2020zji ; Bauer:2020itv ; Kosmas:2015sqa ).

Meanwhile, dark matter (DM) direct detection (DD) experiments have experienced remarkable progress. Current detectors have significantly increased their target size and sensitivity, to the point where they will be able to observe the scattering of solar neutrinos. This constitutes a new background for DM searches, which leads to the so-called neutrino floor (or fog) Billard:2013qya ; Billard:2021uyg ; OHare:2021utq , but it also offers the unique opportunity to probe new physics with these instruments Cerdeno:2016sfi ; Dutta:2017nht ; Gelmini:2018gqa ; Essig:2018tss ; Amaral:2020tga ; Dutta:2020che ; Amaral:2021rzw ; Munoz:2021sad ; deGouvea:2021ymm . Neutrinos can be observed in DD experiment through elastic neutrino-electron scatting (E$\nu$ES) or their coherent elastic scattering with nuclei. Due to their larger target size, liquid noble gas detectors like LUX-ZEPLIN (LZ) Aalbers:2022fxq , PandaX PandaX-4T:2021bab , and XENONnT XENONCollaboration:2022kmb , are better positioned than other DD techniques to carry out this kind of search. Indeed, the larger xenon-based DD experiments have thus far succeeded in placing upper bounds on the $\mathrm{{}^{8}B}$ neutrino flux XENON:2020rca ; PandaX:2022aac . The sensitivity of DD experiments to these processes has already led to studies in which the expected solar neutrino scattering rate has been used as a laboratory for gaining a deeper understanding of the nature of solar physics, neutrino oscillations, and BSM neutrino physics Harnik:2012ni ; Baudis:2013qla ; Billard:2014yka ; Cerdeno:2016sfi ; Cui:2017ytb ; Dutta:2017nht ; Cerdeno:2017xxl ; Newstead:2018muu ; Essig:2018tss ; AristizabalSierra:2019ykk ; Dutta:2019oaj ; Dutta:2020che ; deGouvea:2021ymm ; A:2022acy .

In this work, we set out to exploit the sensitivity of DD experiments to solar neutrino scattering with the aim of exploring their impact on the NSI landscape. In the context of dedicated neutrino experiments, NSI studies are numerous; however, the potential of DD experiments has not been fully investigated. Previous works have pointed out that non-zero NSI parameters could produce appreciable signals for both CE$\nu$NS and E$\nu$ES Dutta:2017nht ; Dutta:2020che in DD, as well as potentially modify the neutrino fog AristizabalSierra:2017joc ; Gonzalez-Garcia:2018dep .

To this end, we will introduce a convenient parametrisation of NSI, extending the framework of Ref. Esteban:2018ppq to include an explicit separation between NSI in the electron and proton directions. This is needed to interpret the results of DD experiments. Ignoring the electron contribution is a valid choice as long as one is mostly interested in matter effects for oscillation experiments111Non-standard matter effects enter the matter Hamiltonian via a contribution from the neutron and an overall charged contribution from both the proton and the electron., but this is a non-general treatment. Striving for greater generality, we allow for the possibility that the ‘charged’ neutrino NSI is shared between both the proton and the electron. While the total charged contribution can be designed to leave neutrino oscillations unchanged, allowing for electron NSI can instead lead to changes in the E$\nu$ES cross section. This, in turn, can affect the bounds set by oscillation experiments, which could instead be dominated by NSI effects at the detection point Esteban:2018azc .

Furthermore, as was recently pointed out in Ref. Coloma:2022umy , when new physics introduces potential flavour-changing neutral currents (FCNC), the full flavour structure of the cross section must be retained when dealing with a neutrino flux composed of an admixture of flavour eigenstates. This is in contrast with the SM, where interactions are diagonal in the flavour basis. Thus, in the general NSI case, it is no longer appropriate to project the neutrino state that reaches Earth onto any one particular flavour state and convolve the result with flavour-pure cross sections, as neutrinos arrive in a superposition of flavour eigenstates.222We stress that the simplified treatment of calculating the number of neutrino scattering events in the presence of new physics, given by $N_{\nu}\propto\sum_{\alpha}P_{e\alpha}\,{\mathrm{d}\sigma_{\nu_{\alpha}T}}/{\mathrm{d}E_{R}}$, where $P_{e\alpha}$ is the transition probability to a neutrino of flavour $\alpha$, is only appropriate in two cases. Firstly, if the flux of neutrinos incident on a target is only composed of one flavour. Secondly, if the new physics contribution is flavour-conserving. Instead, we must consider the full flavour-structure of both the cross section and the density matrix describing the evolution of the initial neutrino state.

The rate of neutrino events in a generic neutrino scattering experiment is then described by the expression Coloma:2022umy ,

[TABLE]

where $N_{T}$ is the number of targets, $\phi_{\nu}$ is the neutrino flux at the source333Of particular relevance to experiments sensitive to solar neutrinos is the fact that electron neutrino production in the Sun proceeds through a series of charged-current interactions. Since we are only considering neutral-current NSI, the electron neutrino flux produced in the Sun is unchanged., $\bm{\rho}$ is the neutrino density matrix at the experiment and $\bm{\zeta}$ is a generalised scattering cross section in the neutrino-flavour space, encoding correlations between scattering amplitudes of neutrinos with different flavours. Here, $E_{\nu}$ is the energy of the incident neutrinos and $E_{\nu}^{\mathrm{min}}$ is the minimum $E_{\nu}$ required to produce a target recoil energy of $E_{R}$.

Using this generalised framework, in this paper we study how DD experiments will constitute a valuable complementary probe of the NSI landscape. To do this, we first explore how previous limits derived from oscillation and spallation source experiments map onto the full NSI parameter space. Then, we derive new limits from the recent data from LZ and XENONnT onto the NSI parameters and make projections for their full exposure runs, as well as for the future DARWIN experiment DARWIN:2016hyl . We demonstrate that xenon-based DD experiments like XENON XENON:2018voc ; XENON:2020kmp ; XENON:2020rca , DARWIN DARWIN:2016hyl , PandaX PandaX-II:2017hlx ; PandaX:2018wtu ; PandaX-4T:2021bab and LZ Mount:2017qzi ; LUX-ZEPLIN:2018poe ; LZ:2021xov will be sensitive to generic NSI couplings in a competitive and complementary way to oscillation, beam, and spallation source experiments. We do this by comparing our results and projections to those derived in Refs. Farzan:2017xzy ; Esteban:2018ppq ; Coloma:2019mbs ; Dutta:2020che ; Coloma:2022umy . Given that DD will be sensitive to both CE$\nu$NS and E$\nu$ES, we explore their limits and projections in our extended parametrisation, emphasising the complementarity of both signals. Indeed, due to the high flux of solar neutrinos and the excellent background reduction of DD experiments, their sensitivity to electron NSI is remarkable and is competitive with that of conventional neutrino oscillation experiments.

This article is organised as follows. In Section II, we introduce the framework of non-standard neutrino interactions, explicitly incorporating interactions with electrons as well as the impact of such NSI on solar neutrino physics. We then derive the relevant formalism for computing the density matrix and the generalised cross section, both required to compute the expected solar neutrino scattering rates. In Section III, we shed light on the current landscape of NSI constraints derived from oscillation and spallation source experiments, as well as their sensitivity to interactions with electrons. We present and discuss the results of our sensitivity study of DD experiments to NSI as the main results of this work in Section IV. Finally, we draw our conclusions in Section V.

II Solar Neutrino Physics and Non-Standard Interactions

In this section, we introduce the framework of neutrino NSI and study their impact on solar neutrino physics, both in propagation effects and in scattering with nuclei and electrons. In doing so, we derive the relevant expressions for the solar neutrino density matrix, $\bm{\rho}$, and the generalised cross section, $\bm{\zeta}$, for CE$\nu$NS and E$\nu$ES entering the rate in Eq. 1. For an explanation of the origin of this rate equation, we refer to Appendix A.

II.1 NSI parametrisation

In order to understand how potential new neutrino interactions enter the scattering rate in Eq. 1, we need to specify a BSM model. Since we want to remain as general as possible about the origin of such new physics, we will work in terms of a low-energy effective theory. Making use of the framework of neutrino NSI Wolfenstein:1977ue ; Miranda:2015dra ; Davidson:2003ha ; GonzalezGarcia:2004wg , we can parametrise new physics effects in the neutrino sector by contact terms of the form444Note that this parametrisation is not $SU(2)_{L}$ invariant and is mainly motivated by the structure of the SM weak current. In order to systematically capture all gauge invariant dimension-six operators modifying neutrino interactions, it is more suitable to consider a complete basis of EFT operators and map them onto the enlarged basis of general neutrino interactions Gavela:2008ra ; Lindner:2016wff ; Bischer:2019ttk ; Terol-Calvo:2019vck .

[TABLE]

where $G_{F}$ denotes the Fermi constant and $P\in\{P_{L},P_{R}\}$. The NSI parameters $\varepsilon^{fP}_{\alpha\beta}$, which are in general flavour-violating, quantify the strength of the interaction between the neutrinos $\nu_{\alpha}$ and $\nu_{\beta}$ and the pair of fermions $f$ relative to the SM weak interaction, characterised by $G_{F}$. In this work, we will not consider any new source of CP-violation and hence assume the parameters $\varepsilon^{fP}_{\alpha\beta}$ to be real. Furthermore, we have assumed that the charged fermions $f$ are identical, resembling the SM neutral current (NC) interaction. However, charged current (CC) NSI could also exist, where the neutrinos couple to two different charged fermions, $f$ and $f^{\prime}$. Since these are, in general, subject to much harsher constraints and DD experiments do not probe CC interactions in NRs, we do not consider them here but rather direct the reader to, for example, Refs. Kopp:2010qt ; Akhmedov:2010vy ; Ohlsson:2012kf ; Miranda:2015dra ; Farzan:2017xzy ; Khan:2021wzy ; Falkowski:2021bkq .

To describe neutrinos interacting with ordinary matter (made up of electrons, protons and neutrons), only interactions with the first generation of SM fermions need to be considered. If we assume that the neutrino flavour structure of the NSI is independent of the charged fermion $f$ that the neutrinos couple to, we can factorise the NSI coupling as Esteban:2018ppq

[TABLE]

where $\xi^{fP}$ describes the relative strength of the interaction with the fermions $f\in\{e,u,d\}$ and $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ denotes the overall strength of the NSI. We further define the vector and axial-vector NSI couplings

[TABLE]

with $\xi^{f}=\xi^{fL}+\xi^{fR}$ and $\tilde{\xi}^{f}=\xi^{fL}-\xi^{fR}$. As matter effects are only sensitive to the vector part of the interaction, we focus only on vector NSI in this work, setting $\tilde{\varepsilon}_{\alpha\beta}$ to zero. Since we are ultimately testing neutrino interactions with matter, it is convenient to parametrise the NSI with quarks in terms of proton and neutron NSI,

[TABLE]

Extending the parametrisation of Ref. Esteban:2018ppq by re-introducing the electron direction via a second angle $\varphi$, the relative strengths of the electron, proton, and neutron NSI are written as555The normalisation factor of $\sqrt{5}$ was originally introduced in Ref. Esteban:2018ppq to have unit vectors $\xi^{u}$ and $\xi^{d}$ if the NSI are entirely in the up- and down-quark direction, respectively. We cohere to this normalisation for comparability of our results with the literature on NSI global fits.,

[TABLE]

In Fig. 1, we illustrate our parametrisation of the three base NSI directions $\hat{\varepsilon}^{p}_{\alpha\beta},\hat{\varepsilon}^{n}_{\alpha\beta}$ and $\hat{\varepsilon}^{e}_{\alpha\beta}$. We define the angle $\eta$ as the angle of general NSI coupling $\bm{\varepsilon_{\alpha\beta}}=(\varepsilon^{p}_{\alpha\beta},\varepsilon^{n}_{\alpha\beta},\varepsilon^{e}_{\alpha\beta})$ with the plane of charged NSI $(\hat{\varepsilon}^{p}_{\alpha\beta},\hat{\varepsilon}^{e}_{\alpha\beta})$. The second angle $\varphi$ is defined as the angle between the general NSI element $\bm{\varepsilon_{\alpha\beta}}$ and the plane of hadronic NSI $(\hat{\varepsilon}^{p}_{\alpha\beta},\hat{\varepsilon}^{n}_{\alpha\beta})$. In order to match our notation with the literature on global NSI fits Esteban:2018ppq ; Coloma:2019mbs , we allow for both positive and negative values for $\varepsilon^{\eta,\varphi}_{\alpha\beta}$. Thus, the azimuthal angle $\eta$ only runs in the interval $[-\pi/2,\pi/2]$ to span the full two-dimensional plane of hadronic NSI $(\varepsilon^{p}_{\alpha\beta},\varepsilon^{n}_{\alpha\beta})$. The second, polar angle $\varphi$ (taken from the hadronic NSI plane) also runs in the interval $[-\pi/2,\pi/2]$ to cover the full sphere. For example, $\eta=0$ and $\varphi=0$ corresponds to NSI only in the proton direction $\hat{\varepsilon}^{p}_{\alpha\beta}$, $\eta=0$ and $\varphi=\pi/2$ to NSI only in the electron direction $\hat{\varepsilon}^{e}_{\alpha\beta}$, and $\eta=\pi/2$ to NSI only in the neutron direction $\hat{\varepsilon}^{n}_{\alpha\beta}$.

II.2 Three-flavour neutrino oscillations in the presence of NSI

With this extended framework, we can describe the evolution of neutrino and antineutrino states during propagation in the Hamiltonian formalism by

[TABLE]

where the standard vacuum Hamiltonian is given by

[TABLE]

with $\Delta m_{ij}^{2}\equiv m_{i}^{2}-m_{j}^{2}$ and $U_{\mathrm{PMNS}}$ being the PMNS matrix, defined as

[TABLE]

Here, $\delta_{\mathrm{CP}}$ is the CP-phase, and $c_{ij}$ and $s_{ij}$ refer to $\cos\theta_{ij}$ and $\sin\theta_{ij}$, respectively.

The matter Hamiltonian, consisting of both the SM charged current and the NSI neutral current contributions, is given by

[TABLE]

with

[TABLE]

where $N_{f}(x)$ is the spatial fermion density in matter. With the definition of the nuclear NSI couplings in Eq. 5 and the fact that in neutral matter $N_{p}(x)=N_{e}(x)$, we can express the dimensionless NSI matter Hamiltonian elements as

[TABLE]

where $Y_{n}(x)={N_{n}(x)}/{N_{e}(x)}$ denotes the fractional neutron density. In studying solar neutrino propagation effects, we take $Y_{n}(x)$ from Ref. Bahcall:2005va .

In the context of solar neutrino physics, it is convenient to switch from the conventional neutrino flavour basis to a new basis $\bm{\hat{\nu}}=O^{\dagger}\bm{\nu}$, which we will refer to as the solar neutrino flavour basis, via the rotation $O=R_{23}\,R_{13}$. In this basis, the full Hamiltonian reads,

[TABLE]

where we have defined the matter potential $A_{\text{cc}}=2\,E_{\nu}\,V_{\text{cc}}=2\,E_{\nu}\sqrt{2}G_{F}N_{e}(x)$. From the structure of the Hamiltonian above, we see that if $\Delta m_{31}^{2}\gg\Delta m_{21}^{2},\,A_{\mathrm{cc}},\,2E_{\nu}\,G_{F}\sum_{f}N_{f}(x)\varepsilon_{\alpha\beta}^{f}$, the Hamiltonian is dominated by the third eigenvalue, $\Delta m_{31}^{2}$. In this case, it is effectively block-diagonal, turning our $3\nu$ problem into a $2\nu$ one. In this rotated basis, the third mass eigenstate decouples from the rest of the system and evolves adiabatically. Throughout its journey from the Sun to the Earth, this third eigenstate can be well-approximated by its vacuum mass eigenstate. Within this approximation, the Hamiltonian in Eq. 13 is transformed to an effective $2\times 2$ picture, where we only have to track the evolution of the two lighter matter mass eigenstates.

The first condition, $\Delta m_{31}^{2}\gg\Delta m_{21}^{2}$, is satisfied by current best-fits to oscillation parameters Esteban:2020cvm . The second is satisfied for solar neutrinos across the full range of solar neutrino energies, $E_{\nu}\lesssim$20\text{,}\mathrm{MeV}$$. The third condition can be interpreted as one on the value of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$,

[TABLE]

Taking the maximum of all these quantities, which occurs at the solar core, and using $E_{\nu}\sim$20\text{,}\mathrm{MeV}$$, we find that $\varepsilon^{\eta,\varphi}_{\alpha\beta}\lesssim 3$. We treat this as an upper bound on the value of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$, and we do not interpret our results above this value throughout our analyses. Ultimately, for $\varepsilon^{\eta,\varphi}_{\alpha\beta}\sim 3$ at these higher neutrino energies, which are relevant for NRs in DD experiments, a full numerical simulation should be performed to more accurately model neutrino oscillations. For the purposes of our sensitivity study, however, our approach is sufficient.

Following the conventions of Ref. Esteban:2018ppq and setting $\delta_{\mathrm{CP}}=0$ in this work, we can write the effective Hamiltonian as $H^{\mathrm{eff}}\equiv H_{\mathrm{vac}}^{\mathrm{eff}}+H_{\mathrm{mat}}^{\mathrm{eff}}$, where

[TABLE]

and

[TABLE]

The coefficients $\varepsilon_{N}^{\eta,\varphi}$ and $\varepsilon_{D}^{\eta,\varphi}$ are related to our parametrisation by

[TABLE]

and

[TABLE]

Diagonalising $H^{\mathrm{eff}}$ then allows us to find the matrix $U^{m}_{12}$ such that $U^{m\dagger}_{12}H^{\mathrm{eff}}U^{m}_{12}=\operatorname{diag}(E^{m}_{1},E^{m}_{2})$. Typically, $U^{m}_{12}$ is parametrised as

[TABLE]

for some matter mixing angle $\theta_{12}^{m}$. We find that the eigenvalues of the effective matter Hamiltonian $H^{\mathrm{eff}}$ are given by

[TABLE]

where we have defined the two quantities

[TABLE]

Thus, the energy difference between the two energy eigenvalues in matter, responsible for the coherent mixing of the two matter mass eigenstates, is given by

[TABLE]

Moreover, we find that the matter mixing angle, $\theta^{m}_{12}$, obeys the relations

[TABLE]

With these expressions, we are in the position to describe the neutrino evolution in the full $3\times 3$ picture. Using the notation of Ref. Coloma:2022umy and the fact that solar neutrinos are relativistic (such that $t\simeq x$), we can write the evolution equation in the solar neutrino flavour basis as

[TABLE]

To obtain the evolved matter Hamiltonian, we can split up the distance of propagation within the Sun into $N$ equidistant slabs of thickness $\Delta x$ with approximate homogeneous matter density. We can then obtain the evolved Hamiltonian by formally taking the limit

[TABLE]

where $\dot{U}^{m}_{12}=(\mathrm{d}/{\mathrm{d}x})\,U^{m}_{12}$ and $x_{n}=x_{n-1}+\Delta x$. From this, we can write the full evolution equation in the conventional vacuum-flavour basis as

[TABLE]

with the rotation matrix $O=R_{23}R_{13}$. In this notation, the full $S$-matrix is thus given by

[TABLE]

with $\Phi_{33}=\Delta m^{2}_{31}L/(2E_{\nu})$, where we evolve the neutrinos from their production point within the Sun, $x_{0}$, to their detection point at an experiment over the distance $L$. Finally, the $2\times 2$ time-evolution matrix is given by

[TABLE]

To simplify our analysis, we make the assumption that the two light matter mass eigenstates of $H^{\mathrm{eff}}$, $\ket{\nu_{1m}}$ and $\ket{\nu_{2m}}$, propagate adiabatically within the Sun. As such, the two eigenstates do not mix with one another as they travel to the surface of the Sun, remaining eigenstates of $H^{\mathrm{eff}}$ throughout their evolution. This assumption is appropriate because the matter density within the Sun, described by $N_{f}(x)$, varies slowly enough to allow the matter eigenstates to adapt to the medium as they propagate through it. The adiabatic approximation is valid if the adiabaticity parameter, $\gamma$, satisfies

[TABLE]

where $\Delta E_{21}^{m}$ is given by Eq. 22. In the adiabatic approximation, the matrix $D(x)$ is thus approximately diagonal, and after a common rephasing of the neutrino matter mass eigenstates, the upper $2\times 2$ block in Eq. 27 can be expressed as

[TABLE]

with $\phi=\int_{0}^{L}\Delta E_{21}^{m}(x)\,\dl x$. Since the neutrinos exiting the Sun will free-stream to the Earth, there is no further evolution effect to be taken into account for the Sun-Earth propagation.

However, in principle, there is a further propagation effect when neutrinos pass thorough the Earth at night, which should be taken into account for a complete treatment. For high-energy $\mathrm{{}^{8}B}$ neutrinos, for which $E_{\nu}\sim$10\text{,}\mathrm{MeV}$$, this effect typically changes oscillation probabilities only at the percent level Bahcall:1997jc ; Akhmedov:2000cs ; Chiang:2000kh . In particular, Super-Kamiokande has determined the day-night asymmetry to about $-3.3\%$ in $\mathrm{{}^{8}B}$ neutrinos Super-Kamiokande:2013mie , while Borexino has found no asymmetry in $\mathrm{{}^{7}Be}$ neutrinos Borexino:2011bhn . Therefore, in this work we neglect Earth matter effects for simplicity.

II.3 Solar neutrino density matrix

From the expression of the $S$-matrix in Eq. 27, we can derive the expression for the full three-flavour density matrix for solar neutrinos reaching the Earth. With the projector onto the electron-neutrino flavour state, $\pi^{(e)}=\mathrm{diag}(1,0,0)$, the density matrix reads,

[TABLE]

Since the density matrix is Hermitian, $\rho_{\alpha\beta}=\rho_{\beta\alpha}^{*}$, the solar neutrino density matrix $\rho^{(e)}$ is completely characterised by the three independent $S$-matrix components,

[TABLE]

where $c_{m}$ and $s_{m}$ refer to $\cos\theta_{12}^{m}$ and $\sin\theta_{12}^{m}$, respectively. Given that we do not know precisely where neutrinos are produced in the solar core, we must average over the neutrino production positions. This effectively removes terms dependent on $\phi$ and $\Phi_{33}$ from the density matrix. The six independent density matrix elements then read,

[TABLE]

where we have defined

[TABLE]

Since neutrinos are produced within a finite volume of the Sun and the matter mixing angle, $\mathop{{}\theta^{m}_{12}}(x)$, depends on position, there is some ambiguity in what to take for the value of $\cos(2\theta^{m}_{12})$. We have taken its spatial average over the radius of the Sun as a representative value, given by

[TABLE]

where $x$ is the fractional solar radius, $p$ denotes a particular solar neutrino population, and $f_{p}(x)$ is the spatial distribution function describing where in the Sun that population is produced. These populations are labelled according to the reaction that generated them, with $p\in\{pp,\,\mathrm{{}^{8}B},\,\ldots\}$. The distributions $f_{p}(x)$ are SSM-dependent; we have used the BP16-GS98 predictions calculated by Ref. Vinyoles:2016djt . We have taken the values for each oscillation parameter from the latest NuFIT results Esteban:2020cvm .

II.4 Generalised neutrino cross sections

Following our discussion of neutrino propagation in the presence of NSI and the relevant formalism needed to derive the neutrino density matrix, $\bm{\rho}$, we move on to find expressions for the generalised scattering cross sections, ${\mathrm{d}\bm{\zeta}}/{\mathrm{d}E_{R}}$, for both neutrino-nucleus and neutrino-electron scattering.

Considering the process of elastic scattering of a neutrino $\nu$ off a target $T$ with mass $m_{T}$ via the matrix element $\mathcal{M}$, the general expression for the cross section reads

[TABLE]

From this, we can define the generalised cross section correlating the matrix elements of the flavour specific scattering processes $\nu_{\alpha}\,T\to f\,T$ and $\nu_{\beta}\,T\to f\,T$ as

[TABLE]

where we have made use of the relations $t=-2\,m_{T}\,E_{R}$ and $s=m_{T}^{2}+2\,m_{T}\,E_{\nu}$ for relativistic neutrino scattering. Note that the diagonal elements of the generalised cross section are the conventional scattering cross sections of a neutrino of flavour $\alpha$ off the target material $T$,

[TABLE]

With the general expression of Eq. 45, we can now derive the corresponding expressions for the generalised CE$\nu$NS and E$\nu$ES cross sections.

II.4.1 CE$\nu$NS cross section

Following Ref. Papoulias:2015vxa , we can derive the expression for the generalised coherent elastic neutrino-nucleus scattering cross section using the NSI formalism introduced in Section II.1. The cross section reads

[TABLE]

where $F(E_{R})$ is the Helm form factor Helm:1956zz ; Lewin:1995rx , and ${Q_{\nu N}=N-(1-4\,\sin^{2}\theta_{W})\,Z}$ is the SM coherence factor. Furthermore, we have used the Hermiticity of the NSI nucleus coupling, defined by

[TABLE]

As observed in Ref. Amaral:2021rzw , the BSM contribution can destructively interfere with the SM one. Thus, there are regions in the NSI parameter space where, despite having a non-zero NSI $\varepsilon_{\alpha\beta}$, the cross-section is the same as for the SM. In these blind spots, the presence of new physics cannot be distinguished from the SM. Since the CE$\nu$NS rate is determined by the trace over density matrix times cross section, the cancellation conditions are non-trivial. We discuss them in greater detail in Section IV.2.

II.4.2 E$\nu$ES cross section

Similarly, following Ref. Coloma:2022umy and by use of Eq. 45, we can derive the expression for the generalised neutrino-electron scattering cross section in the presence of NSI,

[TABLE]

where we have defined the generalised neutrino-electron couplings as

[TABLE]

The SM electroweak neutrino-electron couplings are given by

[TABLE]

with $g^{f}_{P}=T^{3}_{f}-\sin^{2}\theta_{w}\,Q^{\mathrm{EM}}_{f}$. In order to express the generalised neutrino-electrons coupling in terms of their vector and axial-vector components with the parameterisation of Section II.1, we introduce

[TABLE]

Effectively, this means that in Eq. 49 we can make the replacements

[TABLE]

where $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ denotes the vector component of the general NSI as before and $\tilde{\varepsilon}^{\eta,\varphi}_{\alpha\beta}$ denotes the axial-vector component (which does not contribute to matter effects and CE$\nu$NS). Note that if the NSI is only due to a vector interaction, we have $\varepsilon^{L}=\varepsilon^{R}$, such that the axial-vector component vanishes, $\tilde{\varepsilon}^{\eta,\varphi}_{\alpha\beta}=0$. As stated before, we only focus on the vector interaction for electron scattering. We do this because the results from oscillation and coherent experiments will have no impact on $\tilde{\varepsilon}^{\eta,\varphi}_{\alpha\beta}$. Furthermore, to accurately predict the signal from the axial-vector interaction, one would have to use a different ionisation form factor to that of the vector interaction.

III Extending current constraints to the full NSI parameter space

With our extended formalism in place, we are ready to explore how previous NSI results map onto the extended parameter space. Earlier constraints on NSI parameters derived from spallation source Coloma:2017ncl ; Miranda:2020tif ; Coloma:2019mbs ; COHERENT:2020iec ; COHERENT:2017ipa ; COHERENT:2021xmm and neutrino oscillation Coloma:2022umy ; Coloma:2019mbs ; Khan:2019jvr ; Esteban:2018ppq experiments have assumed that the NSI contribution in the charged plane is entirely in either the proton ($\varphi=0$) or the electron ($\varphi=\pm\pi/2$) directions. In the $\varphi=0$ case, the CE$\nu$NS cross section is maximally modified with no change to the E$\nu$ES cross section, leading to the strongest constraints on $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ from spallation source experiments and limits from oscillation experiments that only arise from non-standard propagation effects. In the $\varphi=\pm\pi/2$ case, constraints from oscillation source experiments arise from propagation effects and a maximal change to the E$\nu$ES cross section. However, the evolution of these bounds with variable charged NSI contribution has not yet been studied, and our parametrisation provides a convenient way to visualise this. Since we have no reason to believe that charged NSI would lie preferentially in any direction, a general treatment must be sought.

In this section, we recompute the bounds from spallation source and oscillation experiments, allowing for $\varphi$ to vary along its entire allowed range. In particular, we consider the CENNS-10 LAr COHERENT:2020iec and Borexino Borexino:2008gab experiments as our spallation source and oscillation experiment candidates, respectively. To demonstrate the non-trivial evolution of previously computed constraints with variable $\varphi$, we take inspiration from earlier analyses, showing how the same approaches can lead to very different results.

III.1 The CENNS-10 LAr Experiment

The CENNS-10 LAr experiment COHERENT:2020iec has measured the coherent elastic scattering of neutrinos with nuclei using a liquid argon scintillator target. The neutrino flux has three components: a prompt flux of muon neutrinos generated by the decay of pions, and two delayed fluxes of anti-muon and electron neutrinos, produced in the three-body decay of anti-muons. Their normalised spectra are given by

[TABLE]

where, from kinematics, $E_{\nu}\in[0,m_{\mu}/2]$. The expected neutrino flux is then given by scaling these spectra to account for the total beam luminosity and distance of the liquid argon target from the source. This scaling is given by $\eta\equiv r\,N_{\mathrm{POT}}/(4\pi L^{2})$, where $r$ is the number of neutrinos produced per proton collision, $N_{\mathrm{POT}}$ is the number of protons on target, and $L$ is the length of the experimental baseline. This gives us the total expected neutrino flux, $\phi_{\alpha}(E_{\nu})\equiv\eta\,f_{\alpha}(E_{\nu})$, where $\alpha\in\{\nu_{\mu},\bar{\nu}_{\mu},\nu_{e}\}$. For the CENNS-10 LAr experiment, $r=0.08$, $N_{\mathrm{POT}}=$1.37\text{\times}{10}^{23}$\,\mathrm{yr^{-1}}$, and $L=$27.5\text{,}\mathrm{m}$$ COHERENT:2020iec .

From these fluxes, we can retrieve the expected CE$\nu$NS rate spectrum. Since the neutrino beam does not undergo significant decoherence over the experimental baseline, it can be treated as being composed of independent $\nu_{\mu}$, $\overline{\nu}_{\mu}$, and $\nu_{e}$ parts. This means that the rate is given by the integral of the neutrino flux and the appropriately flavoured cross section, as it is usually written. In our notation, this reads

[TABLE]

where $M_{\mathrm{det}}=$24\text{,}\mathrm{kg}$$ is the mass of the detector, $m_{N}$ is the mass of an $\mathrm{{}^{40}Ar}$ nucleus (for which we assume $100\%$ isotopic abundance), and $\epsilon(E_{R})$ is the energy-dependent efficiency function, which we have taken from Analysis A of Ref. COHERENT:2020iec . Since this function is given in units of electron-equivalent energy ($\mathrm{\mathrm{keV}_{\mathrm{ee}}}$), we convert our spectrum into $E_{ee}$ energies before folding in the efficiency function using the energy-dependent quenching factor COHERENT:2020iec

[TABLE]

Finally, the integral over neutrino energy runs from the minimum neutrino energy required to cause a recoil of energy $E_{R}$, $E_{\nu}^{\mathrm{min}}\approx\sqrt{m_{N}E_{R}/2}$.

To compute the allowed regions for $\varepsilon^{\eta,\varphi}_{\alpha\beta}$, we perform a similar analysis to that of Ref. Miranda:2020tif , with the key difference that we allow for the charged NSI contribution to lie anywhere within the charged plane. Using a $\chi^{2}$ statistic, we compare the number of events measured by CENNS-10 LAr to the theoretical expectation given a particular choice for $\varphi$ and $\varepsilon^{\eta,\varphi}_{\alpha\beta}$. We fix $\eta=\tan^{-1}({1/2})$ in order to match the analysis of Ref. Miranda:2020tif , equivalent to having neutrino NSI with the up-quark only when $\varphi=0$. Our $\chi^{2}$ statistic is given by

[TABLE]

where $N_{\mathrm{exp}}=159$ is the number of measured events and $N_{\mathrm{bkg}}=563$ is the number of background events (primarily from the beam-related neutron rate) COHERENT:2020iec . The nuisance parameter $a$ acts as a pull parameter on the theoretical rate, allowing it to vary around its central value. This accounts for the systematic uncertainties in its calculation, and we take it to be $\sigma_{a}=8.5\%$ COHERENT:2020iec . The quadratic penalty term in Eq. 59 penalises deviations of size much greater than this.

To compute the $90\%$ CL allowed regions, we vary one NSI parameter at a time for a given angle $\varphi$ and find those values of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ for which $\Delta\chi^{2}(\varepsilon^{\eta,\varphi}_{\alpha\beta})\equiv\chi^{2}(\varepsilon^{\eta,\varphi}_{\alpha\beta},\,\varphi)-\chi^{2}_{\mathrm{min}}(\varphi)\leq 2.71$, where $\chi^{2}_{\mathrm{min}}(\varphi)$ is the minimum $\chi^{2}$ optimised over $\varepsilon^{\eta,\varphi}_{\alpha\beta}$. We repeat this analysis over the full range of $\varphi$, drawing the $90\%$ CL allowed regions in Fig. 2. We also show our $\Delta\chi^{2}$ plot for the extremal cases of $\varphi=0$ and $\varphi=\pi/2$ in Fig. 6 of Appendix B. We see that, for $\varepsilon^{\eta,\varphi}_{ee}$ and $\varepsilon^{\eta,\varphi}_{\mu\mu}$, these regions allow for two solutions: one that is consistent with the SM (i.e. $\varepsilon_{\alpha\beta}=0$) and one that is not. This first region is slightly displaced from $\varepsilon^{\eta,\varphi}_{\alpha\beta}=0$ as CENNS-10 LAr observed a slight excess of events over the SM expectation. The second region is due to a cancellation between the interference and NSI-only terms in the cross section, which can be seen by inspecting Section II.4.1 and is discussed in greater detail in the context of DD experiments in Section IV.2. While this second minimum occurs for all $\varepsilon^{\eta,\varphi}_{\alpha\beta}$, the effect is most pronounced for $\varepsilon^{\eta,\varphi}_{ee}$ and $\varepsilon^{\eta,\varphi}_{\mu\mu}$, as can be seen from Fig. 6. Importantly, no bounds can be placed on $\varepsilon^{\eta,\varphi}_{\tau\tau}$ since the CENNS-10 LAr neutrino beam has a negligible $\nu_{\tau}$ component.

Typically, the intervals that would be quoted correspond to the allowed values at $\varphi=0$. This reflects the assumption that the charged NSI lies purely in the proton direction. However, we see in Fig. 2 that these bounds generally worsen for increasing values of $|\varphi|$. While this trend is partially led by our parametrisation (whereby the strength of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ required for a constant contribution should scale as $1/\cos\varphi$ in any one of the proton, neutron, or electron directions), the bounds do not vary via this same scaling. This is particularly evident from the limits drawn for the second minima in the cases of $\varepsilon^{\eta,\varphi}_{ee}$ and $\varepsilon^{\eta,\varphi}_{\mu\mu}$, both of which deteriorate more rapidly than the first minima bounds. Moreover, the constraints on the NSI contribution from the neutron, which is inherently independent of $\varphi$ in our formalism, would worsen for increasing $|\varphi|$ (at fixed $\eta$), reflecting the requirement for a stronger NSI with the neutron to account for the diminishing contribution from the proton.

III.2 The Borexino Experiment

The Borexino experiment, located at the Laboratori Nazionali del Gran Sasso, observes solar neutrinos through their elastic scattering with electrons in its multi-ton scintillator target Borexino:2008gab . The differential scattering rate per target electron is given most generally by Eq. 1. In the case of Borexino, we consider the flux of solar neutrinos, which has contributions from different populations of electron neutrinos depending on where in the $pp$ chain or CNO cycle they are produced. We take the spectrum for each population, $\mathrm{d}\phi_{\nu_{e}}^{p}/{\mathrm{d}E_{\nu}}$, from the predictions of the B16-GS98 SSM Vinyoles:2016djt , where $p\in\{pp,\,\mathrm{{}^{8}B},\,\ldots\}$. In the case of electron recoils, the minimum neutrino energy necessary to cause a recoil of energy $E_{R}$ is given by

[TABLE]

Ignoring experimental effects, such as energy resolution and efficiency functions, the scattering rate due to a particular neutrino population, $p$, is given by

[TABLE]

where $E^{p,\mathrm{max}}_{R}$ is the maximum possible recoil energy for the population $p$ and ${\mathrm{d}R^{p}}/{\mathrm{d}E_{R}}$ is the differential rate calculated from Eq. 1. We take the number of target electrons in the scintillator to be $$3.307\text{\times}{10}^{31}$/(100\,$\mathrm{ton}$)$ Borexino:2017rsf .

We wish to explore the evolution of previous oscillation bounds in the full plane of charged NSI, i.e. with variable angle $\varphi$ at a fixed angle $\eta$. To this end, we perform a similar analysis to that of Ref. Khan:2019jvr . Namely, we consider how Borexino’s Phase-II measurements of the $pp$, $\mathrm{{}^{7}Be}$, and $pep$ solar neutrino rates Borexino:2017rsf can be used to constrain neutrino NSI with our more general formalism. In conducting our analysis, we determine the bounds on the off-diagonal matrix elements $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ ($\alpha\neq\beta$), which were not computed in Ref. Khan:2019jvr . This is only possible through the correct treatment of the differential rate in Eq. 1 using the density matrix formalism. We note that a direct comparison between our results and those of Ref. Khan:2019jvr is particularly difficult due to our different treatments of the NSI Lagrangian.

The results from the Borexino’s Phase-II run, along with the results of our calculations for the theoretical rate for each respective neutrino population, are shown in Table 1. As was done in Ref. Khan:2019jvr , we assume that the fractional uncertainties in the theoretical rates for each solar neutrino population are the same in our calculation as those reported by Borexino. Our results are in good agreement with the measured rate and the rate predicted by the collaboration Borexino:2017rsf .

To perform our statistical analysis, we construct a $\chi^{2}$ function similar to that of Section III.1:

[TABLE]

where the sum is taken over each considered solar neutrino population, $p\in\{pp,\mathrm{{}^{7}Be},\,pep\}$. The rates $R^{p}_{\mathrm{Borexino}}$ are the measured rates from the Phase-II run, with statistical uncertainties $\sigma_{\mathrm{stat}}^{p}$, while $R^{p}_{\mathrm{Theo}}$ are our calculated rates given a choice of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ and $\varphi$. We have also introduced the pull parameters $a^{p}$ for each rate, the values of which we show in Table 1. To compute our $\chi^{2}$, we profile over the nuisance parameters $\bm{a}\equiv(a^{pp},\,a^{\mathrm{{}^{7}Be}},\,a^{pep})^{\mathrm{T}}$, whose standard deviations are given in the last column of Table 1. The $90\%$ CL regions are computed via the same prescription as in Section III.1.

We show the $90\%$ CL allowed regions in Fig. 3 and the corresponding $\Delta\chi^{2}$ values for the extremal cases of $\varphi=0$ and $\varphi=\pi/2$ in Fig. 7 in Appendix B. The shapes of these regions can be understood by expanding the trace of Eq. 1 using the E$\nu$ES cross section of Eq. 49 when only one $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ is turned on at a time. The resulting formula for the rate then contains three types of terms: a propagation-only term, which contains NSI effects only at the level of neutrino propagation; a term linear in $\xi^{e}\varepsilon^{\eta,\varphi}_{\alpha\beta}$, which can be understood as an inteference term between the SM and NSI; and a positive-definite term quadratic in $\xi^{e}\varepsilon^{\eta,\varphi}_{\alpha\beta}$, which encodes the pure NSI effect in the cross section. These terms can be explicitly seen in Eqs. 71 and 72 in the context of our DD analysis.

At $\varphi=0$, the NSI effect is purely due to a change in the matter potential experienced by neutrinos on their way out of the Sun, altering neutrino propagation as per the description of Section II.2. This leads to constraints that are only due to propagation effects, with the neutrino-electron cross section unchanged. Around $\varphi=0$ and for $\varepsilon^{\eta,\varphi}_{\alpha\beta}\lesssim 1$, NSI effects remain dominated by propagation-only effects, but non-standard cross section terms linear in $\varepsilon^{\eta,\varphi}_{\alpha\beta}\xi^{e}$ begin to contribute. While the impact on the expected rate due to propagation-only effects in this regime is approximately symmetric under the exchange $\varepsilon^{\eta,\varphi}_{\alpha\beta}\rightarrow-\varepsilon^{\eta,\varphi}_{\alpha\beta}$, the effect due to the term linear in $\varepsilon^{\eta,\varphi}_{\alpha\beta}\xi^{e}$ is approximately symmetric under the combined exchange $\{\varphi,\,\varepsilon^{\eta,\varphi}_{\alpha\beta}\}\rightarrow\{-\varphi,\,-\varepsilon^{\eta,\varphi}_{\alpha\beta}\}$. This means that, depending on the sign of $\varphi$, cross section and propagation-only effects will either positively or negatively interfere with one another. For larger values of $|\varphi|$, NSI effects are predominantly due to changes in the scattering cross section, which are dominated by the term quadratic in $\xi^{e}\varepsilon^{\eta,\varphi}_{\alpha\beta}$ in the E$\nu$ES cross section for large values of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$.

These effects are perhaps best evidenced by the lower-right panel of Fig. 3. For $\varphi=0$, propagation effects alone lead to rates that are irreconcilable with the data, allowing us to constrain $\varepsilon^{\eta,\varphi}_{\mu\tau}$ without alterations to the E$\nu$ES cross section. For small negative values of $\varphi$, both cross section and propagation effects suppress the expected neutrino rate, leading to a large overall predicted deficit and a more constrained allowed region. Beyond this, terms quadratic in $\xi^{e}\varepsilon^{\eta,\varphi}_{\mu\tau}$ begin to dominate, reducing this deficit to momentarily retrieve the SM expectation, but ultimately leading to a large predicted excess. This results in the valley at negative $\varphi$ values. On the other hand, for small positive values of $\varphi$, these two effects destructively interfere with one another, resulting in a larger allowed region. For negative values of $\varepsilon^{\eta,\varphi}_{\mu\tau}$, the oscillation-only and quadratic terms reinforce one another, such that this enlarged region quickly shrinks for large $\varepsilon^{\eta,\varphi}_{\mu\tau}$ and $\varphi$. For positive $\varepsilon^{\eta,\varphi}_{\mu\tau}$, these two terms instead cancel one another out.

Additionally, Fig. 3 contains a grey region for $\varepsilon^{\eta,\varphi}_{e\mu}$ where the adiabatic approximation used to model matter effects in the Sun may be inappropriate Friedland:2000rn ; Giunti:2007ry . Within this region, the adiabaticity parameter, defined in Eq. 29, takes values $\gamma<100$, where we have calculated $\gamma$ at $E_{\nu}=$1\text{,}\mathrm{MeV}$$, approximately corresponding to the highest energy reached by $\mathrm{{}^{7}Be}$ neutrinos. We conservatively interpret these values to be in violation of the adiabaticity condition, $\gamma\gg 1$, such that a full numerical calculation of the density matrix elements would be required for an accurate analysis. This would be beyond the scope of our work, and since the allowed NSI regions in Fig. 3 are almost entirely within the adiabatic regime, we do not believe a numerical treatment is necessary. We have checked that the we fulfil the adiabatic criterion for all other $\varepsilon^{\eta,\varphi}_{\alpha\beta}$.

Our analysis also shows that constraints are strongest for off-diagonal NSI. This is because the trace in Eq. 1 leads to two terms (which are equal in our case as $\delta_{\mathrm{CP}}$ is set to [math]) contributing to the total NSI rate in the off-diagonal case, as opposed to a single contribution arising from diagonal NSI. Thus, the allowed regions for off-diagonal NSI are generally tighter than those for diagonal NSI. The exception to this is in the bounds for $\varepsilon^{\eta,\varphi}_{ee}$, which are highly constrained due to the enhanced E$\nu$ES cross section arising from the additional CC contribution via the $W$-boson exchange.

This enhanced cross section not only leads to much tighter bounds for $\varepsilon^{\eta,\varphi}_{ee}$ but also allows for a finely tuned second minimum at non-zero $\varepsilon^{\eta,\varphi}_{ee}$. With a non-zero $\varepsilon^{\eta,\varphi}_{ee}$, the differential rate spectrum is modified such that the total rate, given by integrating over all recoil energies, coincidentally retrieves the SM expectation. The spectrum itself, however, is significantly modified, and incorporating spectral information into our analysis would ultimately prohibit this second solution.

Our results should not be taken as a dedicated Borexino analysis. Though we have attempted to capture the variation in the calculated solar neutrino rates by introducing pull parameters, as was done in Ref. Khan:2019jvr , this only accounts for the theoretical uncertainty in these rates. Ultimately, a more sophisticated analysis would require a spectral fit of Borexino’s data and allow for multiple NSI parameters to vary at a time. Such a fit should then allow for the various background components inherent in the data to float and permit correlations between all fit parameters. Such an analysis was recently done in the context of neutrino NSI by Ref. Coloma:2022umy assuming charged contributions only from the electron ($\varphi=\pi/2$). Our results should instead be taken as a demonstration of how bounds on NSI can vary dramatically depending on what one takes as the underlying charged NSI contribution.

IV Direct Detection Experiments

Finally, we turn to the main motivation of this paper: determining the potential of DD experiments within the NSI landscape using the extended framework introduced in Section II.

DD experiments will provide a unique probe of neutrino interactions because they have access to both CE$\nu$NS and E$\nu$ES. This is due to the properties of the solar neutrino flux. To produce a recoil in the energy range detectable in DD experiments ($\sim 10$ eV$-100$ keV) E$\nu$ES and CE$\nu$NS probe different solar neutrino populations. In particular, the main contribution to NRs comes from 8B neutrinos, whereas for ERs it is $pp$ and $\mathrm{{}^{7}Be}$ that contribute the most. Thus, although the scattering cross-section for CE$\nu$NS is significantly larger, the much larger $pp$ neutrino flux compensates for the smaller E$\nu$ES cross section.

As we saw in Section III.2, experiments such as Borexino have optimal sensitivity when NSI occur purely with the electron, but they rapidly lose their constraining power as $\varphi\rightarrow 0$. On the other hand, CE$\nu$NS experiments such as CENNS-10 LAr, which we explored in Section III.1, have excellent sensitivity when the charged NSI contribution is wholly in the proton direction, but they lose this sensitivity as this contribution turns to the electron. Moreover, having no $\nu_{\tau}$ component, they are completely insensitive to $\varepsilon^{\eta,\varphi}_{\tau\tau}$. Not only can DD experiments probe $\varepsilon^{\eta,\varphi}_{\tau\tau}$, but their ability to measure and discriminate between NRs and ERs means that they retain their constraining power across the full range of $\varphi$. When applied to specific BSM models, this can be crucial to identify the underlying nature of the new physics scenarios (see e.g., Ref. Amaral:2021rzw ).

In this work, we focus on the xenon-based DD experiments LZ Mount:2017qzi ; LUX-ZEPLIN:2018poe ; LZ:2021xov , XENON XENON:2018voc ; XENON:2020kmp ; XENON:2020rca , and DARWIN DARWIN:2016hyl . More specifically, we derive exclusions from the data reported by the recent LZ WIMP search LZ:2022ufs and the XENONnT electron-recoil excess search XENONCollaboration:2022kmb . We also determine the expected sensitivities of LZ, XENONnT and DARWIN, based on their projections for their final experimental configurations. Similar results can be obtained with PandaX PandaX-II:2017hlx ; PandaX:2018wtu ; PandaX-4T:2021bab .

IV.1 Expected Number of Events and Statistical Procedure

To calculate our sensitivities, we consider the differential rate of Eq. 1 and incorporate detector effects, such as efficiency and energy resolution. The expected recoil rate from neutrino scattering is then given by

[TABLE]

Here, ${\mathrm{d}R}/{\mathrm{d}E_{R}^{\prime}}$ is given by Eq. 1. When computing this rate, we use the solar neutrino fluxes predicted by the B16-GS98 model Vinyoles:2016djt , as we did for our Borexino analysis in Section III.2. The integral over the expected energy, $E_{R}^{\prime}$, is the convolution that describes the effect that the detector resolution, $\sigma$, has on the observed signal; we assume this to be equivalent to a Gaussian smearing. This resolution is typically reported in terms of the measured, electron-equivalent energy, so we first convert the CE$\nu$NS differential rate and NR efficiency functions into electron-equivalent energies when considering NRs. We do this by applying an energy-dependent quenching factor, which relates the two energy scales via $E_{\mathrm{ee}}=Q(E_{\mathrm{nr}})\,E_{\mathrm{nr}}$. We take this to be the Lindhard factor Lindhard1963INTEGRALEG with $k=0.1735$, reflecting the $k$-value found in the fit performed by the LUX collaboration LUX:2016ezw . Finally, $\epsilon$ is the energy-dependent efficiency function.

The differential rate depends on the number of targets per unit mass of the detector, $n_{T}=N_{T}/m_{\mathrm{det}}$. For NRs, we take this to be the number density of atoms in the detector, $n_{T}=1/m_{N}$, where $m_{N}$ is the nuclear mass of the relevant xenon isotope. For ERs, this corresponds to the number of ionisable electrons given a recoil of energy $E_{R}$, scaled to agree with the ab initio calculations from the relativistic random-phase approximation for xenon Chen:2016eab . This takes into account the many-body dynamics involved in such collisions, and they have shown a consistent suppression of the rate at low recoil energies Chen:2016eab ; Hsieh:2019hug .

To compute the number of expected events within the $i^{\mathrm{th}}$ bin, we integrate the differential rate in the energy window defined by the edges of the bin, $[E^{i}_{1},\,E^{i}_{2}]$, and sum the contributions from each nuclear isotope $A$ multiplied by its corresponding relative isotopic abundances, $X_{A}$:

[TABLE]

Here, $\varepsilon$ is the experimental exposure and ${\mathrm{d}R_{A}}/{\mathrm{d}E_{R}}$ is the differential rate in Eq. 63 due to isotope $A$.

We determine our sensitivities using a series of log-likelihood-ratio tests in which we vary only one NSI parameter at a time, fixing all others to zero666A global analysis in which all parameters are allowed to vary is beyond the scope of this article, where our aim is to motivate next-generation and far-future DD experiments to be included in future such analyses.. We construct our likelihoods from a Poisson part and a Gaussian part, which we use to capture the effect of uncertainties on nuisance parameters. For this latter part, we consider Gaussian distributed pull parameters serving to scale the number of expected neutrino events, as we did for CENNS-10 LAr and Borexino, and the number of expected background events. We respectively label these parameters as $a$ and $b$, having standard deviations $\sigma_{a}$ and $\sigma_{b}$. Given some number of observed events in bin $i$, $N^{i}_{\mathrm{obs}}$, we define the likelihood function

[TABLE]

where $N^{i}_{\mathrm{bkg}}$ is the number of expected background events in the $i^{\mathrm{th}}$ bin. The product is over $N_{\mathrm{bins}}$ bins. The number of observed events for each of our analyses depends on whether we compute exclusions based on data or derive projected limits. If the former, we take the number of observed neutrino events reported in each bin. If the latter, we assume an Asimov data set Cowan:2010js , such that $N^{i}_{\mathrm{obs}}$ is set to the number of expected SM neutrino events ($\varepsilon^{\eta,\varphi}_{\alpha\beta}=0$) in the $i^{\mathrm{th}}$ bin.

Finally, to derive our limits, we use Eq. 65 to define the test statistic

[TABLE]

where hatted variables indicate quantities that maximise the likelihood given the parameter $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ and double-hatted variables indicate those quantities that maximise the overall, unconstrained likelihood. By fixing the angles to take some values $\eta=\eta_{0}$ and $\varphi=\varphi_{0}$, we constrain only the parameter of interest, $\varepsilon^{\eta,\varphi}_{\alpha\beta}$, in each of our analyses. In practice, we perform our analysis by only considering the most dominant nuisance parameter (either $a$ or $b$), which depends on the experiment and is discussed in detail below. The $90\%$ CL limits are then calculated by finding the value for $q_{\varepsilon}$, $q^{\mathrm{lim}}_{\varepsilon}$, for which

[TABLE]

where $\mathop{{}f}(q_{\varepsilon})$ is the distribution of the test statistic. In the limit of high statistics and that the true parameter value does not lie on the boundary of our parameter space, Wilks’ theorem tells us that this distribution asymptotically follows a $\chi^{2}$-distribution with number of degrees of freedom $k=1$777We have checked that this holds true in all of our analyses.. This leads to $q_{\varepsilon}^{\mathrm{lim}}=2.71$, and our limits then follow from finding that $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ which yields this value for the test statistic.

To implement the analysis in the sections that follow, we make use of SNuDD SNuDD-code 888https://github.com/SNuDD/SNuDD.git (Solar NeUtrinos for Direct Detection), a novel code-base that we have developed. SNuDD is a Python package that calculates the generalised cross section of Section II.4 and the density matrix elements of Section II.3, combining them to compute the trace of Eq. 1 and arrive at a prediction for the expected solar neutrino rate at a DD experiment while folding in detector efficiency and resolution effects. We will release it in a separate publication. We hope that SNuDD facilitates future DD analyses in the NSI landscape.

Lastly, we have checked that our limits are not greatly impacted by uncertainties in the neutrino oscillation parameters when these uncertainties are computed using standard neutrino oscillations and interactions. While including NSI can significantly impact the best-fit values of and uncertainties in these parameters, especially those derived from solar neutrino experiments Esteban:2018ppq , we note that the future medium-baseline reactor experiment JUNO JUNO:2015sjr will be largely insensitive to the effects of NSI. In particular, the fits to the solar neutrino parameters, $\theta_{12}$ and $\Delta m_{12}^{2}$, and the reactor neutrino angle, $\theta_{13}$, will be robust in the presence of NSI Martinez-Mirave:2021cvh . Moreover, the atmospheric angle ($\theta_{23}$), which impacts the muon and tau neutrino fractions, is not greatly changed by the inclusion of NSI Esteban:2018ppq . Thus, we take the liberty of ignoring these systematics in this first study, taking the best-fit values for these parameters from Ref. Esteban:2020cvm ; a more comprehensive analysis that includes these systematics is left for future work.

IV.2 Sensitivities in the nucleon NSI plane

From the many existing DD NR constraints, we consider only the recent, leading LZ result LZ:2022ufs to derive a bound in the NSI landscape. We take the efficiency function given in Ref. LZ:2022ufs and we model the energy resolution according to Ref. LUX:2016rfb . For our signal region, we use the 90$\%$ quantile of the nuclear recoil band in the S1$c$ (scintillation) and S2$c$ (ionisation) event reconstruction space as shown in Fig. 4 of Ref. LZ:2022ufs . For our analysis, we focus only on the low-energy region ($[5,\,15]\,{\rm keV}$), which is sensitive to solar neutrinos, and we integrate over it to constitute one bin. We take the number of expected background and observed events to be 1 and 0, respectively (see Ref. Chang:2022jgo ). We have validated our procedure by reproducing the WIMP-nucleon cross-section limit reported by the collaboration (Fig. 5 of Ref. LZ:2022ufs ) to a good agreement in the low and high DM mass values.

To assess the future prospects of detecting an NR signal, we plot the projected sensitivities for LZ, XENONnT, and DARWIN with the full exposures of $15.34$, $20$, and $200$ $\mathrm{ton}\text{\,}\mathrm{yr}$, respectively. We have taken a background-free scenario in the NR search, which is consistent with experimental aims (for our analysis, the neutrino signal is not considered a background). With no backgrounds, we then conservatively perform a one-bin analysis, using the total number of expected solar neutrino NR events as our observation.Additionally, since the expected background is so low, the nuisance parameter that will have the largest impact on our sensitivities will be associated with the solar neutrino flux, and hence we set $b=0$ in this case. For our remaining pull parameter in Eq. 65, we assume an uncertainty of $\sigma_{a}=12\%$, reflecting the $12\%$ uncertainty in the theoretical value of the total $\mathrm{{}^{8}B}$ flux in the B16-GS98 SSM Vinyoles:2016djt . We take the resolution function for LZ at full exposure to be the same as that of their first result LZ:2022ufs ; LUX:2016rfb , whereas for XENONnT and DARWIN we use the resolution function given in Ref. XENON:2020rca . The NR efficiency functions, as presented by the collaborations, reach $50\%$ at $3.8\,\mathrm{keV}_{\mathrm{nr}}$ for LZ LZ:2018qzl and $5.7\,\mathrm{keV}_{\mathrm{nr}}$ for XENONnT and DARWIN XENON:2018voc . However, to explore how DD experiments could feasibly probe NSI in the future, we take the liberty of further lowering these thresholds. This is to take advantage of the higher $\mathrm{{}^{8}B}$ rate at these energies.

In particular, we augment the efficiency functions such that, for each future experiment, the efficiency instead reaches $50\%$ at $3\text{\,}\mathrm{\mathrm{keV}_{\mathrm{nr}}}$, which we consider to be a feasible future goal. For instance, the xenon-based LUX experiment was able to reach thresholds as low as $1.1\text{\,}\mathrm{\mathrm{keV}_{\mathrm{nr}}}$ while retaining NR/ER discrimination LUX:2015abn . The LUX collaboration has also developed techniques allowing for single-photon sensitivities, resulting in sensitivities to much lower recoil energies at the cost of a lower overall detection efficiency. Finally, the XENON1T collaboration has recently performed a dedicated $\mathrm{{}^{8}B}$ search by lowering their threshold to $1.6\text{\,}\mathrm{\mathrm{keV}_{\mathrm{nr}}}$, achieved by relaxing the necessity for a three-fold S1 coincidence in the PMTs to a two-fold one XENON:2020gfr . Furthermore, taking the systematic 8B uncertainty to be 12%, we find that lowering the threshold further provides little-to-no benefit in terms of NSI sensitivity. For each future experiment, we take $E^{\max}_{R}=30\,{\rm keV}$.

We show our results in Fig. 4. The shaded areas represent the 90% CL limits set by the different experimental configurations. From less constraining (smaller areas) to more constraining (larger areas), we show the limits derived from the first LZ results (turquoise with solid boundary) and the expected sensitivities of the full exposure of LZ (baby blue, dashed), XENONnT (dark blue, dashed) and the proposed DARWIN (purple, dashed). For comparison, we also show with red bars the NSI limits derived from the global study of Ref. Coloma:2019mbs , which included the results from COHERENT and a variety of oscillation experiments, when NSI take place purely with the proton ($\eta=0$), up-quark ($\eta=\tan^{-1}(1/2)$) and down-quark ($\eta=\tan^{-1}(2)$). We extract the limits we have computed using the LZ WIMP search data and tabulate them in Table 2, contrasting them with the results from the global fits of Coloma:2019mbs . We see that, currently, DD experiments are not sensitive to globally allowed NSI values, but our projections indicate that they will be in the near future.

Like in Fig. 3 of our Borexino analysis in Section III.2, Fig. 4 contains grey regions, indicating those points in the parameter space where the adiabatic approximation may be invalid. Within these regions, $\gamma<100$, where in this case we have calculated $\gamma$ at $E_{\nu}=$16\text{,}\mathrm{MeV}$$, approximately corresponding to the highest energy reached for $\mathrm{{}^{8}B}$ neutrinos. As current global fits show that the allowed values of the NSI parameters are firmly within the adiabatic regime, we believe that our analytical approach is valid for the regions of interest. However, it is important to keep in mind that our sensitivities may be inaccurate within the grey bands.

Our limits exhibit many interesting non-trivial features. Specifically, we see that there are regions in each NSI parameter space where every DD experiment loses sensitivity. The two most remarkable of these regions are, firstly, the strong cancellation in the angle $\eta$ occurring at $\eta\approx-35^{\circ}$ and, secondly, the band of insensitivity in $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ across the entire range of $\eta$ values (made manifest by the gaps in the projected sensitivity areas). These blind spots present a challenge for DD experiments, and they should be understood if DD experiments are to maximise their constraining power in the NSI landscape.

We first consider the cancellation in $\eta$, which occurs at the same point regardless of the nature of the NSI. From Section II.4.1 and Section II.4.1, we see that the non-standard contribution to the CE$\nu$NS cross section vanishes when $\xi^{p}Z+\xi^{n}N=0$ , recovering the SM cross section regardless of the value of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$. For a given nuclear isotope, this occurs when999This cancellation was used in direct dark matter detection to argue that dark matter particles might escape detection in some specific targets (see e.g., Ref. Feng:2013vod ).

[TABLE]

This condition depends on the choice of target material. For composite materials or targets, or non-mononuclidic targets, the cancellation is not exact, as the contributions from the different isotopes must be added up in Eq. 64. Yet, for xenon the ratio $Z/N$ is very similar in all its natural isotopes, and the observed rate is greatly reduced for $\eta\approx-35^{\circ}$ when $\varphi=0$. Since stable nuclei tend to have similar $Z/N$ fractions, the position of the blind spot does not vary greatly for different target choices. Interestingly, argon (a target employed in current detectors DarkSide:2014llq and planned tonne-scale ones DarkSide-20k:2017zyg ) leads to a considerable shift, with $\eta\approx-39^{\circ}$. This could lead to a noticeable effect if a full spectral analysis is performed of the observed signal, thus strengthening the notion of complementarity among DD targets. This would be even more noticeable in light nuclei, such as He, F, Na, or Si, for which $Z/N\sim 1$ and the cancellation takes place for $\eta\rightarrow 45^{\circ}$ (although a large detector would still be needed).

The second blind spot occurs at intermediate values of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$, stretching over the full range of $\eta$ values. These insensitivity bands arise due to interference effects, where the NSI contribution is cancelled and thus the SM CE$\nu$NS differential rate is restored. The exact location of these bands differs for flavour-conserving and flavour-violating NSI.

In the case of flavour-conserving NSI, we can derive a simple analytical formula for the values of the NSI parameters leading to a non-trivial realisation of the SM differential rate. This relation, which defines the centres of each of these insensitivity bands where the NSI contribution exactly cancels, has previously been pointed out in Ref. Dutta:2020che and in our framework is given by

[TABLE]

The dependence on $\eta$, encoded in $\xi^{p}$ and $\xi^{n}$, gives us the band over different values of $\varepsilon^{\eta,\varphi}_{\alpha\alpha}$ as a function of $\eta$. Note that, as with the first blind spot, the locations of these bands depend on the choice of the target material due to the dependence on $Z$ and $N$. For $\eta=0$, for instance, Eq. 69 gives $\varepsilon^{\eta,\varphi}_{\alpha\alpha}\approx 0.6$ for xenon, whereas it yields the lower $\varepsilon^{\eta,\varphi}_{\alpha\alpha}\approx 0.5$ for argon. Since the non-trivial cancellation occurs at different values for different targets, this could be important in determining whether global minima are driven by data or just artefacts of the blind spot that nuclei have, see for instance Ref. Chaves:2021pey . Considering different materials thus gives us one possible avenue to mitigate this particular loss of sensitivity, though we note that the blind spots of xenon and argon move closer together as $\eta\rightarrow\pi/2$.

In the case of flavour-changing NSI, the cancellation condition becomes more complicated and is only retrieved in the (correct) basis-independent formulation of the scattering rate in terms of the trace $\mathrm{Tr}\left[\bm{\rho}\ {\mathrm{d}\bm{\zeta}}/{\mathrm{d}E_{R}}\right]$ in Eq. 1. Due to the flavour-coherence effects, we still expect regions where the SM-NSI interference term cancels the NSI-only term; however, these regions now also depend on the density matrix elements. To investigate this behaviour, we consider, as a simplification, the values of $\varepsilon_{\alpha\beta}^{\eta,\varphi}$ for which the differential rate spectrum returns to its expected SM value for a given recoil energy $E_{R}$. This prescription removes the need to integrate over $E_{R}$ to find the number of events. From Eqs. 1 and II.4.1, we find that in general, for $\alpha\neq\beta$, the condition for restoring the SM rate reads

[TABLE]

The difference in the forms of the relations in Eqs. 69 and 70 is why the positions of these bands are different for flavour-conserving and flavour-violating NSI. In particular, we note that, unlike in the case of the former, for the latter we have a flavour dependence through the appearance of the density matrix elements. This is why, for instance, we see a sign flip of the bands in the case of $\varepsilon^{\eta,\varphi}_{e\tau}$ and $\varepsilon^{\eta,\varphi}_{\mu\tau}$ with respect to $\varepsilon^{\eta,\varphi}_{e\mu}$, as the relevant off-diagonal density matrix elements $\rho_{e\tau}$ and $\rho_{\mu\tau}$ are negative in contrast to $\rho_{e\mu}$.

As we can see from Fig. 4, the insensitivity bands for the off-diagonal NSI elements $\varepsilon^{\eta,\varphi}_{e\mu}$ and $\varepsilon^{\eta,\varphi}_{e\tau}$ exhibit some interesting features at $\eta\approx-5/16\,\pi$, where they seem to expose a kink. The origin of these kinks in the off-diagonal NSI insensitivity bands can be traced back to the appearance of the off-diagonal density matrix elements, $\rho_{\alpha\beta}$, in the CE$\nu$NS cancellation condition. These kinks arise because the last term in Eq. 70 proportional to $\rho_{\alpha\beta}$ undergoes a qualitative change of behaviour at $\eta\approx-5/16\,\pi$.

We will discuss the behaviour of the kinks using the insensitivity band for $\varepsilon^{\eta,\varphi}_{e\mu}$ in the top right plot of Fig. 4 as an example. At very negative angles $\eta\approx-\pi/2$, the off-diagonal interference term in Eq. 70 proportional to $\rho_{\alpha\beta}$ has an extremum for positive $\varepsilon^{\eta,\varphi}_{\alpha\beta}$. In this regime, the behaviour of the cancellation line (which occurs at negative $\varepsilon^{\eta,\varphi}_{\alpha\beta}$) is entirely dominated by the NSI-only term proportional to the diagonal density matrix elements, $\rho_{\alpha\alpha}$. However, at larger angles, $\eta\approx-5/16\,\pi$, the extremum in the off-diagonal term shifts from positive values of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$ to negative values, and it hence begins to dominate the behaviour of the cancellation bands in Eq. 70. This change of behaviour in the cancellation integral leads to the appearance of the kinks in the insensitivity band.

For $\varepsilon^{\eta,\varphi}_{e\tau}$, the same effect leads to the appearance of a kink, but with opposite signs of $\varepsilon^{\eta,\varphi}_{\alpha\beta}$. In principle, the same reasoning holds for $\varepsilon^{\eta,\varphi}_{\mu\tau}$; however, in this case, the off-diagonal term exhibits an almost negligible extremum, such that there is no visible kink. Finally, the fact that the behaviour of the off-diagonal density matrix element $\rho_{\alpha\beta}$ is responsible for the appearance of these kinks also explains why they are absent for the diagonal NSI elements, $\varepsilon^{\eta,\varphi}_{\alpha\alpha}$, since there is no contribution from $\rho_{\alpha\beta}$.

IV.3 Sensitivities in the charged NSI plane

While the E$\nu$ES cross section will only be modified when $\varepsilon^{\eta,\varphi}_{\alpha\beta}\neq 0$ and $\varphi\neq 0$, even in the case of pure nuclear NSI couplings ($\varphi=0$), propagation effects within the solar medium can still alter the expected ER rate in DD experiments. Thus, since the charged plane contains proton NSI modifications, both ER and NR signals must be considered. There is only one direction in which a non-zero NSI will not affect the NR signal, and that is precisely the electron-only direction, $\varphi=\pm\pi/2$ and $\eta=0$. For this reason, in Fig. 5 we include both the NR and ER analysis to show the projected sensitivities of DD when NSI lie in the $(\varepsilon^{p}_{\alpha\beta}\varepsilon^{e}_{\alpha\beta})$-plane.

Currently, the world-leading DD constraint on ERs comes from XENONnT XENONCollaboration:2022kmb , which has thus far reached an exposure of $1.16\text{\,}\mathrm{ton}\text{\,}\mathrm{yr}$. We have replicated this analysis by taking the efficiency function, expected backgrounds, and observed number of events from Ref. XENONCollaboration:2022kmb , where, for the background model $B_{0}$, we subtract their expected solar neutrino background (Fig. 4 of Ref. XENONCollaboration:2022kmb ). In the signal region of $[0-140]\,${\rm keV}$$, the SM counts predicted by SNuDD is $274$, which is lower than the quoted $300\pm 30$ XENONCollaboration:2022kmb . This could be due to the fact that the neutrino signal in Ref. XENONCollaboration:2022kmb uses a simplified modelling of the neutrino spectrum. Specifically, SNuDD uses the relativistic random phase approximation (RRPA) studied in Ref. Chen:2016eab along with a series of step-functions to model the effect of electron binding energies in xenon. This reduces the overall rate and introduces discontinuous jumps in the spectra when more electrons can be ionised above certain energies. One such discontinuity is around $\sim 30\,{\rm keV}$, which does not appear to be present in Fig. 4 of Ref. XENONCollaboration:2022kmb . When we remove both RRPA and the step function approximation in SNuDD to determine the expected number of solar neutrino events with the same setup as Ref. XENONCollaboration:2022kmb , we predict $299$ solar neutrino events.

Unlike in the NR case, the sizeable backgrounds for ERs mean that spectral information should be used to harness greater sensitivity. Following XENONnT, our analysis uses $2\,{\rm keV}$-width bins from $[0-30]\,{\rm keV}$. We have refrained from using the entire signal region because the backgrounds are at their lowest at low energies. Additionally, the backgrounds below $30\,{\rm keV}$ are dominated by one source, 214Pb, which has an associated uncertainty that we treat as a nuisance parameter. If we consider higher recoil energies, backgrounds such as 124Xe, 83mKr and 136Xe become important, all of which have different associated uncertainties. A dedicated experimental analysis would include all backgrounds and their uncertainties to perform a multivariate fit to the observed events. We believe such an in-depth study should be done in consort with the collaboration. The uncertainty we take for the 214Pb dominated background is $\sigma_{b}=12.5\%$ XENONCollaboration:2022kmb . We consider this to be our dominant nuisance parameter and find that if we instead perform the fit assuming the $pp$ neutrino flux is the dominant nuisance parameter ($\sigma_{a}=1\%$ Coloma:2022umy ) our limits for XENONnT see a substantial improvement.

As mentioned above, the potential for future ER analyses relies primarily on the anticipated background reduction. For the full XENONnT run, we take the backgrounds from Ref. XENON:2020kmp , for LZ we use Ref. LZ:2021xov , and for DARWIN we use the predictions given in Ref. Baudis:2013qla . Unlike with the NR signal, the ER neutrino spectrum does not fall off sharply at $E_{R}\sim\rm{keV}$, so we do not extend the ER efficiency functions to lower energies as we did for our NR projections. We use the efficiency functions given in Ref. LZ:2018qzl for LZ, and in Ref. XENON:2020rca for XENONnT and DARWIN, where they reach $50\%$ at $1.46\,\rm{keV}_{ee}$ and $1.51\,\rm{keV}_{ee}$, respectively. For these projections, we also perform spectral analyses, binning with $2\,{\rm keV}$-width bins in the energy ranges $[0-60]\,{\rm keV}$ for XENONnT and DARWIN, but in the range $[0-30]\,{\rm keV}$ for LZ. LZ’s maximum is limited by their reported efficiency function LZ:2018qzl . Additionally, we have assumed that these experiments will have a greater understanding of their backgrounds and therefore consider the dominant nuisance parameter to be the $pp$ neutrino flux, $\sigma_{a}=1\%$. We believe this is an achievable goal for future DD experiments and see our projected sensitivities as an additional motivation for improved understanding and reduction of backgrounds. We also considered a far-future xenon detector with an exposure of $10^{3}~{}$\mathrm{ton}\text{,}\mathrm{yr}$$ as in Ref Dutta:2020che and found that the $pp$ flux uncertainty drives the projected sensitivity to the extent that even with five times the exposure of DARWIN and no backgrounds, only marginal improvements are made to the sensitivities.

We show our results in Fig. 5, where we have filled contours for ER (pink/red colours, dotted lines for projections) and NR (blue/purple colours, dashed lines for projections) analyses. In order to place these sensitivities in the wider experimental context, we take the recent results of Ref. Coloma:2022umy , which used the spectral data from Phase-II of the Borexino experiment Borexino:2017rsf to constrain $\varepsilon^{e}_{\alpha\beta}$. As Ref. Coloma:2022umy does not mention the potential impact of either proton or neutron NSI on neutrino oscillations, we assume that they have only considered NSI with the electron, with no contribution from either the proton or the neutron. As a result, we set $\eta=0$, and we place their bounds at $\varphi=\pm\pi/2$, corresponding to electron-only NSI.

We note that, while previous studies have also constrained electron NSI, most of them place individual bounds on the left- and right-handed components of the interaction Bolanos:2008km ; Agarwalla:2012wf ; Gonzalez-Garcia:2011vlg ; Davidson:2003ha ; TEXONO:2010tnr . For comparison with the previous literature, we tabulate many of these results alongside the corresponding allowed intervals derived in this work from the XENONnT ER data in Table 3. In this context, it is worth noting that Ref. Dutta:2020che previously considered the potential impact of including E$\nu$ES data from DUNE and a theoretical high-exposure DD experiment on global fit results of electron NSI. We investigate this impact in more detail by computing the solar neutrino scattering rate via the coherent treatment of oscillation effects in the density matrix formalism in Eq. 1, considering both CE$\nu$NS and E$\nu$ES, and studying the non-trivial behaviour of DD sensitivities in the full plane of charged NSI $(\varepsilon^{e}_{\alpha\beta}\varepsilon^{p}_{\alpha\beta})$ by means of our parametrisation in Fig. 5. Finally, comparing our XENONnT limits and future projections to the limits derived using Borexino data in Refs. Dutta:2020che and Coloma:2022umy , from Table 3 we see that while current DD data sets are not able to yield competitive bounds, next-generation and far-future DD experiments will be able to improve on current limits.

As in Figs. 3 and 4, we show the points of parameter space where the adiabatic limit may no longer be valid for neutrino propagation in the Sun. Since $\gamma$ is energy-dependent, this region is different for the values of $E_{\nu}$ probed by NR and ER analyses. In light grey, we show the regions relevant for NRs ($E_{\nu}=16\,{\rm MeV}$ as in Fig. 4) and in dark grey we show the regions for ERs ($E_{\nu}=1\,{\rm MeV}$), which roughly correspond to the highest energy of 7Be neutrinos. This is not the primary neutrino source for the E$\nu$ES signal ($pp$), but it does contribute at higher energies. Taking this value is a conservative choice since higher energies correspond to a greater violation of adiabaticity. This is reflected in the fact that the dark grey regions, if present, are contained within the light grey regions in Fig. 5.

Fig. 5 demonstrates that next-generation and far-future DD experiments will form powerful probes of electron NSI, with almost all of our projections cutting into portions of the bounds placed with the Borexino experiment. DARWIN can give us considerably more sensitivity to all NSI parameters, showcasing its excellent potential in searching for new physics in the neutrino sector. However, such potential exhorts substantial efforts in background modelling and reduction, something which is already well underway in the respective collaborations.

As in the CE$\nu$NS case, the limits for the E$\nu$ES case exhibit blind spots where the predicted rate is indistinguishable from the SM expectation. Once again, this weakens the limits at certain values of $\varphi$ and a series of bands where DD experiments appear to lose sensitivity. Here, the complementarity with the NR analysis can be seen explicitly, since at precisely $\varphi=0$, the NSI effect on the NR signal is maximal. However, there are two notable physical differences between the E$\nu$ES case and the CE$\nu$NS case, arising from both the different CE$\nu$NS and E$\nu$ES cross sections and the way in which non-standard matter effects enter.

Firstly, we have no strong cancellation in the ER limits. While one might expect a complete loss of sensitivity when $\varphi=0$, where the E$\nu$ES cross section is unchanged by the presence of NSI, neutrino oscillations are still impacted by the NSI contribution to the matter Hamiltonian from the nucleons (in this case only the proton since $\eta=0$). Thus, for high enough values of $\varepsilon_{\alpha\beta}^{\eta,\varphi}$, the effect of NSI on the neutrino flavour fractions is large enough to give us an observable deviation from the SM expectation. This is analogous to what we saw in our Borexino analysis of Section III.2. Consequently, while we do lose sensitivity in ERs as $\varphi$ approaches zero, our limits ultimately reach a finite value.101010Note that this is only possible as the cross section for electron neutrinos contains the extra CC contribution, making it different from that of the muon and tau neutrinos. Changes in the electron neutrino fraction then lead to measurable changes in the total number of CC interactions in the detector; the NC interactions from all flavours, on the other hand, remain equal.

Secondly, we have fewer bands of insensitivity in ER over $\varphi$ than we did for the NR case over $\eta$. The location of these bands can be calculated through identical arguments to the CE$\nu$NS case, whereby those values of the NSI parameters where the NSI-augmented rate is equal to the expected SM rate, $\mathrm{d}R/\mathrm{d}E_{R}-\mathrm{d}R/\mathrm{d}E_{R}|_{\mathrm{SM}}=0$, are found. The derivation of these cancellation conditions again crucially hinges on the coherent treatment of the neutrino propagation via the density matrix formalism in Eq. 1. Critically, the condition for cancellation in the off-diagonal NSI elements $\varepsilon_{\alpha\beta}$ is completely missed in the simplified treatment of the rate as the sum over the oscillation probabilities times scattering cross section, $\sum_{\alpha}P_{e\alpha}\,\mathrm{d}\sigma_{\nu_{\alpha T}}/\mathrm{d}E_{R}$.

In the case that only one diagonal NSI element, $\varepsilon_{\alpha\alpha}^{\eta,\varphi}$, is active, the cancellation equation for the differential rate for E$\nu$ES reads,

[TABLE]

On the other hand, for off-diagonal NSI, $\varepsilon_{\alpha\beta}^{\eta,\varphi}$ with $\alpha\neq\beta$, the cancellation condition can be expressed as

[TABLE]

where the last term in the second line is only present for $\alpha\beta=\mu\tau$. In the above expressions, $\rho^{\mathrm{SM}}$ refers to the density matrix obtained in the SM case (i.e. $\varepsilon_{\alpha\beta}=0$) and $\rho$ to the one obtained with non-zero NSI elements. As can be seen in Fig. 5, for E$\nu$ES we obtain insensitivity bands similar to those in the CE$\nu$NS case for the projected limits; however we only see this for $\varepsilon_{ee}^{\eta,\varphi}$ and $\varepsilon_{e\mu}^{\eta,\varphi}$. In most cases, DD sensitivities are not good enough to reach these cancellation regions where the SM scattering rate is recovered.

Finally, it is worth noting that, for very small but non-zero $\varphi$, the line of exact cancellation has a very sharp zero-transition from very large (positive) to very small (negative) values of $\varepsilon_{\alpha\beta}^{\eta,\varphi}$ (or vice versa) making it seem like there exists an asymptote at $\varphi=0$. This effect is similar to that exhibited by Fig. 3 in our analysis of Borexino data. In this region of parameter space, the NSI-only term is negligible since it is quadratic in $\xi^{e}$, and thus $\varphi$. The reason for this change in sign is a rapid flattening of the SM-NSI interference terms in Eqs. 71 and 72, which are linearly proportional to $\xi^{e}$, and thus $\varphi$. This flattening leads to a rapid change in the value of $\varepsilon_{\alpha\beta}^{\eta,\varphi}$ where the interference term cancels off the residual SM-like term proportional to $\rho_{ee}-\rho^{\mathrm{SM}}_{ee}$ and hence restores the SM neutrino rate.

The bounds from the global analysis of Ref. Coloma:2019mbs are shown in Fig. 5 as they were in Fig. 4, but now only the proton direction ($\eta=0,\varphi=0$) is visible. By plotting the NR analyses, we can see how the sensitivities of Fig. 4 extend into the $\varphi$ direction. We observe the same regions of nonzero $\varepsilon^{\eta,\,\varphi}_{\alpha\beta}$ where our xenon-based DD experiments lose sensitivity. For the diagonal elements, this region simply follows from Eq. 69 and, for the off-diagonal elements, the more complicated behaviour is expressed in Eq. 70. Furthermore, the off-diagonal elements again exhibit some non-trivial behaviour in the form of ‘kinks’ (see for example $\varepsilon^{\eta,\,\varphi}_{e\mu}\approx 1.0$ and $\varphi\approx-5/16\,\pi$). The appearance of these kinks is analogous to those observed in the NR sensitivities in the nucleon plane, as described at the end of Section IV.2.

We re-iterate that the limits presented in Figs. 4 and 5 have been calculated by switching on only one NSI parameter at a time. Due to potential interference effects between different NSI parameters, a global analysis that allows all NSI parameters to vary, before marginalising to compute the limits on any one parameter, would generally lead to weaker bounds Dutta:2020che ; Coloma:2022avw . However, the point of our study is to illustrate the potential of DD experiments in this direction. Our study makes a strong case for their inclusion in future global analyses.

IV.4 Final Remarks

Having access to both the NR and ER signals in one experiment makes DD incredibly powerful as a probe for NSI. As far as we are aware, this is the only experimental technology that is able to perform such simultaneous analyses. For example, if a signal inconsistent with the SM is detected in the future, both channels would be pivotal for exploring the possible values of $\eta$ and $\varphi$, or, equivalently, the relative strength of NSI with electrons, protons and neutrons. This will come in tandem with other more traditional searches for new physics in the neutrino sector. However, given the number of parameters one is trying to constrain or fit, the addition of DD will provide important input complementary to that of oscillation and spallation source experiments.

Above, we have treated the NR and ER signals in DD experiments as separable. Indeed, in the name of background discrimination for DM searches, DD experiments are capable of this for large parts of the signal region. Taking into account experimental inputs, as described in Eq. 64, we are able to model NR and ER spectra accurately without resorting to a full Monte Carlo simulation of the detector responses in terms of S1 (scintillation) and S2 (ionisation) signals. Since detector responses from the point of interaction, be it NR or ER, have been well studied and calibrated within experimental collaborations, we are confident that introducing nonzero NSI will not alter the expectation that future experiments will be able to resolve S1 and S2 signals.

Interestingly, many DD collaborations also perform S2-only analyses, which has the benefit of lowering the experimental threshold $E_{\rm th}$, increasing the sensitivity to lighter DM values. However, this comes at the cost of losing NR/ER discrimination. In our NR projections for future experiments, we took the liberty of lowering $E_{\rm th}$ in a modest way, assuming that NR/ER discrimination was still possible, and indeed S2-only analyses boast much lower thresholds. As this choice implies, reducing $E_{\rm th}$ is beneficial for the NR signal, but not necessarily for the ER signal. This is because low-energy neutrinos are unable to impart sufficient energy to excite the bound electrons. It is likely then that an S2-only analysis will only improve the prospects of the NR signals, but one would then have to account for larger background rates.

Furthermore, in our analysis of DD experiments, we have not included argon-based liquid experiments. This direction is not without its potential, but we leave incorporation of such experiments for future work. As can be seen in Ref. Amaral:2020tga , the prospects for argon detectors are not as promising as those for xenon. However, Ref. Amaral:2020tga only considered the implications for specific BSM scenarios. The limiting factor for argon detectors seemed to be the experimental threshold and increased ER backgrounds, both of which tend to be much higher than their xenon counterparts. Recent progress from the DarkSide collaboration indicates that argon detectors may be able to provide competitive bounds in the future DarkSide:2018kuk ; DarkSide:2022dhx ; DarkSide-50:2022qzh .

To our knowledge, this work is the first to derive dedicated limits on NSI from both CE$\nu$NS and E$\nu$ES in DD experiments from the recent first results XENONnT XENONCollaboration:2022kmb and LZ LZ:2022ufs . Our analysis of future multi-ton LXe detectors in this section has exposed the huge potential that DD experiments have in fully exploring the parameter space of NSI, especially due to their increased sensitivity to E$\nu$ES. While there have been some initial studies considering the E$\nu$ES signals for non-zero NSI Dutta:2017nht ; Dutta:2020che , we take a comprehensive approach by considering the solar CE$\nu$NS and E$\nu$ES signals, modelling the experimental setups as close to the experimental collaborations as possible, and treating solar neutrino propagation in the coherent density matrix approach. Combined with our convenient parametrisation of the NSI parameter space, this allows us to derive an accurate overview of DD sensitivities and blind spots in the entire NSI parameter space. Moreover, as we pointed out in the previous section, since the blind spots for CE$\nu$NS and E$\nu$ES do not coincide in the charged NSI plane, DD experiments using a combination of both signatures can effectively avoid these and remain sensitive for most of this region.

When presenting our extended parametrisation in Section II, we introduced the axial-vector NSI coupling $\tilde{\varepsilon}_{\alpha\beta}^{f}$, only to set it to zero because it does not contribute to matter effects. Similarly, the effect of $\tilde{\varepsilon}_{\alpha\beta}^{f}$ on NRs will be minimal because of the coherent nucleon number enhancement that the vector current receives over the axial-vector. This enhancement is no longer present when one considers ERs, and would constitute an additional set of parameters that one could probe. The sensitivity of DD experiments to such axial-vector NSI is an interesting direction to be studied in future work.

Finally, we comment on the particle physics interpretation of the most promising projections we report in this work. Namely, the potential for a DARWIN-like experiment to probe NSI at the level of $\varepsilon^{\eta,\,\varphi}_{ee}\sim 10^{-2}$. Reinterpreting this value in a more canonical EFT approach implies $\Lambda_{\rm NP}/\sqrt{C}\sim 1.8\,{\rm TeV}$, where $C$ is the Wilson coefficient of the four-fermion operator. We see here that for $C>1$ these experiments have the capability to probe new physics above the TeV scale. In this context, it would presumably make the most sense to embed NSI analyses within the more general SMEFT framework. It would be interesting to study whether this can already be done in a consistent way. SMEFT observables are typically at collider scales, while NSI studies are much below, so the effective approach is appropriate for a greater range of $\Lambda_{\rm NP}$.

V Conclusions

We have demonstrated that direct detection experiments will soon become powerful probes of neutrino non-standard interactions, testing the parameter space in a complementary way to spallation source and oscillation experiments. This owes to their simultaneous sensitivity to nuclear and electron recoils and their unique capability to test tau neutrinos from the solar neutrino flux.

To do so, we have developed an extension of an earlier NSI parametrisation, allowing for non-standard interactions with nucleons and electrons simultaneously. Our parametrisation captures the rich phenomenology that arises when one allows for NSI to impact both neutrino propagation and neutrino scattering. We have shown that previous NSI limits from spallation source experiments, such as CENNS-10 LAr, and oscillation experiments, such as Borexino, map non-trivially to this extended parameter space, demonstrating the importance of allowing for a variable NSI contribution from the proton and the electron.

We have derived current direct detection constraints and projected the sensitivities of future direct detection experiments on the NSI landscape using the expected solar neutrino rate. We have thoroughly studied the resulting bounds in the different NSI directions by taking into account both CE$\nu$NS and E$\nu$ES signals, accurately modelling the experimental setups and consistently treating the coherent neutrino propagation via the density matrix. Furthermore, we have identified the potential blind spots where sensitivity is lost due to cancellations in the expected rate.

While current leading constraints from LZ and XENONnT are not competitive in this landscape yet, we have shown that those from future experimental runs and the projected DARWIN detector will cut into new regions of the NSI parameter space. We believe that the conclusion is clear: upcoming multi-ton, LXe-based DD experiments are poised to make a considerable impact in the neutrino NSI landscape. We therefore recommend that they be included in future global NSI studies, incorporating a more complete treatment of the systematics.

Acknowledgements

We want to thank Michele Maltoni for many insightful discussions on neutrino oscillations and non-standard interactions. We also thank Felix Kahlhoefer for the invaluable help with the LZ implementation, as well as Christropher Tunnell and Aarón Higuera for valuable discussions regarding XENONnT. Finally, we are also grateful to Pilar Coloma, Enrique Fernández Martínez, Danny Marfatia, Iván Martínez-Soler, Pablo Martínez-Miravé, Yuber F. Pérez-González and Thomas Schwetz for helpful discussions during the preparation of this manuscript.

DA is supported by the National Science Foundation under award 2209444. DGC acknowledges support from the Spanish Ministerio de Universidades under grant SI2/PBG/2020-00005. AC is supported by the grant “AstroCeNT: Particle Astrophysics Science and Technology Centre” carried out within the International Research Agendas programme of the Foundation for Polish Science financed by the European Union under the European Regional Development Fund. PF would like to express special thanks to the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149), for its hospitality and support. The work of PF was partially supported by the UKRI Future Leaders Fellowship DARKMAP. This work is partially supported by the Spanish Agencia Estatal de Investigación through the grants PID2021-125331NB-I00 and CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033.

APPENDICES

Appendix A Solar neutrino transition rate

In writing Eq. 1, we automatically retain the full phase correlation of the different solar neutrino flavour states reaching the detector. The way to understand how this formula comes about is to consider the amplitude for the combined propagation $\nu_{\alpha}\to\nu_{\gamma}$ of solar neutrinos from the point of production to the detector, and the scattering $\nu_{\gamma}\,T\to\nu_{f}\,T$ of the propagated neutrino $\nu_{\gamma}$ with the target material $T$ into any final state neutrino $\nu_{f}$, which is then lost. In quantum field theory, asymptotic in- and out-states are definite momentum states of the free theory, and as such they are on mass-shell. Hence, they describe well-defined mass eigenstates in the far past and far future. For our scattering of the neutrino $\nu_{\gamma}$ on the target $T$, we thus have to sum over all possible mass eigenstates $\nu_{i}$ in the out-state.111111We would like to express our special thanks to Thomas Schwetz for clarifying discussions on this point. Hence, the amplitude square of the process $\nu_{\alpha}\to\sum_{i}\nu_{i}$ can be written as

[TABLE]

where we have factored out the nuclear part of the elastic scattering process. Here, $S$ is the $S$-matrix for the full process, and we have decomposed the neutrino mass eigenstate $\nu_{i}$ into its flavour components, $|\nu_{i}\rangle=\sum_{\beta}U_{\beta i}|\nu_{\beta}\rangle$, in terms of the PMNS matrix $U_{\beta i}$. Then, the amplitude can be decomposed as

[TABLE]

where in the second to last line we have used the unitarity of the PMNS matrix, $\sum_{i}U^{*}_{\beta i}U_{\lambda i}=\delta_{\beta\lambda}$. Here, $\pi^{(\alpha)}$ is the projector onto the neutrino-flavour state $|\nu_{\alpha}\rangle$. In the second line, we have separated the $S$-matrix into $S_{\mathrm{prop}}$, describing the propagation of the initial neutrino $\nu_{\alpha}$ from the source to the detector, and $S_{\mathrm{int}}$, describing the interaction with the detector material. Thus, decorating the expression Eq. 76 with the relevant phase-space factors for the generalised cross section (cf. Eq. 45), we finally find that

[TABLE]

Appendix B $\Delta\chi^{2}$ Plots for CENNS-10 LAr and Borexino