In the Donnelly-Walecka theory [67] all semi-leptonic nuclear processes at low and intermediate energies may be described by an effective interaction Hamiltonian through the leptonic jμlept and hadronic Jμ current densities,

where G = G_F is the Fermi coupling constant for neutral-current processes and G = G_F cos θ_c (θ_c is the Cabbibo angle) for charged-current processes. For partial scattering rates, the evaluation of the transition amplitudes 〈f | Ĥ_eff | i〉 are treated via a multipole decomposition analysis of the hadronic current (see the Appendix). Then, for a given set of an initial |J_i〉 and a final |J_f〉 nuclear state, the double differential SM cross section becomes [68]

with ε_f (|k_f|) denoting the final energy (momentum) of the outgoing lepton, while ω = ε_i − ε_f stands for the excitation energy of the nucleus where ε_i is the initial lepton energy. For charged-current processes, the Fermi function F(Z, ε_f), takes into account the final state interaction of the outgoing charged particle, while for neutral-current processes such as coherent and incoherent neutrino-nucleus scattering it is F(Z, ε_f) = 1.

The individual cross sections in Equation (2) receive contributions from the so-called Coulomb M^, longitudinal L^, transverse electric Tel^ and transverse magnetic Tmag^ operators for both vector and axial vector components (see the Appendix). The cross sections σCLJ and σTJ are expressed in terms of the reduced matrix elements of the eight basic irreducible tensor operators [67]

Here, the + (−) sign refers to neutrino (antineutrino) scattering and θ represents the scattering angle, while the parameters a, b, d are expressed as

The 4-momentum transfer is trivially obtained from the kinematics of the process and in natural units reads

while for later convenience the magnitude of the 3-momentum transfer is defined as

For sufficiently small momentum transfer, i.e., q ≤ 1/R where R is the inverse nuclear radius¹, CEνNS dominates (see Figure 1). In this case, only ground state to ground state (g.s. → g.s.) transitions occur and lead to the following simplifications: the kinematics of the reaction imply m_f = 0 and |k_f| = ε_f so that a = 1 and d = 0, while the momentum transfer can be cast in terms of the incoming neutrino energy E_ν in the simple form

where the usual notation ε_i = ε_f ≡ E_ν has been adopted. Note also that the excitation energy in this case is ω = 0 and κ=-q2=Q2, while angular momentum conservation implies that for CEνNS processes the only non-vanishing operator is the Coulomb, T10≡M^00 (see the Appendix for the definition of the operators TiJ). Then the corresponding differential cross section is further simplified and takes the form

where the matrix element for g.s. → g.s. transitions is explicitly written in terms of the nuclear form factors for protons Fp(Q2) and neutrons Fn(Q2), as

At CEνNS experiments the detection mechanism is sensitive to the tiny nuclear recoils generated in the aftermath of the scattering process. It is therefore reasonable to express the differential cross section with respect to the nuclear recoil energy T_N, which in the low energy approximation T_N ≪ E_ν, reads

where TN=Q2/2M and M is the mass of the nuclear isotope. The calculations of Papoulias and Kosmas [61] involved the BCS form factors for protons (neutrons)

with N_n = Z or N and vp(n)j represents the occupation probability amplitude of the j-th single nucleon level.

The method described above, involves realistic nuclear structure calculations making it more reliable compared to the use of phenomenological form factors, especially for accelerator neutrino sources (see the discussion in section 2.2). For the reader's convenience, Equation (11) is also expressed through the vector weak nuclear charge QWV in the approximation of equal proton and neutron form factors, as [69]

where the vector QWV weak charge is given by [70]

with the left- and right-handed couplings of u and d quarks to the Z-boson being

The latter expressions include the radiative corrections from the PDG [71]: ρνNNC=1.0082, κ^νN=0.9972, λ^{u, L} = −0.0031, λ^{d, L} = −0.0025 and λ^{d, R} = 2λ^{u, R} = 3.7 ×10⁻⁵ while concerning the weak mixing-angle the adopted value is ŝZ2≡sin2θW=0.2382. Regarding the incoherent neutrino-nucleus cross section and for the sake of completeness we note that apart from the Donnelly-Walecka method given in Equation (2) a usefull formalism has been recently given in Bednyakov and Naumov [62].

The differential cross sections dσ/dT_N and dσ/dcosθ are shown in the upper left and upper right panel of Figure 2, from where it can be seen that large differences appear if the form factor dependence is neglected. On the other hand at low neutrino energies, i.e., E_ν ≤ 20MeV (relevant for reactor and solar neutrinos), the agreement of these two approximations is rather good. It is worth mentioning that forward scattering (θ = 0) leads to maximum dσ/dcosθ, as well as that for this particular case the form factor is by definition equal to unity due to the zero momentum transfer, see Equation (8). Finally the bottom panel illustrates a comparison of the CEνNS cross section by incorporating the nuclear form factors and by assuming F = 1.

Electron scattering data provide high precision measurements of the proton charge density distribution [73]. The absence of similar data for neutron densities, restricts us to rely on the approximation of ρ_p(r) = ρ_n(r) and thus assume F_p(Q) = F_n(Q) ≡ F(Q) (see Equation 13). In the context of nuclear theory, it is possible to treat separately the proton and neutron nuclear form factors by employing non-trivial techniques. The most reliable methods for this purpose include the large-scale Shell-Model [74, 75], the QRPA [76], Microscopic Quasiparticle Phonon Model (MQPM) [77], the deformed Shell-Model (DSM) method [39] and others. Recently, crucial information on important nuclear parameters has been extracted from the analysis of the recent COHERENT data in Cadeddu et al. [32], Ciuffoli et al. [33], Huang and Chen [34], and Papoulias [78].

The point-nucleon charge density distribution ρ(r), is defined as the expectation value of the density operator [79]

where the + (−) sign refers to the point-proton (neutron) charge density distribution. Assuming the nuclear ground state to be approximately described by a Slater determinant constructed from single-particle wavefunctions, the distributions of Equation (16) are given by summing in quadrature the point-nucleon wavefunctions. According to Kosmas and Vergados [79], for closed (sub)shell nuclei the charge density distribution is assumed to be spherically symmetric while the interesting radial component (r = |r|) of the proton charge density distribution, ρ_p(r), can be cast in the form

where R_nlj(r) denotes the radial component of the single-particle wavefunction with quantum numbers n, l, and j. The nuclear form factor depends on the three momentum transfer squared q² ≡ |q|² and can be obtained via a Fourier transformation

where j₀(x) = sinx/x denotes the zero-order Spherical Bessel function of first kind. The nuclear form factors lead to a suppression of the CEνNS cross section and subsequently to a suppression of the expected event rates (see [32] for a comparison with the COHERENT data). The uncertainties of the nuclear form factors are explored in Aristizabal Sierra et al. [35] where it is pointed out that studies looking for physics beyond the SM can be seriously affected by the uncertainty of the neutron form factor [80]. It is therefore important to treat with special care the nuclear form factors since new physics could be claimed or missed, if their uncertainties are not properly taken into account. In addition to the form factors obtained in the framework of the nuclear BCS method of Equation (12), below we present a summary of various form factor approximations widely considered in the literature.

The proton nuclear form factors Fp(q2), may be evaluated through a model independent analysis (e.g., by employing a Fourier-Bessel expansion) of the electron scattering data [66], having however the disadvantage of assuming Fp(q2) =Fn(q2).

For Harmonic Oscillator (h.o.) wavefunctions the nuclear form factor Fp(q2) for protons can be expressed in polynomial form [81, 82]

with N_max = (2n + l)_max denoting for the number of quanta of the highest occupied proton (neutron) level. In a similar manner, the radial nuclear charge density distribution ρ_p(r) is written in terms of the polynomials Π(r/b, Z) in the following compact form [81, 82]

with the definition χ = r/b (b stands for the h.o. size parameter). The explicit expressions for calculating the coefficients θ_λ and f_λ are given in the Appendix.

The occupation probabilities entering Equations (18) and (19) are assumed equal to unity (zero) for the states below (above) the Fermi surface, e.g., the filling numbers of the states for closed (sub)shell nuclei are those predicted by the simple Shell-Model. Going one step further, Kosmas and Vergados [79] introduced depletion/occupation numbers to describe the occupation probabilities of the surface levels, which satisfy the relation

In this framework, there is a number of active surface nucleons (above or below the Fermi level) with non-vanishing occupation probability 0 ≤ α_nlj ≤ 1 and a number of core levels with α_nlj = 1. The parameters are properly adjusted so that a high reproducibility of the experimental data is achieved [66]. By introducing four parameters α_i, i = 1, 2, 3, 4 in Equation (21) the polynomial Π(χ, Z) of Equation (20) reads

with λ = α₁ + α₂ − α₃ − α₄ (see [61] for the fitted values).

Besides calculations in the spirit of a nuclear structure model, a reliable description of the nuclear form factors (at least for low-energy reactor and solar neutrinos) may be obtained through the use of phenomenological approximations of the charge density distribution. The Helm-type density distribution is a convolution of a uniform nucleonic density with cut-off radius R₀ (accounting for the interior density) with a Gaussian falloff with folding width s (surface thickness). The corresponding Helm form factor takes the analytical form as [83]

where j1(x)=sinxx2-cosxx is the 1st-order Spherical Bessel function. The first three moments can be analytically expressed as [84]

The surface thickness parameter is fixed to s = 0.9 by fitting muon spectroscopy data [85], having also the advantage of improving the matching between the Helm and the symmetrized Fermi (SF) distributions [86]. The SF approximation follows from a Woods-Saxon charge density distribution and is expressed through the half density radius c and the diffuseness parameter a, as [87]

with

while the surface thickness is written as t = 4aln 3 [32]. The corresponding first three moments of the SF form factor read [84]

The Klein-Nystrand (KN) distribution is obtained from the convolution of a Yukawa potential with range a_k = 0.7 fm over a Woods-Saxon distribution (hard sphere with radius R_A). The resulting KN form factor reads [88]

and is adopted by the COHERENT Collaboration, while in this case root mean square (rms) radius reads

Figure 3 presents the charge density distribution in the top panel and the corresponding nuclear form factors for ⁴⁰Ar (interesting for LAr CEνNS detectors) and ⁴⁸Ti (interesting for μ⁻ → e⁻ conversion in nuclei) in the lower panel, while the results are compared for the various methods used. A comparison of the form factors for ¹²⁷I and ¹³³Cs that are of interest for COHERENT, evaluated with the DSM method (not covered here), with those of the Helm, SF and KN parametrizations, is given in Papoulias et al. [36]. By incorporating realistic nuclear structure calculations on the basis of the BCS method, the SM CEνNS cross section is given in Figure 4 for a set of different isotopes throughout the periodic table. For heavier isotopes the form factor suppression is more pronounced and therefore the cross sections flatten more quickly, since the nuclear effects become significant even at low-energies.

The left panel of Figure 7 shows the expected number of events at the CsI[Na] detector of COHERENT and gives a comparison with the experimental data, from where it can be seen that a good agreement is reached. In Papoulias and Kosmas [21] the authors analyzed the CEνNS data and obtained a low-energy determination of the weak mixing angle, as illustrated in the right panel of Figure 7. The obtained constraint at 90% C.L. reads [21]

An interesting analysis combining atomic parity violating (APV) and CEνNS data was performed in Cadeddu and Dordei [101], while the prospects regarding the future reactor-based CEνNS experiments such as those presented in Table 1, have been extracted in Cañas et al. [102]. On the other hand, an improved determination of the CsI[Na] quenching factor can in principle lead to a significantly better agreement between the experimental results and the theoretical simulations [103], as well as to an improved sensitivity on the weak mixing angle [78].

The discussion made in the previous section emphasized how the CEνNS cross section depends on the nuclear physics effects which are incorporated through the momentum variation of the relevant nuclear form factor. The authors of Aristizabal Sierra et al. [27] demonstrated how the intrinsic nuclear structure uncertainties may have a significant impact to searches beyond the SM such those regarding NSIs, sterile neutrinos and neutrino generalized interactions (GNIs). Starting from the form factor of Equation (18) and expanding in terms of even moments of the charge density distribution one arrives to a model independent expression [104]

with the k-th radial moment defined as

allowing the study of contributions of higher-order moments to nuclear form factors [33]. A sensitivity analysis of the two first moments with current and future COHERENT data is depicted in Figure 8 where the allowed regions are presented at 1σ, 90% and 99% C.L. The calculation in this case was restricted in the physical region [0,6] fm in order to obey the upper limit on R_n(²⁰⁸Pb) =5.75 ± 0.18 fm from the PREM experiment [105]. The future scenarios considered assume improved statistical uncertainties and more massive detectors in accord with the next generation COHERENT experiments [93] (see [36] for details), while as demonstrated in Patton et al. [104] multi-ton scale detectors will provide significant improvements.

The average CsI neutron rms radius has been explored in Cadeddu et al. [32], Huang and Chen [36], and Papoulias et al. [34] using the energy spectrum of the available CEνNS data. The corresponding sensitivity profiles are presented in Figure 9, leading to the best fits at 90% C.L. [36]

while the potential of improvement through a more accurate determination of the QF is promising (see e.g., [78]). An independent analysis combining APV and CEνNS data was performed in Cadeddu et al. [106] leading to essentially similar results. Finally, it is worthwhile to mention the reported upper bound on the neutron skin ΔRnp=ΔRn-ΔRp=0.7-1.1+0.9 fm [32].

Non-standard interactions (NSI) [107] appear in several appealing SM extensions [108] involving four-fermion contact interaction, various seesaw realizations [109–111], left-right symmetry [112], gluonic operators [113], etc., constituting an interesting model independent probe of new physics. NSIs may have implications to SN [114], neutrino oscillations [57] and CEνNS [70, 115], while recently NSI terms were explored in the context of GNI [116] and effective field theory (EFT) operators [117, 118]. Finally the RG issue has been partly addressed in the context of NSI in Davidson and Gorbahn [119].

For sufficiently low energies vector-type NSIs arise from the effective four-fermion operators [56]

leading to new contributions to the CEνNS rate from exotic processes of the form

where α, β = {e, μ, τ} (α ≠ β), q denotes a first-generation quark q = {u, d} and P = {L, R} is the chiral projection operator. For the case of CEνNS the new interactions are taken into account through the NSI charge with the substitution QWV→QNSIV in Equation (14). The latter contains flavor-preserving non-universal (ϵααqV) and flavor changing (ϵαβqV) terms and is expressed as

implying that the NSI CEνNS cross section becomes flavor dependent.

There is a reach literature on NSI investigations with the recent COHERENT data. Assuming one nonvanishing coupling at a time, the authors of Papoulias and Kosmas [21] focused on the non-universal terms and obtained the sensitivity profiles shown in the left panel of Figure 10, while the corresponding allowed regions resulting from a combined analysis of the NSI couplings are illustrated in the upper panel of Figure 11 at 90% C.L. Regarding the future prospects of MINER, Ricochet, NUCLEUS and CONNIE, similar studies were conducted concentrating on the non-universal [120] and flavor-changing [121] terms. Indeed, a multitarget strategy can break degeneracies involved between up and down flavor-diagonal NSI terms that survives analysis of neutrino oscillation experiments [14]. Constraints on the corresponding parameters arising from leptoquarks [120], GNI [116], and EFT [117, 118] have been also reported. NSI constraints from CEνNS place meaningful constraints excluding a large part of the existing CHARM constraints and overlap with results coming out of LHC monojet searches (see [120] for a usefull comparison). Regarding the near-term future, a potential improvement on determination of the QF [103] may yield severe constraints. For example, updated bounds are possible by analyzing the number of events [78], the energy spectrum [122] as well as through a combined analysis of both energy and timing COHERENT data [123].

Novel tensor-type interactions are predicted in the general context of NSI [124] and GNI [116] which induce terms of the form

Such interactions violate the chirality constraint allowing for a wide class of new interactions, e.g., relevant to neutrino EM properties (see [125, 126]). Contrary to the vector NSI case, for tensorial interactions there is absence of interference with the SM interactions. In the approximation of a vector-type translation the corresponding tensor NSI charge has been expressed as [124]

while a more systematic interpretation has been carried out in Aristizabal Sierra et al. [116].

To account for the new contributions in the presence of tensorial NSI, the CEνNS cross is written [21]

with the tensor NSI factor defined as

From the analysis of the COHERENT data, the sensitivity profiles accounting to tensor NSIs, assuming one non-zero coupling at a time, are illustrated in the right panel of Figure 10 (see also [21]). The corresponding allowed regions coming out from a two parameter analysis are presented in the lower panel of Figure 11 at 90% C.L. The result is more stringent as compared to the analysis carried out in the framework of GNI for reasons discussed above. On the other hand, comparing with the vector NSI case the absence of SM-tensor NSI interference causes the allowed regions to appear with more narrow bands.

Theories beyond the SM with an additional U(1)′ symmetry have been comprehensively investigated. Regarding CEνNS related studies a novel massive mediator predicted in these concepts is expected to induce a detectable distortion to the nuclear recoil spectrum, provided that its mass is comparable to the momentum transfer. The study of models with new vector or scalar interactions that involve hidden sector particles may be also accessible at CEνNS experiments [127]. Such frameworks are interesting since they may play a central role in explaining anomalies with regards to B-meson decays at the LHCb experiment [128] and at DM searches [129].

We first examine the case of a new massive vector boson Z′. Restricting ourselves to the neutrino sector with only left-handed neutrinos the Lagrangian reads [130]

The arising cross sections imply a re-scaling of the SM one according to the expression

with the Z′ factor taking the form

Here, gZ′νV denotes the neutrino-vector coupling, while the respective Z′ charge reads [129]

However, in the general case the ν − Z′ coupling is flavor dependent (gZ′νV)αβ. Denton et al. [16] has explored this possibility and concluded that for a sufficiently small momentum transfer with respect to MZ′, the new physics contributions can be addressed in the form of NSIs

where the Z′ has been integrated out. Unlike the NSI case that can only modify the energy spectrum by a global factor, the additional momentum dependence expected due to the new light mediators discussed here can be well encoded in the detected signature and subsequently lead to conclusive indications of the new physics nature.

We now turn our attention on new interactions induced by the presence of a CP-even mediator. In particular, we consider a new real scalar boson ϕ with mass M_ϕ, based on the Lagrangian [130]

with gϕqS and gϕνS representing the scalar-quark and scalar-neutrino couplings, respectively. In this framework the SM CEνNS cross section acquires an additive contribution due to the boson exchange that can be quantified in terms of the respective cross section

with the scalar factor Gϕ being

Analogously to the previous case, the corresponding scalar charge is defined as [130]

where the form factors fT,q(N) capture the effective low-energy coupling of ϕ to the nucleon N={p,n} (mN is the nucleon mass) for the quark q.

As discussed previously, different new physics signatures are expected to leave different imprints on the event and recoil spectrum. Contrary to the Z′ scenario discussed above, the Dirac structure of the ϕνν̄ vertex accounting for the scalar mediator is different (chirality-flipping) with respect to the SM one (chirality-conserving). Indeed, there is no interference between vector (or axial-vector) neutrino interactions and (pseudo-)scalar, tensor neutrino interactions [29, 131]. Therefore, the absence of interference between SM-scalar interactions gives rise to an overall modification of the expected CEνNS spectrum (see [129]). Moreover, comparing the vector and scalar cross sections it becomes evident that the corresponding signals are expected to be well distinguishable. The scalar effects are not pronounced at eV-thresholds, while on the contrary they are expected to be stronger at recoil energies of the order of keV. A thorough classification of the new physics signatures with respect to vector and scalar interactions is given in Aristizabal Sierra et al. [132] providing also key information on the possibility of breaking isospin-related degeneracies from combined measurements with different detector material.

Assuming universal couplings, one finds the equalities [130]

Using the COHERENT data, bounds have been put on the parameter planes (gZ′2,MZ′) and (gϕ2,Mϕ) for the vector and scalar mediators, respectively [21]. In the left panel of Figure 12, the limits are shown at 90% C.L., where a degenerate area appears that cannot be excluded by the data is found due to the cancellations involved in Equation (44). However as shown in Liao and Marfatia [13] this degeneracy can be broken in the context of NSI, while for heavy mediator masses, MZ′≫2MTN ~50 MeV, it remains unbroken and depends on the ratio

For the case of light mediator masses MZ′≪2MTN, it holds

which implies that the bound is only sensitive to the coupling. The latter could be drastically improved by combining data from different detectors [133]. The case of a scalar field mediating CEνNS is explored in the right panel of Figure 12 where the extracted bounds in the (Mϕ,gϕ2) plane are depicted at 90% C.L. Significant improvements are possible through powerful analyses that are based on both energy and timing COHERENT data [17] as well as by taking into account improved quenching factors [78]. The future of CEνNS experiments will offer a complementary probe to various existing limits in the low- and high-energy regime. The currently best results for a low-energy light vector mediator of MZ′<10 MeV have been recently reported by the CONNIE Collaboration [134]. The attainable sensitivity is expected to be competitive with existing bounds from neutrino-electron scattering, dark photon searches at BaBar and LHCb results (see [30, 120]). Before closing our discussion it is important to note that, very recently CP violating effects have been also analyzed with the current and future COHERENT data in the context of light vector mediator scenarios [27]. The latter have been also found to be applicable to reactor or solar/atmospheric neutrino searches with important implications on multi-ton dark matter detectors.

Non-trivial electromagnetic (EM) properties of massive neutrinos constitute an interesting probe to look for physics beyond the SM and at the same time they are crucial for distinguishing between the Dirac or Majorana nature of neutrinos [135]. The two main phenomenological parameters observable at a neutrino experiment are the effective neutrino magnetic moment and the neutrino charge radius (the possibility of a neutrino millicharge is explored in [23, 106]). Assuming Majorana neutrinos, the EM neutrino vertex is described by the electromagnetic field tensor F_αβ of the effective Hamiltonian [3, 136].

while for the case of Dirac neutrinos one has

It is important to note that for Majorana (Dirac) neutrinos λ = μ−iϵ is an antisymmetric complex (general complex) matrix. The two imaginary matrices μ (magnetic moment) and ϵ (electric dipole moment) obey the respective properties μ^T = −μ (μ = μ^†) while ϵ^T = −ϵ (ϵ = ϵ^†). It thus becomes evident that, unlike the Dirac case, for Majorana neutrinos the diagonal moments are vanishing μiiM= ϵiiM = 0.

For a low-energy neutrino scattering experiment the observable neutrino magnetic moment (flavor dependent) is in fact a combination of the neutrino transition magnetic moments (TMMs) discussed above. In the mass basis it reads [137]

In Equation (56) the following transformations have been introduced

where the 3−vectors 𝔞₊ and 𝔞₋ denote positive and negative helicity states, respectively, while the magnetic moment matrix λ (λ~) in the flavor (mass) basis is written as [138]

with Λα=|Λα|eiζα and Λi=|Λi|eiζi being the TMMs in the flavor and mass basis, respectively, where ζ_α and ζ_i denote the corresponding CP-phases.

The potential EM neutrino properties appear in the form of an effective neutrino magnetic moment that is conveniently expressed in the mass basis according to Equation (56) in terms of fundamental parameters (TMMs, CP-violating phases and neutrino mixing angles). The latter induce an additive contribution to the SM cross section [139]

where the EM factor reads [21]

Here, the EM charge QEM is written in terms of the fine structure constant a_EM and the effective neutrino magnetic moment, as [115]

Moreover, the effect of a neutrino charge radius can be taken into consideration through a shift in the definition of the weak mixing angle [140]

where by sin2θW¯ it is denoted the low energy value of the weak mixing angle, e.g., ŝZ2=0.2382.

The presence of a neutrino magnetic moment is expected to yield a distortion in the recoil spectrum during the CEνNS process, i.e., a detectable excess of events especially for low recoil energies. The left panel of Figure 13 shows the χ² profile of the effective neutrino magnetic moment extracted by the first light of COHERENT data. A similar analysis has been performed in order to quantify the sensitivity of COHERENT on the neutrino charge radius as shown in the right panel of the same figure. Note that, an essential improvement due to a more accurate treatment of the QF is possible (see [78] for more details).

The authors of Miranda et al. [22] performed a comprehensive analysis on the sensitivity of various existing and future CEνNS experiments and extracted constraints on the different components Λ_i of the neutrino TMM matrix. In particular, their study focused on existing and next generation experimental setups of COHERENT as well as on the expected data from the future reactor-based experiments: CONUS, CONNIE, TEXONO, MINER, and RED100. In a similar manner, Billard et al. [120] extracted constraints focusing on the NUCLEUS and Ricochet detectors at the Chooz NPP, however assuming the effective neutrino magnetic moment. Miranda et al. [22] performed a systematic combined analysis with regards to the TMMs exploring also the effects of the CP violating phases of the complex matrix given in Equation (58). As a concrete example, Figure 14 shows the contours in the (|Λ_i|–|Λ_j|) parameter plane for the case of current and next generation reactor-based CEνNS experiments. It is worth mentioning that these bounds are comparable to existing ones from low energy solar neutrino data at Borexino phase-II [141]. Figure 15 shows the sensitivity contours in the (〈rνα2〉,〈rνβ2〉) plane that resulted from the COHERENT data. A similar analysis is performed in Papoulias [78] and Khan and Rodejohann [122], while for a comprehensive fit including energy and timing data the reader is referred to Cadeddu et al. [106].

The three-neutrino paradigm has been put in rather solid grounds from the interpretation of solar and atmospheric oscillation data. On the other hand, controversial anomalies such as those coming from recent reactor data as well as existing anomalies implied by the LSND and MiniBooNE experiments inspired a reach phenomenology beyond the three-neutrino oscillation picture, based on the existence of a fourth sterile neutrino state with eV-scale mass (m_{1, 2, 3} ≪ m₄) and tiny mixing angles. To accommodate sterile neutrinos the lepton mixing matrix is minimally extended so that the flavor eigenstates ν_α, (α = e, μ, τ, s, ⋯ ) are related to the mass eigenstates ν_i (i = 1, 2, 3, 4, ⋯ ) by the unitary transformation να=∑iUαiνi. Then, for the short-baseline CEνNS experiments the survival probability of an active neutrino at a distance L is written [26]

with the abbreviation Δ≡Δm2L/4Eν and the mass splittings under the approximation Δm412≈Δm422≈Δm432≡Δm2. At this point it should be stressed that neutrino-nucleus scattering experiments are favorable facilities to probe sterile neutrinos being complementary to dedicated experiments such as MINOS/MINOS+ [142], MiniBooNE [143], Daya-Bay [144], Juno [145], and NEOS [146]. Indeed, due to the purely neutral-current character of the CEνNS process it is not necessary to disentangle between active-sterile neutrino mixing [147].

The possibility of investigating sterile neutrinos in the simplest (3+1) scheme through the CEνNS process was examined for the first time in Anderson et al. [148], relying on an SNS source. A combined sterile neutrino analysis was performed in Kosmas et al. [24] highlighting the complementarity between accelerator and reactor neutrino sources, by focusing on COHERENT and TEXONO experiments, respectively (see Figure 16). Moreover, a detailed study of various reactor-based CEνNS proposals has been carried out in Cañas et al. [25], showing how such future measurements can be exploited to solve the reactor antineutrino anomaly. After the first observation of CEνNS by the COHERENT experiment, Papoulias and Kosmas [21] reported the first constraints under the assumption of a universal new mixing angle, extracting the conclusion that the current sensitivity is rather poor. By exploiting timing data the potential of a future measurement at the next generation of COHERENT with a 100 kg CsI detector has been demonstrated in Blanco et al. [26], concluding that the prospects of probing the exclusion regions in the (Δm412–sin22θeμ) plane from the latest MiniBooNE [143] and LSND [149] are promising. Finally, focusing at CONUS in Berryman [150] it was shown that the complementarity between terrestrial-cosmological experiments may resolve the tension raised by astrophysical observations regarding the existence of sterile neutrinos.

Emphasis has been put on the physics beyond the SM by devoting a great part to the past and current research efforts and by concentrating on the various channels contributing to CEνNS processes and their interpretation. Through a χ² sensitivity analysis, based on the recoil or timing spectra of the COHERENT data, the current limits are listed at 90% C.L. in Table 2. For a given parameter set S, the best fit is found through the minimum value χmin2(S). The limits involve electroweak (weak-mixing angle), nuclear (nuclear radius), and physics beyond the SM (NSIs and EM neutrino properties). Significant improvements are expected through a more accurate determination of the QF and from a better control of the systematic uncertainties. The reported constraints on electroweak and NSIs have been extracted with various analysis methods, i.e., by combining existing APV measurements or global oscillation constraints with the recent COHERENT data, emphasizing the complementarity of CEνNS data in the low energy regime.