Coherent elastic neutrino-nucleus scattering: EFT analysis and nuclear responses

The cross section for coherent elastic neutrino-nucleus scattering (CE$\nu$NS) depends on the response of the target nucleus to the external current, in the Standard Model (SM) mediated by the exchange of a $Z$ boson. This is typically subsumed into an object called the weak form factor of the nucleus. Here, we provide results for this form factor calculated using the large-scale nuclear shell model for a wide range of nuclei of relevance for current CE$\nu$NS experiments, including cesium, iodine, argon, fluorine, sodium, germanium, and xenon. In addition, we provide the responses needed to capture the axial-vector part of the cross section, which does not scale coherently with the number of neutrons, but may become relevant for the SM prediction of CE$\nu$NS on target nuclei with nonzero spin. We then generalize the formalism allowing for contributions beyond the SM. In particular, we stress that in this case, even for vector and axial-vector operators, the standard weak form factor does not apply anymore, but needs to be replaced by the appropriate combination of the underlying nuclear structure factors. We provide the corresponding expressions for vector, axial-vector, but also (pseudo-)scalar, tensor, and dipole effective operators, including two-body-current effects as predicted from chiral effective field theory. Finally, we update the spin-dependent structure factors for dark matter scattering off nuclei according to our improved treatment of the axial-vector responses.

A crucial input in interpreting the measured cross section is the response of the nucleus. If BSM constraints are to be extracted, the nuclear structure has to be provided from elsewhere. In fact, since the weak charge of the proton is small, the SM CEνNS process mainly probes the nuclear neutron distribution, which is significantly less constrained experimentally than the electromagnetic charge distribution. Apart from CEνNS, the only direct information of the neutron distribution comes from parity-violating electron scattering (PVES) off lead [15,16]. Accordingly, assuming the absence of a significant BSM signal, the measured CEνNS cross section could be used to constrain the neutron distribution instead [17][18][19][20][21]. Currently, the nuclear input used in the interpretation of CEνNS experiments is mainly derived * E-mail: hoferichter@itp.unibe.ch † E-mail: menendez@fqa.ub.edu ‡ E-mail: schwenk@physik.tu-darmstadt.de from relativistic mean-field methods (RMF) [22,23], even though results based on nonrelativistic energy-density functionals are also available [24][25][26]. For argon, there is a recent first-principles calculation based on coupledcluster theory [27].
Here, we provide nuclear structure results for CEνNS, extending large-scale nuclear shell model calculations developed in the context of nuclear structure factors in direct-detection searches for dark matter [28][29][30][31][32][33][34][35][36]. First, the level of agreement between the shell-model, the RMF, and coupled-cluster results suggests that the form factor uncertainties are not as severe as claimed in Ref. [37], but in addition the shell-model approach also allows us to address the spin-dependent (SD) responses, which are similar, but somewhat different to the ones in SD dark matter searches.
To this end, we first derive the decomposition of the cross section into Wilson coefficients of effective operators, hadronic matrix elements, and nuclear structure factors. In the SM the effective operators just parameterize the Z-boson exchange, but this approach can be conveniently extended to include BSM effects. The hadronic matrix elements determine the hadronization of these operators at the single-nucleon level, and finally the nuclear structure factors take into account the many-body nuclear matrix element of these single-nucleon currents. As a first step, we demonstrate how the standard weak form factor emerges when combining all these ingredients into a single object. However, this analysis shows that even for the coherent part of the nuclear response four different underlying structure factors contribute to the cross sec-arXiv:2007.08529v1 [hep-ph] 16 Jul 2020 tion. Therefore, in principle the weak form factor needs to be modified as well when allowing for BSM effects in the Wilson coefficients, since their contribution does not factorize.
While for the dominant vector operators corrections beyond the single-nucleon currents are small [38,39], since the magnetic-moment form factors happen to be kinematically suppressed, this is no longer true for the axial-vector [29] or for scalar currents, see [32-34, 36, 40-45] for the analogous effects in the case of dark matter scattering off nuclei. As long as the SM dominates, such effects will only become relevant for CEνNS once experiments become sensitive to SD responses. Otherwise, mainly the limits on scalar operators would be affected, but in CEνNS such contributions are suppressed due to the need for right-handed neutrinos and lack of interference with the SM.
The paper is organized as follows: in Sec. II we first review the necessary formalism, regarding effective operators, hadronic matrix elements, and nuclear responses, with some details of the multipole expansion and nuclearstructure calculations summarized in the appendix. In Sec. III we first introduce the charge and weak form factors in the context of electron scattering, before discussing the application to CEνNS. In particular, we present an improved treatment of the axial-vector responses both for CEνNS and dark matter. In Sec. IV we discuss how the nuclear responses need to be adapted when considering SM extensions, before concluding with a summary in Sec. V.

A. Effective Lagrangians
As a first step, we review the operator basis for CEνNS [36,46] 1 L (5) = C Fν σ µν P L νF µν , L (6) = q C V qν γ µ P L νqγ µ q + C A qν γ µ P L νqγ µ γ 5 q + C T qν σ µν P L νqσ µν q , where we adopted the following conventions: neutrino indices are suppressed throughout, indicating that oscillation effects are usually negligible at the scale of CEνNS experiments, so that incoming and outgoing flavors are understood to be identical. In comparison to the case of a dark-matter spin-1/2 particle [36] the number of operators is reduced by a factor of 2 when assuming that the neutrino is left-handed. This is implemented in Eq. (1) in terms of projectors P L = (1 − γ 5 )/2, and given that observing chirality-violating effects would require right-handed neutrino beams (suppressed by tiny neutrino masses), we will not consider the opposite chirality in the following. With these conventions the set of operators is then similar to the dark-matter case: at dimension-5 level there is a single (dipole) operator involving the photon field strength tensor F µν . At dimension 6 we have the vector and axial-vector operators already present in the SM, as well as the tensor operator. Introducing quark masses for renormalization-group invariance, the scalar and pseudoscalar operators would be counted as dimension-7. The operator involving the QCD trace anomaly θ µ µ would also be of dimension 7. We have already rewritten the gluon term G a µν G µν a in terms of this operator (we will not consider the operators involving the dual field strength ten-sorG a µν or the photon field strength). In particular, we already integrated out the heavy quarks [47] and absorbed their effect into where C S g is the original coefficient of thē νP L ν α s G a µν G µν a gluon operator and we used the relation in the transition. Elsewhere, the sum over q in principle refers to all quark species, but in practice we will restrict the analysis to the light quarks q = u, d, s. Finally, there are operators with derivatives acting on the neutrino fields (in analogy to the spin-2 operator for dark matter), but we will concentrate on the more frequently considered operators in Eq. (1). We note that the dimensional counting is not unambiguous regarding the quark mass m q , e.g., sometimes the tensor operator is introduced at dimension-7 by adding a factor m q in this operator as well [46] (the notation in Eq. (1) follows Ref. [48]). Finally, we stress that the chirality-flipping operators, with scalar and tensor operators on the neutrino bilinear, require the presence of (final-state) righthanded neutrinos. In SMEFT [49] such operators are suppressed beyond dimension-6 level. In addition, the dipole operator leads to a new long-range interaction, and therefore C F is strongly constrained by astrophysical observations [50,51]. However, such operators have been suggested as a potential BSM explanation of the excess of electronic recoil events observed by XENON1T [52].
In the SM all Wilson coefficients except for C V q and C A q vanish, with Z exchange leading to the matching relations with Fermi constant G F = 1.1663787(6) × 10 −5 GeV −2 [53,54]. In the notation of Refs. [2,3,14], the deviations from these SM values are often expressed as where the sign of qA ee has been chosen in accordance with Ref. [14]. Finally, we can define dimensionless Wilson coefficientsC i = C i /Λ n , where Λ either corresponds to the respective BSM scale, or, in the SM, to the Higgs vev

B. Dimension-5 matrix elements
Having defined the operator basis (1), the second step concerns the nonperturbative input required to define amplitudes at the hadronic level. We will largely follow the conventions of Refs. [32,36], but for completeness review here the respective hadronic matrix elements.
For the dimension-5 operator only the electromagnetic form factors of the nucleon are required, without the need for a flavor decomposition. With N = {p, n}, we thus have the usual Dirac and Pauli form factors F 1 and F 2 , where j µ em = q=u,d,sq Q q γ µ q, Q = diag(2, −1, −1)/3, and q = p−p . For small momentum transfer t = (p −p) 2 , it is sufficient to consider the expansion around t = 0 with charge Q N , anomalous magnetic moment κ N , and in terms of the charge radius r 2 E N . We will use the numerical values given in Table I. In particular, we will use the proton charge radius from muonic atoms r 2 E p = 0.84087 (39) fm [57,58], in line with most recent electron spectroscopy measurements [61][62][63], the PRAD electron scattering data [64], and the expectation from dispersion relations [65,66]. For the neutron, we use the charge radius from Ref. [54], but note that a recent extraction from the deuteron points to a slightly smaller value [67].
C. Dimension-6 matrix elements At dimension 6 we first need the vector matrix elements for each quark flavor separately To perform the flavor decomposition, we will assume isospin symmetry (see Ref. [68] for corrections), which leads to Information on the strangeness form factors has traditionally been extracted from PVES, but the uncertainties are sizable [69]. More recently, it has been shown in lattice QCD that the strangeness contribution is very small, in Table I we quote the naive average of Refs. [70,71]. The second dimension-6 operator requires input on the axial-vector form factors, as they appear in the decomposition where, for completeness, we included the second-class current G q,N T (t) [72], but will not further consider its contribution in the following. The normalization is determined by the axial-vector charges Assuming again isospin symmetry these couplings are constrained by in terms of the axial-vector coupling of the nucleon g A = 1.27641(56) [73] and the SU (3) couplings D and F that can be extracted from semileptonic hyperon decays. In combination with the singlet combination from Ref. [74],   this leads to the couplings listed in Table I. These values are in reasonable agreement with lattice QCD [75] N f = 2 + 1 + 1 [76] : g u,p A = 0.777 (39), but in view of the present uncertainties we adopt the phenomenological determination. However, while part of the difference to phenomenology could be due to the scale dependence of the singlet combination, both lattice calculations point to a smaller strangeness coupling than extracted from the spin structure functions. The triplet and octet components of the induced pseudoscalar form factor G P (t) are constrained by Ward identities, whose manifestation at leading order in the chiral expansion becomes Finally, for the triplet component there is also experimental information on the momentum dependence. Defining the axial radius by a simple dipole ansatz with mass scale M A around 1 GeV [78], implies r 2 A = 12/M 2 A ∼ 0.47 fm 2 . The central value agrees with Ref. [79], a global analysis of muon capture and neutrino scattering, but the uncertainties are substantial, see Table I. To ensure that the Ward identity is satisfied up to higher orders, the pseudoscalar form factor needs to be modified according to [80] when including the radius corrections (17). The full πN coupling constant g πN N has been introduced as a convenient way to capture all chiral corrections at O(1). In the numerical analysis we will use F π = 92.28(9) MeV [54] and g 2 πN N /(4π) = 13.7(2) [81][82][83][84][85]. With this input, the Goldberger-Treiman discrepancy becomes demonstrating that the chiral corrections are rather small. The final matrix elements in the dimension-6 Lagrangian concern the tensor operatorqσ µν q, for which we use the decomposition The tensor charges g q,p T = F q,p 1,T (0), as given in Table I, are taken from lattice QCD [86]. The other form-factor normalizations come from Ref. [87] (with strangeness input updated to Ref. [71]).

D. Dimension-7 matrix elements
At dimension 7 we first need the scalar matrix elements of the nucleon To separate the momentum dependence we define The scalar couplings f N u,d are largely determined by the pion-nucleon σ-term σ πN [88], up to isospin-breaking corrections that can be extracted from the protonneutron mass difference [89][90][91][92][93]. The numbers given in Table I follow from σ πN as extracted from data on pionic atoms [82,83,[94][95][96] when used as input for a dispersive analysis of pion-nucleon scattering [97,98]. This result has been confirmed using independent input from scattering data [99], but there is a persistent tension with lattice QCD that still has not been resolved [100][101][102][103][104][105]. Accordingly, we have increased the error in f N s given that in this case all phenomenological extractions are subject to large SU (3) uncertainties. The heavy-quark couplings f Q effectively describe the matrix element of the trace anomaly at O(α s ) with normalization θ N 0 (0) = m N , while perturbative corrections especially for the charm quark lead to additional uncertainties [106]. For the momentum dependence,σ andσ s are taken from Refs. [33,107,108], and [75,109] which also determines the scalar couplings of the pion according to (14), (14).
These matrix elements arise in two-body corrections to scalar currents and are also included in Table I. Finally, we parameterize the pseudoscalar matrix element as Sij SN not coherent axial responses For the nonsinglet component the new form factor is related to the axial-vector ones by the Ward identity but in the singlet case this relation is broken by the anomaly contribution from G a µνG µν a , similarly to the gluonic contribution to the trace anomaly. For a consistent treatment of singlet effects one would thus have to extend the operator basis in Eq. (1) accordingly. In the past, the singlet pseudoscalar matrix element has often been estimated by assuming [110,111] N | q=u,d,sq but the analogous relation for the axial-vector singlet combination q=u,d,s ∆q does not display the expected 1/N c suppression. The matrix element of the gluon anomaly could be extracted with similar techniques as used for lattice calculations of the QCD θ term [112].

E. Nuclear responses
As a final step, the nucleon-level matrix elements need to be convolved with the nuclear states. Formally, the decomposition into distinct nuclear responses requires a multipole decomposition, see Refs. [113][114][115][116][117][118], which in full generality becomes very complex. Here, we concentrate on the most relevant contributions, with the main features summarized in Table II, and review some of the details needed later in App. A.
By far the most important response is related to the charge operator, it is denoted by the structure factors F M ± (q 2 ) normalized by This is the only response that is fully coherent. In addition to F M ± , we also need the so-called F Φ ± structure factor, which can be interpreted in terms of spin-orbit corrections. This response vanishes for q 2 = 0, but it interferes with F M ± and receives some coherent enhancement, especially for heavy nuclei. This is because in the relevant nuclei nucleons tend to occupy orbitals with spin parallel to the angular orbital momentum (lowered in energy by the nuclear spin-orbit interaction) and highenergy orbitals with antiparallel spin, which would cancel F Φ ± , remain mostly empty. The interference with F M and partial coherence make the Φ response the most relevant correction. In principle, both F M ± , F Φ ± may contribute beyond the leading L = 0 multipole, but such effects are not coherent, vanish at q 2 = 0, and without interference with the leading multipole effectively become negligible. Due to this we will continue to identify both responses with their L = 0 multipole.
Finally, there are several responses that emerge from the axial-vector operator. Their contribution again is not coherent, but remains finite at q 2 = 0. In these cases, several multipoles and responses become relevant, but we will continue to use a notation in which these effects are subsumed into structure factors S ij (with indices i, j = 0, 1 corresponding to isoscalar/isovector combinations). We keep the induced pseudoscalar form factors G q,N P , whose contribution is enhanced by the presence of the pion pole, but do not consider any other subleading noncoherent responses. Further aspects of the multipole decomposition are discussed in Sec. III whenever necessary to introduce the nuclear responses for a given process.

III. NUCLEAR RESPONSES IN THE STANDARD MODEL
In this section we will collect the nuclear responses as they appear in electron-nucleus and neutrino-nucleus scattering. In particular, we demonstrate how the traditional charge and weak form factors emerge in the formalism established in Sec. II. In either case, the kinematics are defined by and invariants Here, m A refers to the mass of the nucleus and lepton masses are neglected throughout. We also define P = p + p and write t = q 2 = −Q 2 .
A. Parity-conserving electron scattering For electron scattering, the invariants (34) are conventionally replaced in favor of where θ is the scattering angle in the laboratory frame. In this frame the relation of the spin-averaged scattering amplitude |M| 2 to the cross section becomes where the last relation defines the Mott cross section The incoming and outgoing electron energies are given by For the parity-conserving case, the amplitude becomes and at the single-nucleon level the hadronization follows from Eq. (6). The leptonic trace fulfilling k µ L µν = k µ L µν = 0, needs to be contracted with the nuclear amplitude, which we express in terms of multipoles according to Sec. II E, see Ref. [117] and App. A. The leading terms can be read off from the nonrelativistic expansion of the single-nucleon current where we dropped the remaining two-component spinors and the space-like components do not contribute to the M and Φ responses. After the multipole decomposition, the first line of Eq. (41) will contribute to F M , the second to F Φ , and the combination to an interference term between the two responses. Concentrating on the L = 0 multipole, the result takes a very compact form and is typically expressed as with the charge form factor defined by The proton/neutron combinations are related to the isospin components by and we have replaced the full form factors in Eq. (41) by the first terms in the expansion (7). The charge form factor fulfills the normalization F ch (0) = 1, and the corresponding representation (42) is exact for spin-0 nuclei. For nonvanishing spin, there are further form factors, e.g., the magnetic form factor for spin-1/2 in analogy to the nucleon, but for the reasons given in Sec. II E these contributions are small in heavy nuclei. In addition, twobody effects only enter at loop level, so that in contrast to the magnetic form factor two-body modifications of the charge form factor are also small. Finally, we give the corresponding expansion for the charge radius where R 2 p is the so-called point-proton radius defined as and r 2 so represents the spin-orbit contribution encoded in Φ [119]. In the case of Eq. (43) the matching of matrix elements and Wilson coefficients is trivial, since so far only the long-range contribution in the SM has been taken into account. A potential modification would be provided by the electron analog of L (5) given in Eq. (1). In the next step, we extend the discussion to short-range contributions from Z exchange, which are responsible for PVES in the SM.

B. Parity-violating electron scattering
The central observable in PVES is the asymmetry where the cross sections involve left-or right-handed initial-state electrons, respectively. The corresponding lepton traces are are the vector and axial-vector weak charges of the electron in the normalization of Ref. [54]. The terms in Eq. (48) involving an tensor will lead to SD corrections, which we will study below in the context of CEνNS, while the remainder follows in close analogy to the parityconserving case, the only difference being that the electromagnetic form factor needs to be replaced by its weak analog. With quark-level Wilson coefficients as in the SM and matrix elements from Eq. (10), the result takes the simple form where the weak charge has been separated from the weak form factor F w (q 2 ). However, we note that, in contrast to the electromagnetic charge and F ch (q 2 ), Q w does not actually factorize. The explicit definition reads structure factor 133 Cs where we have used that in the SM the Wilson coefficients for d-and s-quarks are identical to write the strangeness contribution in terms of Q n w . The corresponding weak radius reads with spin-orbit contribution [120] Finally, radiative corrections lead to the shifts [54,129,130] and accordingly for Q w .
As an example, Fig. 1 shows the M and Φ responses for 133 Cs. The coherent and partially coherent characters of M and Φ , respectively, are clearly observed at q = 0, where about 20%-25% of the nucleons contribute coherently for F Φ N (0). The minimum of F M n at lower |q| compared to F M p indicates a larger neutron than proton radius. Explicit parameterizations of all nuclear structure factors are provided in App. E.
We obtain the charge and weak radii given in Table III. In addition, Table III also shows the so-called neutron skin, defined as the difference between neutron and proton point radii, R n − R p . Calculated charge radii are in good agreement with experiment, similar to other approaches [17,25,27]. The disagreement between calculations increases for predictions of the weak radii and neutron skin. The shell model generally predicts larger weak radii and especially larger neutron skins than other many-body approaches [17, 23-25, 27, 127, 128].
The corresponding results for the weak form factors are shown in Figs. 2-5. In each case, we show the shell-model results for the modulus of the weak form factor including all corrections given in Eq. (52). Coherence is kept until larger momentum transfers in lighter nuclei with smaller neutron radius, see In the case of 40 Ar, Fig. 3 compares our results to the RMF calculation of Ref. [23] as well as ab initio results from coupled-cluster theory [27]. All calculated weak form factors give similar results, within the uncertainty band estimated in Ref. [27]. This suggests that uncertainties in the neutron distribution are relatively small, in contrast to the assumptions in Ref. [37]. We stress that apart from the nuclear structure, minor differences in the weak form factor arise from the precise input for the hadronic matrix elements and weak charges, primarily the proton charge radius, for which Refs. [23,27]

C. Neutrino scattering
The dominant contribution to the CEνNS cross section in the SM involves the same nuclear form factor as in the case of PVES, since apart from overall prefactors the combination of Wilson coefficients, hadronic matrix elements, and nuclear structure factors remains unchanged. This dominant piece of the differential cross section takes the form where E ν is the energy of the incoming neutrino and the nuclear recoil   is the leptonic trace whose components determine the spin sums The spherical components are defined with respect to the direction of q = k − k, e.g., In particular, the combination L 33 is strongly suppressed by T /m A 2E 2 ν /m 2 A , while L 03 or the additional terms in L 00 and L ii are only suppressed by T /E ν 2E ν /m A . In consequence, the longitudinal multipoles in Eq. (A1) can be safely neglected. The interference with the Coulomb multipoles could in principle become relevant, but the longitudinal multipoles involve an additional suppression by q 0 /|q| = − T /(2m A ) −E ν /m A from the application of current conservation, see Eq. (A3). Accordingly, all potentially relevant subleading kinematic effects can be taken into account by in Eq. (56). Next, there could be interference terms between the vector and axial-vector responses. The vector contributions to the transverse multipoles vanish for T → 0 and are not coherent, so the only potentially relevant interferences arise from the longitudinal and Coulomb multipoles. However, all such interferences vanish due to Eq. (A3).
Therefore, the dominant correction to Eq. (56) comes solely from the axial-vector part of the interaction. This contribution becomes relevant for precision studies of nuclei with nonvanishing spin, especially, because in contrast to other less relevant corrections their contribution remains finite in the limit T → 0. The SD structure factors are obtained by adapting the formalism from Ref. [29], most notably, by only keeping the transverse electric multipoles, due to the strong suppression of the longitudinal ones (transverse magnetic multipoles do not contribute to elastic scattering due to time reversal). Collecting the kinematic factors, the resulting contribution to the CEνNS cross section takes the form where the structure factors S T ij (q 2 ) are the same as for dark matter except that longitudinal multipoles need to be omitted, see Sec. III E as well as Apps. A and B for the precise definitions. In particular, the normalizations are related to S N , the nucleon (proton and neutron) spin expectation values: 2  Table IV, see also Refs. [28,29]. The isoscalar/isovector coefficients are , where g N A = q C A q g q,N A . In the SM we have, using Eqs. (4), (13), and (14), so that the full expression for the cross section becomes and the dominant SD correction arises, as expected, from the isovector component. The induced pseudoscalar form factor G P (t) only contributes to the longitudinal multipoles, see Eq. (A3). Since g A factorizes, the radius corrections from Eq. (17) are usually absorbed into the structure factors, as are corrections from two-body currents, to which we will turn in the next subsection.

D. Improved treatment of axial-vector two-body currents
Axial-vector currents are responsible for SD scattering. In the nonrelativistic limit the leading one-body (1b) currents are given by so that axial responses are driven by the nucleon spin S i = σ i /2, as indicated by Table II. A sizable correction to the leading one-body terms comes from subleading axial-vector two-body currents [32]. In medium-mass and heavy nuclei, these contributions have been evaluated in previous studies of β decays [132,133] and WIMP-nucleus scattering [28,29]. However, the studies of SD WIMP scattering off nuclei focus on pion-exchange two-body currents proportional to the low-energy couplings c 3 , c 4 , and c 6 [28,29] and neglected the contact two-body axial-vector current proportional to the couplings d 1 , d 2 [32], which is only included in the |q| = 0 limit in β decay [132,133].
Here we improve previous studies by including all pion-exchange, pion-pole, and contact terms derived in Ref. [32]: where k i = p i − p i , q = −k 1 − k 2 , and (1 ←→ 2) applies to the entire expression except for the last line. Relativistic 1/m N corrections to Eq. (68), besides the term proportional to p 1 + p 1 , can be absorbed into c 4 → c 4 + 1/(4m N ), c 6 → c 6 + 1/m N , where we use a dimensionful c 6 for consistency with the previous literature on the axial current (note that our choice of c 6 corresponds to c 6 /m N in the conventions of Ref. [134]). In the counting of Refs. [32,135] these relativistic corrections are formally of higher order, but we keep them both for consistency with Ref. [133] and in analogy to our treatment of higher-order effects in the c i , see below.
Following Refs. [28,29] we approximate the twonucleon currents by a normal-ordering approximation with respect to spin-isospin symmetric reference state with density ρ = 2k 3 F /(3π 2 ) (k F is the Fermi momentum) where P ij is the exchange operator and the sum is performed over the second nucleon j.
As a result, axial-vector two-body currents transform into effective one-body currents [29,136]: These two effective currents have the same structure as the two terms in the leading one-body current, Eq. (67), so they can be treated in the same way. The currents in Eqs. (70)-(71) depend on the nuclear density ρ, the momentum transfer q, and the combined momentum P. Because the dependence on P is small [29], in practice we evaluate the expressions taking P = 0. Likewise, we neglect additional effective one-body currents proportional to P and P · σ i . The functions I σ 1 (ρ, K), I σ 2 (ρ, K), I P (ρ, K), and I c6 (ρ, K) appear due to the summation over occupied states in the exchange terms in Eq. (69). They can be expressed as integrals, with analytical expressions given in Ref. [29].
In the P = 0 approximation, the combined effective currents can be written in analogy to Eq. (67) where For β decays q 0, and axial-vector two-body currents have been studied beyond the normal-ordering approximation in Eq. (69) [133]. The approximation for J eff i,2b was found to be very good when taking ρ ∼ 0.10 fm −3 , which is a typical value for the density of the nuclear surface. Based on this, for our evaluation of the nuclear structure factors we consider the density range ρ = 0.09 . . . 0.11 fm −3 . This range includes slightly lower densities, but is consistent with the one considered in Refs. [28,29].
The contributions from two-body currents in Eqs.    to be used should in principle be given by the nuclear interaction used to solve the many-body problem for the nucleus of interest. However, accurate many-body calculations using chiral interactions that depend explicitly on c i , c D are still not available for all nuclei discussed in this work. Instead, our results are based on many-body calculations that use shell-model Hamiltonians, which, despite being based on nucleon-nucleon interactions, are modified by phenomenological adjustments in order to improve their description of the nuclear structure of selected regions of nuclei. Therefore, we cannot use consistent c i , c D couplings between the nuclear interactions and the two-nucleon currents given in Eqs. (70)- (71). Our strategy is as follows. First, we use the values for c 1 , c 3 , and c 4 determined in the Roy-Steinerequation analysis of πN scattering [98,137]. This improved determination of the c i values allows us to obtain results with reduced theoretical uncertainties compared to Ref. [28,29], which considered a broad range of c 3 and c 4 (the smaller c 1 contributions are included for the first time in this work). In fact, at a given chiral order the uncertainties in the c i are now negligible, with the main uncertainty arising from the chiral expansion. Strictly speaking, one should use the next-to-leading-order values from Refs. [98,137] to be consistent with the chiral order we use for the axial-vector current, but this assumes that the latter is affected by large loop corrections in the same way as πN scattering, which is known not to be the case. Instead, we make use of the fact that the two-nucleon axial-vector current is matched to the three-nucleon force [135], in such a way that the leading loop corrections in the axial-vector current coincide with the ones in the three-nucleon force [138,139]. These corrections can be represented by a simple shift δc i [140] The values shown in Table V are then obtained as the combination of the next-to-next-to-leading-order values from Refs. [98,137] in combination with these δc i (as well as the relativistic correction for c 4 ), and the uncertainties represent the shifts between the two chiral orders. The value of c 6 is related to the isovector magnetic moments via [134] where we have indicated the leading loop correction. Similarly to the other c i , this correction is large despite being formally of higher order (in part due to the enhancement by a factor of π [141]). However, similar corrections arise from chiral loops in the axial-vector current [135,142], the dominant of which can again be represented as a shift in c 6 and cancels the matching correction in Eq. (76). Including the relativistic corrections discussed before, we will thus equate c 6 = (κ p − κ n + 1)/m N = 5.01, as given in Table V.
We then fix the value of the contact coupling c D , while at the same time correcting for the shortcomings of our phenomenological calculations. Shell-model nuclear matrix elements involving the axial-vector current typically overestimate experiment [143] by about 20% to 30%. Recently, Ref. [133] showed for β decay (where it is sufficient to take |q| = 0) that this is because of a combination of missing two-body axial-vector currents, see Eq. (68), and additional nuclear correlations that are beyond the standard shell-model approach. In order to account for this, we adjust the value of c D so that our shell-model calculations receive a contribution from two-nucleon currents such that, at |q| = 0, Eq. (70) reduces the leading term in Eq. (67) in the range 20% to 30%. The q dependence of the effective two-body currents is the one predicted by Eqs. (70)- (71). Since the leading contribution from twobody axial-vector currents comes from the pion-exchange part proportional to c 3 and c 4 , the part considered in Refs. [28,29], our results are consistent with these previous calculations.
The values of c i and c D used in this work are summarized in Table V, where the extreme values c D = −6.08 (c D = 0.30) only correspond to the low density ρ = 0.09 fm −3 (high density ρ = 0.11 fm −3 ). In practice, we neglect the remaining uncertainties in the c i due to effects from higher chiral orders not captured here, as those are subleading compared to the uncertainty in the range of c D values, which also depend on the nuclear density ρ. Ultimately, our uncertainty depends on the range imposed on the impact of the two-body currents at |q| = 0, 20%-30%, as estimated from β decay [133,143].

E. Spin-dependent responses for CEνNS and dark matter
The nuclear responses for CEνNS and SD dark matter scattering off nuclei can be expressed in terms of the transverse and longitudinal SD structure factors F Σ L ± (q 2 ) and F Σ L ± (q 2 ), respectively. For CEνNS, only the transverse component contributes, while for dark matter scattering both longitudinal and transverse parts need to be taken into account.
The expressions are given by which can be expressed in terms of the proton/neutron instead of the isoscalar/isovector basis as where the proton/neutron combinations are related to the isospin ones analogously to Eq. (44) The terms δ (q 2 ), δ (q 2 ) encode the corrections beyond the leading SD coupling to the transverse and longitudinal SD responses, respectively. They capture the combined effect of the pseudoscalar form factor, radius corrections, and two-body currents. They are given by where the two-body current contributions δa(q 2 ) and δa P (q 2 ) are defined in Eqs. (73)- (74).
Note that currents proportional to (q · σ i )q only contribute to the longitudinal multipoles. Moreover, their contribution can be treated similarly to terms proportional to σ i because where the second term is perpendicular to q and vanishes for longitudinal multipoles.
As a first application we show the results for the structure factors S N (q 2 ) for xenon, in comparison to our previous work from Ref. [29], see Fig. 6. There is good consistency within the earlier theoretical band. As expected, recent progress in the understanding of low-energy constants and two-body currents in β decays allows us to reduce the theoretical uncertainties. Figure 6 shows that for xenon this is especially the case for S p , as this response is dominated by two-body contributions. In general, uncertainty bands are reduced most for the smaller structure factors corresponding to the species with even number of nucleons.
Second, we show the variant of the SD structure factors required for CEνNS, see Figs. 7 and 8. As discussed in Sec. III C, only the transverse multipoles contribute to the final expression in Eq. (66), but unless the strangeness contribution is neglected all isospin components enter. The figures show our shell-model results, including two-body currents and form factor corrections represented by δ (q 2 ), δ (q 2 ) in Eq. (84). For a given nucleus, the shape of the isovector and isoscalar responses is similar because all of them are ultimately dominated by either S p (q 2 ), if the nucleus has an unpaired proton, or S n (q 2 ), for nuclei with odd number of neutrons. A comparison between the 131 Xe structure factors in Figs. 6 and 8 shows that the shape of the transverse component may differ significantly from the total structure factor (dominated by the longitudinal component in that case, see Ref. [29]). According to Eq. (63), the normalization of the transverse contribution differs by 2/3 from the sum. Moreover, as can be seen from Figs. 7 and 8, the isovector combination S T 11 , which is most relevant for Eq. (66), is the smallest of the isospin components. This is partly because of the reduction caused by axial-vector two-body currents, which are isovector, as one-body S 11 and S 00 structure factors are of similar size.
, for xenon. The dark bands refer to the results from this work, the light bands to the ones from Ref. [29].

A. Vector and axial-vector operators
As a first step, we generalize Eq. (66) to include scenarios in which still only vector and axial-vector operators are present, but whose Wilson coefficients are allowed to deviate from the SM. Especially the case with BSM contributions only to the vector operators is a frequently studied scenario [2,3].
To collect the combination of Wilson coefficients and hadronic matrix elements, we define as well as the short-hand notation wherė   Fig. 7, for the two odd-mass xenon isotopes. and the neutron equations follow by g p For the strangeness contribution we have introduced the "baryon-number" coupling In the SM, where C V d = C V s , this new coupling coincides with g n V and was therefore not needed in Eq. (52). Collecting all terms, the generalization of Eq. (66) becomes where the isoscalar and isovector couplings for the axialvector part are defined as in Eq. (64). In the vector contribution, the new "weak charge," reduces to −G F / √ 2 Q w in the SM, see Eq. (51), and the new "weak form factor" becomes Modifications due to BSM physics thus affect the CEνNS cross section in two ways: the normalization at q 2 = 0  changes, visible as the change in the weak charge, but in addition the weak form factor changes as well, which is due to the fact that Q w does not actually factorize, but emerges as a sum of different underlying nuclear responses. Only in special cases in which the shifts in the Wilson coefficients are aligned with the SM, i.e., all coefficients are modified by the same relative factor, would F w (q 2 ) remain unaltered.
To quantify the changes with respect to F w (q 2 ), the new form factor is shown in Fig. 9 for several points in the BSM parameter space. These contributions to the u-and d-quark vector Wilson coefficients, defined as in Eq. (5), are large but realistic in view of current bounds from CEνNS [2,3]. By definition, the deviations vanish at |q| = 0, and they become most visible in the vicinity of the zeroes. The second point is illustrated in Fig. 10, which shows that sufficiently far away from the zeros the changes are at the few-percent level, while the relative deviations are enhanced once the process becomes less coherent.

B. Operators not present in the Standard Model
Next, we turn to the operators in Eq. (1) not present in the SM. At dimension 5 there is only the dipole operator, leading to the lepton trace where we dropped terms that vanish upon contraction with the nuclear matrix element due to gauge invariance.
Since the interference terms with the SM contribution FIG. 10: Relative changes in the weak form factor for 133 Cs, for the same scenarios shown in Fig. 9.
vanish, the presence of a dipole contribution would manifest itself as a new, long-range interaction One power of 1/t from the photon propagator in the squared matrix element cancels with the lepton trace in Eq. (93), but the second remains and leads to the divergence for T → 0, due to the relation between momentum transfer and nuclear recoil given in Eq. (57). Next, the lepton trace for the scalar operator is The diagonal term in the cross section can be expressed as This expression vanishes for T → 0, but otherwise there is no kinematic suppression compared to the vector contribution due to the scaling m A T /(2E 2 ν ) 1. We have collected all the relevant couplings and form factors in the scalar combination F S , which is defined as with F M N given in Eq. (44), the two-body contributions F π (q 2 ), F b (q 2 ) from Ref. [36], and the following combi-nations of Wilson coefficients and hadronic couplings Again, there is no interference with the SM, but the scalar contribution does interfere with the dipole, leading to For the pseudoscalar operator there is also no interference with the SM, and due to the SD nature of the nucleon matrix elements such a response should be even further suppressed than in the scalar case. To corroborate that expectation we rewrite the operator by means of the axial Ward identitȳ so that we can define a leptonic trace to be contracted with the same nuclear responses already studied for the axial-vector case. The relevant spin sums are given by L 33 = L ii = t 2 /4, leading to a kinematic suppression with respect to the axial-vector contribution that scales as The scale m N emerges assuming that the formal difference between the dimension-7 and dimension-6 operators is mainly due to hadronic scales (as is manifest for the matrix elements of the scalar operator, see Eq. (22)), and for higher scales the suppression would be even stronger.
In either case we conclude that pseudoscalar contributions to CEνNS are negligible. For the tensor operator, the most relevant contributions are expected from the space-like components σ ij , because only those are momentum independent and not suppressed by 1/m N in the nonrelativistic expansion. For the same reason, the induced terms in Eq. (21) are subleading. The result of the multipole decomposition for tensor currents, see App. D, then leads to the following expressions: defining the couplings via and the cross section becomes Contrary to the axial-vector response, there is now also a contribution from the longitudinal multipoles,S L ij (q 2 ). These response functions are identical to the ones derived for the axial-vector case only at leading order, i.e., the two-body corrections for the tensor current would take a different form and likewise the corrections from the induced pseudoscalar and the axial-vector radius need to be removed .
(106) There are again no interference terms with the SM, but the lepton traces do allow for potential interference terms with scalar, pseudoscalar, and dipole operators. In addition, there would be additional contributions from the σ 0i components of the tensor current as well as the induced form factors in Eq. (21). In case such contributions became relevant, the formalism could be extended accordingly.

V. SUMMARY
In this paper we have provided a detailed account of the CEνNS cross section both within the SM and beyond. To this end, we started from a decomposition into effective operators, hadronic matrix elements, and nuclear structure factors, including both the vector and axialvector operators already present in the SM, but also considering the effects of (pseudo-)scalar, tensor, and dipole operators. Light BSM degrees of freedom could be included along similar lines.
As first step, we introduced the charge and weak form factors as typically defined in electron scattering, to exemplify their decomposition in terms of underlying nuclear structure factors, but also hadronic matrix elements and Wilson coefficients. The analogous decomposition for CEνNS is then used to address the question how, e.g., the weak form factor needs to be modified once BSM contributions are permitted, and to derive master formulae for the cross section in the various cases.
Our results for the nuclear structure factors are based on the large-scale nuclear shell model. In addition to the coherent part of the response, which is largely determined by charge operators, radius and relativistic corrections, as well as spin-orbit contributions, we have also performed a detailed study of the typically neglected axialvector responses. While the general formalism is similar to the spin-dependent responses for dark matter scattering off nuclei, there are key differences. Most notably, only the transverse multipoles contribute to CEνNS due to the lepton trace. We have also calculated updates for the structure factors relevant for spin-dependent dark matter scattering. 3 Our calculation of the spin-dependent responses takes advantage of several developments in recent years that allow us to improve the treatment of two-body currents as predicted from chiral EFT. These include improved determinations of the relevant low-energy constants from pion-nucleon scattering, the calculation of one-loop corrections to the nuclear axial-vector current, and insights from ab initio studies of two-body effects in mediummass and heavy nuclei. While the nuclear interactions used in this work are still phenomenological, this strategy allows us to incorporate as many constraints from chiral EFT as possible, including, for the first time, the effect of contact operators and pion-pole contributions to the two-body currents.
Finally, we provide further details of the multipole expansion of the nuclear responses, tailored towards the aspects relevant for the CEνNS application and making the connection to the notation in the nuclear-physics literature. Together with the fits of the resulting nuclear responses as well as the EFT decomposition of the cross section, this defines general CEνNS responses for a wide range of isotopes and effective operators. SFB 1245, and the Max Planck Society.

Appendix A: Multipole expansion
In this appendix we review the main features of the multipole expansion, following closely Refs. [113,118]. The starting point is the leptonic current l µ , which is decomposed into the temporal component l 0 and the spatial, spherical components l λ , λ = ±, 3 with respect to the reference vector q, where the latter index is chosen to avoid confusion with the temporal component. The spin sum takes the form where the reduced matrix elements refer to the longitudinal (L), Coulomb (M), transverse electric (T el ), and transverse magnetic (T mag ) multipoles. The latter can be simplified to The single-nucleon contributions, obtained by nonrelativistic expansion of Eqs. (9) and (11), can then be expressed in terms of fundamental multipole operators according to where we dropped the quark labels for the form factors, terms suppressed by q 0 /m N , and several subleading multipoles in the axial-vector contribution. The explicit expressions for the multipoles in harmonic oscillator basis are given in Ref. [113], where an additional operator Ω L = ∆ L − Φ L is introduced. Not all multipoles will be needed in the analysis, the most important ones are M and Φ for the vector responses and Σ , Σ for the SD ones. The nuclear responses Σ L , ∆ L , as well as the combinations ∆ L − 1 2 M L , Ω L + 1 2 Σ L vanish for elastic scattering. Excitation energy (keV) 133  state with angular momentum and parity J P = 5/2 + , the difference with the 7/2 + state is only 10 keV. The angular momentum and parity of the lowest energy levels is predicted well, even though the energy of the calculated second 3/2 + state is lower than in experiment. Overall the agreement with experiment is similar as in other odd-mass nuclei with similar mass number.

Appendix D: Multipole decomposition for tensor currents
Including the tensor operator from Eq. (1) into the analysis requires a generalization of the multipole decomposition reviewed in Sec. A. Here we follow closely the original derivation in Refs. [144][145][146], including the lepton trace L µνλσ = Tr / k σ µν P L / kσ λσ P R = 2 g µλ g νσ − g µσ g νλ k · k + i µνλα k σ k α − i µνσα k λ k α − i µλσα k ν k α + i νλσα k µ k α − g µλ k ν k σ + k ν k σ + g µσ k ν k λ + k ν k λ + g νλ k µ k σ + k µ k σ − g νσ k µ k λ + k µ k λ , and then specify the spin sums relevant for CEνNS. The key idea in the generalized multipole expansion is then that the antisymmetric tensor current j µν essentially admits two vectorial components, j where we dropped the distinction between the two parities in each multipole. Since the nonrelativistic reduction of σ 0i only starts at O(1/m N ) and depends on momenta, the most interesting tensor contribution originates from the σ ij → ijk σ k components, contained in j (1) . The relevant spin sum reads with projections spins l (1) In contrast to Eq. (59), the longitudinal multipole is no longer kinematically suppressed, but instead the interference term between electric and magnetic multipoles can be dropped. In our normalization the hadronic current starts with − i √ 2 ijk σ jk → −i √ 2σ i , so that, up to the prefactor and the different lepton traces, the remainder of the calculation follows along the same lines as for the axial-vector response.  i=0 c Φ ± i u i , with u = q 2 b 2 /2. These forms correspond to the analytical solution in the harmonic-oscillator basis [115,147], with nM and n Φ as implied by the table. Our results for xenon are given in Ref. [33], the ones for germanium in Table VII