Constraining New Physics with Borexino Phase-II spectral data

We present a detailed analysis of the spectral data of Borexino Phase II, with the aim of exploiting its full potential to constrain scenarios beyond the Standard Model. In particular, we quantify the constraints imposed on neutrino magnetic moments, neutrino non-standard interactions, and several simplified models with light scalar, pseudoscalar or vector mediators. Our analysis shows perfect agreement with those performed by the collaboration on neutrino magnetic moments and neutrino non-standard interactions in the same restricted cases and expands beyond those, stressing the interplay between flavour oscillations and flavour non-diagonal interaction effects for the correct evaluation of the event rates. For simplified models with light mediators we show the power of the spectral data to obtain robust limits beyond those previously estimated in the literature.


Introduction
The scattering among electrons and neutrinos played an important historical role in the establishment of the Standard Model (SM) allowing for precise tests of the V-A structure of the weak interactions. Being a purely weak, purely leptonic, two-body reaction, makes both the theoretical cross-section calculation and the experimental signatures particularly clean. The scattering of ν µ,τ on electrons is a purely weak neutral-current (NC) process, while ν e e − scattering proceeds via both weak NC and weak charged-current (CC) interactions. The latter is of particular interest because it is one of the few reactions for which the SM predicts a large destructive interference between these two channels. For the same reasons it is a sensitive probe of physics beyond the Standard Model (BSM). Paradigmatic examples include neutrino electromagnetic properties, non-standard neutrino-matter interactions (NSI), and new interactions mediated by light bosons. Most precise results on ν µ e − → ν µ e − were obtained by the CHARM-II experiment [1], while ν e e − → ν e e − was first observed with an accelerator beam in the experiment performed at LAMPF [2], followed by the measurement at the LSND experiment [3]. Since the pioneering work of Reines and collaborators [4] the scattering processν e e − →ν e e − has been studied at a greater level of precision using neutrinos produced in nuclear reactors [5][6][7], most recently in the TEX-ONO [8] and GEMMA experiments [9,10]. Neutrino electron elastic scattering (ES) can also be searched for at experiments designed to detect Coherent Elastic neutrino-Nucleus Scattering (CEνNS): while the signal in the SM would be too low to yield a significant number of events, the presence of BSM physics may lead to much larger contributions. In fact, CEνNS searches using reactor neutrinos at the CONUS [11,12], CONNIE [13], and Dresden-II reactor experiments [14,15] lead to very strong constraints on some of the BSM scenarios outlined above. Similar searches can be performed using pion decay-at-rest (DAR) sources: for example, the data of the COHERENT experiment on CsI [16] should be sensitive to an ES signal as well [15] since it would pass their signal selection cuts.
Neutrino electron elastic scattering is also key to the study of solar neutrinos. Solar neutrinos were first detected in real time through ES scattering in water in Kamiokande [17]. In fact, the observed rate of ES versus CC interactions allowed for the first model-independent evidence of flavour transitions in solar neutrinos in Super-Kamiokande and SNO [18,19]. For the last 15 years the Borexino [20] experiment has been detecting ES of solar neutrinos in liquid scintillator performing precision spectroscopy of their energy spectrum at Earth [21][22][23]. The main purpose of the experiment was the independent determination of the fluxes produced by the different thermonuclear reactions in the Sun. Nevertheless, the precision of their high-statistics spectral data also allows for precision tests of ES and, as such, it has been used by the collaboration itself [24,25]. Borexino results have also been analyzed/adapted by phenomenological groups to impose bounds on BSM scenarios (see, e.g., Refs. [26][27][28][29][30][31][32][33][34][35] for an incomplete list), in general using either total event rates or early spectral data.
In this work we present a detailed analysis of the spectral data of Borexino Phase II with the aim of exploiting its full potential to constrain new physics which can significantly affect ES (such as neutrino magnetic moments, NSI, and simplified models with light scalar, pseudoscalar, and vector mediators). In general BSM scenarios, neutrino neutral-current interactions may be flavour changing. This brings the issue of the interplay between neutrino oscillations and interactions which do not conserve the neutrino flavour, which has not always been properly treated in previous literature on this topic. Section 2 presents the formalism to correctly evaluate the event rates, accounting for flavour oscillations in presence of flavour non-diagonal interactions (details are also presented in Appendix A). This section also introduces the interaction cross sections for the BSM models considered. Section 3 contains the description of our calculation of the expected spectra at Borexino and the dataset used, as well as our treatment of backgrounds, uncertainties and further details on the statistical analysis. It is complemented with Appendix B.
The main results of our work are presented in Sec. 4. As mentioned above, the Borexino collaboration has already performed two analyses on BSM searches using a spectral fit of the Phase-II data set: a first analysis which is used to set constraints on the neutrino magnetic moment [24], and a later search for NSI with electrons [25]. As we will show in Sec. 4, our results are in excellent agreement with those obtained in Ref. [24] for the magnetic moment scenario: this serves as validation for our χ 2 construction and, in particular, for the implementation of systematic uncertainties (a key point in this analysis, given the large number of events collected by Borexino). Regarding NSI, our analysis expands over the results obtained in Ref. [24] in two main ways: (i) we include all NSI operators in neutrino flavour space at once in the fit, thus allowing for potential correlation effects and cancellations in the cross section; and (ii) we consider four different cases regarding the structure of the operators (vector, axial vector, left-handed only, and right-handed only). Finally, our work includes the derivation of new constraints for simplified models with new light vector, scalar, and pseudoscalar mediators, which to the best of our knowledge have not been reported by the collaboration yet. Finally, we summarize our conclusions in Sec. 5.

Theoretical Frameworks
In the SM, the differential cross section for neutrino-electron ES (ν β e − → ν β e − ) can be expressed as [36] dσ SM where T e stands for electron (kinetic) recoil energy, E ν is the neutrino energy, and y ≡ T e /E ν , and f + , f − , f ± are loop functions given in Ref. [36], while α stands for the finestructure constant, G F is the Fermi constant, and m e is the electron mass. The effective couplings c Lβ and c Rβ contain the contributions from the SM NC, which are equal for all flavours, plus the CC contribution which only affects ν e e − scattering: with ρ and κ β (T e ) departing from 1 due to radiative corrections of the gauge boson selfenergies and vertices [36]. The minimum neutrino energy required to produce an electron with a given recoil energy T e is E ν,min = T e + T 2 e + 2 m e T e 2, while the maximum kinetic energy for an electron produced by a neutrino with given energy E ν is T e,max = E 2 ν (E ν + m e /2). The corresponding cross section forν β e − ES can be obtained from (2.1) with the exchange c Lβ ↔ c Rβ .
We assume that scattering off electrons is incoherent, so for the calculation of the number of events it is enough to multiply the total cross section σ SM β (E ν ) -obtained by integrating Eq. (2.1) over T e -by the overall number of electrons in the fiducial volume. In principle, a correction factor taking into account that the target electrons are actually bound to atoms have to be considered for recoil energies comparable to atomic binding energies. However, for the solar neutrino energies considered here this effect can be safely neglected.
Solar neutrinos produced in the core of the Sun are purely ν e . Propagation through vacuum and matter induces a non-trivial evolution of their flavour, so that neutrinos reaching the detector are in general not definite flavour eigenstates but rather a coherent superposition of them, which can be described in terms of the elements S βe of an unitary matrix S (see Appendix A). As discussed above, in the SM the ES cross section includes contributions from both CC interactions (which are diagonal in the flavour basis) and NC interactions (which are diagonal in any basis), so that the flavour basis is the most adequate to describe them. The common approach consists in projecting the neutrino state at the detector over the flavour eigenstates, thus introducing the ν e → ν β transition probabilities P eβ ≡ |S βe | 2 , and then convoluting them with the cross sections for pure flavour states given in Eq. (2.1). So the number of interactions of solar neutrinos in Borexino will be proportional to the probability-weighted electron scattering cross section: This well-known expression relies on the special role played in the SM by the flavour basis, which as we have seen diagonalizes both CC and NC interactions. However, through this work we will consider several BSM phenomenological scenarios leading to modifications of the interaction cross section in Eq. (2.1). When the BSM is not lepton-flavour conserving, the CC-defined flavour basis may no longer correspond to the NC eigenstates, and extra care must be taken with the interplay of oscillations and scattering. To address this issue, it is convenient to rewrite Eq. (2.3) in a basis-independent form: where ρ (e) ≡ S Π (e) S † is the density matrix of the neutrino state at the detector, with Π (e) being the projector on the electron neutrino flavour at the source. Here σ SM is a generalized cross section, i.e., a matrix in flavour space which contains enough information to describe the ES interaction of any neutrino state without the need to explicitly project it onto the interaction eigenstates. Such projection is now implicitly encoded into σ SM , which in the flavour basis takes the diagonal form σ SM γβ = δ γβ σ SM β . In the same basis we have ρ (e) βγ = S βe S * γe , so that ρ (e) ββ = P eβ and Eq. (2.3) is immediately recovered. Let us stress that the indexes γβ in σ SM γβ do not make reference to "initial" and "final" states of the scattering process but they are both initial flavour indexes, just as the indexes of ρ (e) βγ both refer to the final state of the neutrino as it reaches the detector.

Neutrino Non-Standard Interactions with Electrons
The so-called non-standard neutrino interaction (NSI) framework consists on the addition of four-fermion effective operators to the SM Lagrangian at low energies. For example, the effective Lagrangian would lead to new NC interactions with the rest of the SM fermions. Here P can be either a left-handed or a right-handed projection operator (P L or P R , respectively) and the matrices ε f,P αβ are hermitian: ε f,P βα * = ε f,P αβ . The index f refers to SM fermions and in what follows we will be focusing on f = e. Such new interactions may induce lepton flavour-changing processes (if α = β), or may lead to a modified interaction rate with respect to the SM result (if α = β). For later convenience we also define the vector and axial SM couplings: as well as the corresponding NSI coefficients: The presence of flavour-changing effects implies that the SM flavour basis no longer coincides with the interaction eigenstates. In such case the number of events in Borexino will be given by Eq. (2.4) in terms of a generalized cross section σ NSI (E ν ) defined as the integral over T e of the following matrix expression: 1 which is formally identical to Eq. (2.1) except that the c Lβ and c Rβ real coefficients have been superseded by the 3 × 3 hermitian matrices C L and C R : It is immediate to see that, if the NSI terms ε e,L αβ and ε e,R αβ are set to zero, the matrix dσ NSI dT e becomes diagonal and the SM expressions are recovered. The corresponding cross section forν β e − ES can be obtained from Eq. (2.8) with the exchange C L ↔ C * R . In addition to scattering effects, the presence of NSI also affects neutrino propagation from production to detection point through its contribution to the matter potential both in the Sun and in the Earth. The formalism follows the approach described in Sec. 2.3 of Ref. [37]. The details of the calculation of the density matrix ρ (e) are provided in Appendix A.

Neutrino Magnetic-Moments
In the presence of a neutrino magnetic moment (µ ν ) the scattering cross sections on electrons get additional contributions which do not interfere with the SM ones. Neutrino magnetic moments arise in a variety of BSM models and, in particular, they do not need to be 1 Eq. (2.8) is already summed over the outgoing neutrino state, which is assumed to be undetected.
The expression for a specific final state ν f can be obtained by inserting the corresponding projector Π (f ) between each pair of CL,R matrices, i.e., by replacing C 2 L → CLΠ (f ) CL, C 2 R → CRΠ (f ) CR and CL, CR → CLΠ (f ) CR + CRΠ (f ) CL. Since f Π (f ) = I, Eq. (2.8) is recovered upon summation over f . flavour-universal. Therefore, in what follows we will allow different magnetic moments for the different neutrino flavours -albeit still imposing that they are flavour-diagonal. Under this assumption the ES differential cross section for either neutrinos or antineutrinos, up to order O(y 2 ), takes the form [38] dσ µν Since no flavour-changing contributions are present in this case, the expression for the number of interactions in Borexino will be just proportional to the probability weighted total cross section, in analogy with Eq. (2.3).

Models with light scalar, pseudoscalar, and vector mediators
We also consider simplified models describing the neutrino interaction with fermions of the first generation through a neutral scalar φ, pseudoscalar ϕ, or a neutral vector Z . Concerning the scalar, the Lagrangian we consider is [39] where q j φ are the individual scalar charges and j = {ν α , e}. The corresponding cross section for neutrino-electron (or antineutrino-electron) scattering is where we have included all orders in T e /E ν and T e /m e (beyond the leading-order terms given in Ref. [39]), since these are non-negligible for solar event rates at Borexino. For the pseudoscalar scenario, we consider the Lagrangian where q j ϕ are the individual pseudoscalar charges and j = {ν α , e}. The corresponding cross section for neutrino-electron (or antineutrino-electron) scattering is (2.14) As for the vector mediator, our assumed Lagrangian is: with the individual vector charges q j Z . Unlike for the scalar and magnetic-moment cases, a neutral vector interaction interferes with the SM contribution. The differential cross section for neutrino-electron scattering reads [40] Model q e q νe q νµ q ντ Universal/leptonic scalar (or pseudoscalar) g 2 Z (q tions in Borexino will be just proportional to the probability-weighted total cross section as in Eq. (2.3).
In Sec. 4 we study the constraints on four specific models with light mediators. For scalar and pseudoscalar mediators we will focus on a model in which the mediator couples universally to fermions. For vector mediators we will show the results for three anomalyfree 2 models coupling to electrons, namely B − L, L e − L µ and L e − L τ . For convenience Tab. 1 lists the relevant charges.

Analysis of Borexino Phase-II Spectrum
In order to constrain the BSM scenarios for neutrino interactions described in Sec. 2 we have performed an analysis of the Borexino Phase-II data collected between December 14th, 2011 and May 21st, 2016, (corresponding to a total exposure of 1291.51 days × 71.3 tons) following closely the details presented by the collaboration in Ref. [22]. In particular we use the spectrum data presented in Fig. 7 of such work which we analyze using the total number of detected hits (N h ) on the detector photomultipliers (including multiple hits on the same photomultiplier) as estimator of the recoil energy of the electron scattered by the incoming neutrinos. We statistically confront the spectral data with the expected event rates in the different scenarios, which we compute as follows.
Borexino phase-II detects interactions of solar neutrinos produced in the pp, pep, 7 Be, and 8 B reactions of the pp-chain, as well as from all reactions in the CNO-cycle. The produced total flux is described by the differential spectrum: where E ν is the neutrino energy and f refers to each of the main components of the solar neutrino flux 3 (pp, pep, CNO, 7 Be, and 8 B). In their journey to the Earth the produced ν e will change flavour and interact upon arrival via elastic scattering with the e − of the detector. Thus, the expected number of solar neutrino events corresponding to N h hits is given by the convolution of the oscillated solar neutrino spectrum with the interaction cross section and the energy resolution function. Quantitatively we compute the solar neutrino signal S i in the i-th bin (i.e., with N h ∈ [N i h,min , N i h,max ]) as: Here, dS dT e is the differential distribution of neutrino-induced events as a function of recoil energy of the scattered electrons (T e ) where X ∈ {SM, NSI, Z , φ, ϕ, µ ν } for the different scenarios discussed in Sec. 2. Note that our notation explicitly reflects that dσ dT e may not be diagonal in flavour space (as in the case of NSI, see the discussion in Sec. 2). In our analysis we fix the oscillation parameters to the best-fit values from Ref. [41,42]. Here N tgt is the number of e − targets, that is, the total number of electrons inside the fiducial volume of the detector (corresponding to 71.3 ton of scintillator), while T run = 1291.51 days is the data taking time, and E cut = 98.5% is the overall efficiency. In addition Eq. (3.2) includes the energy resolution function dR dT e for the detector which gives the probability that an event with electron recoil energy T e yields an observed number of hits N h . We assume it follows a Gaussian distribution

4)
whereN h is the expected value of N h for a given true recoil energy T e . We determinē N h (T e ) and σ h (T e ) via the calibration procedure described in Appendix B.1. The data contains background contributions from a number of sources. The main backgrounds come from radioactive isotopes in the scintillator 14 C, 11 C, 10 C, 210 Po, 210 Bi, 85 Kr, and 6 He. The collaboration identifies two additional backgrounds due to pile-up of uncorrelated events, and residual external backgrounds (see Ref. [22] for a complete description of these backgrounds).
Among the considered backgrounds, the one coming from 11 C is particularly relevant for the analysis. This isotope is produced in the detector by muons through spallation on 12 C. In Ref. [22] the collaboration uses a Three-Fold Coincidence (TFC) method to tag the 11 C events, which are correlated in space and time with a muon and a neutron. Thus, they divide the Phase II data set in two samples: one enriched (tagged) and one depleted (subtracted) in 11 C events. The separation of the solar neutrino signal (as well as most of the other backgrounds) into these two samples is uncorrelated from the number of hits N h . Concretely, the tagged sample picks up 35.72% of the solar neutrino events, while the subtracted sample accounts for the remaining 64.28%.
In our analysis we include both sets of data, denoted in what follows by s = "tagged" or "subtracted". As for the background, we have read the contribution B c s,i for each component c, in each bin i, and for each data set s from the corresponding lines in the two panels in Fig. 7 of Ref. [22]. They are shown for the best-fit normalization of the different background components as obtained by the collaboration, and we take them as our nominal background predictions in the absence of systematics. Altogether the nominal number of expected events T 0 s,i in some bin i of data sample s is the sum of the neutrino-induced signal and the background contributions: where the index f ∈ {pp, 7 Be, pep, CNO, 8 B} runs over the solar fluxes, while the index c ∈ { 14 C, 11 C, 10 C, 210 Po, 210 Bi, 85 Kr, 6 He, pile-up, ext} runs over the background components. In our calculation of S f s,i we use the high-metallicity (HZ) solar model for simplicity. 4 In Fig. 1 we show the predicted event distributions for some examples of the new physics models we consider. Along with these spectra, we also include the Borexino data points, after subtracting the best-fit background model and SM signal. For the sake of concreteness we show the results for the "subtracted" event sample. As seen in the figure, the BSM scenarios considered here yield larger contribution in the low energy bins. Therefore the derivation of robust bounds requires a detailed treatment of the background subtraction and careful accounting of all systematic uncertainties, in particular those associated with the energy reconstruction. This was also the case for the analysis performed by the Borexino collaboration in Ref. [22] for the determination of the solar fluxes, in particular the pp flux. We follow closely such analysis in our construction of the χ 2 function as we outline next. We provide a detailed description of our treatment of backgrounds and systematic uncertainties in Appendix B.
Our statistical analysis is based on the construction of a χ 2 function based on the described experimental data, neutrino signal expectations and sources of backgrounds, as well as the effect of systematic uncertainties. We introduce the latter in terms of pulls ξ r , which are then varied in order to minimize the χ 2 . The analysis includes a total of 27 pulls to account for the uncertainties in the background normalizations, detector performance, priors on the solar fluxes, and energy shifts (see Appendix B for details). Assuming a linear dependence of the event rates on the pulls, the minimization can be performed analytically provided the χ 2 can be assumed to be Gaussian. Borexino collaboration bins their data in each sample using 858 bins, from N h = 92 to N h = 950. Such thin binning does not provide enough events in each bin to guarantee Gaussian statistics. So for convenience we rebin the data using coarser bins: the first 42 bins are rebinned in groups of 2, the next 210 bins in groups of 5, the next 260 bins in groups of 10 and the last 346 bins in groups of 50. This yields a total of 96 bins k for each data sample, s, each of them having enough events to guarantee Gaussian statistics. Altogether the χ 2 we employ reads: where O s,k is the observed number of events in bin k of sample s and T s,k is the corresponding predicted number of events in the presence of systematics D r s,k stands for the derivatives of the total event rates with respect to each pull, which are given in Appendix. B.
The last two terms in Eq. (3.6) are Gaussian bias factors which we include for those pulls for which some prior constraint should be accounted for. The first of those last terms correspond to pulls associated to uncorrelated ("unc") uncertainties while the last one to correlated ("corr") uncertainties. Here, Σ stands for the covariance matrix used for the latter and it is given in Appendix B.

Results
This section summarizes our results from the spectral analysis of Borexino data, as outlined in Sec. 3, for the phenomenological scenarios defined in Sec. 2. We will first present our analysis for the magnetic moment scenario introduced in Sec. 4.1. As we will see, our results are in excellent agreement with those obtained in Ref. [24] and therefore serve to validate our χ 2 implementation (and, in particular, the implementation of systematic uncertainties) discussed in Sec. 3. We then proceed to the NSI scenario in Sec. 4.2, while the cases of light vector, scalar, and pseudoscalar mediators are shown in Secs. 4.3 and 4.4.
At this point it is worth stressing the technical challenges associated to the χ 2 minimization procedure. The analysis, a priori, depends on 27 pulls, 6 standard oscillation parameters, plus a number of non-standard oscillation parameters (ranging from 2 to 6, depending on the particular scenario considered). In order to make the problem tractable from the computational point of view, we minimize analytically over the 27 pulls as outlined in the previous section, and we fix the oscillation parameters to the present best-fit values from Ref. [41,42] (we have checked that this has a negligible impact on the results). The exploration of the BSM parameter space is performed numerically, using the Diver software [43,44]. We have also cross-checked that the resulting two-dimensional confidence regions agree with the result obtained using a simplex method algorithm.
Throughout this section, we will follow a frequentist approach in order to determine the allowed confidence regions for the parameters in each case. In other words, for a given BSM scenario characterized by a set of parameters ω, we define a ∆χ 2 function as: where χ 2 min refers to the global minimum of the χ 2 function, found after minimization over all the model parameters ω. We stress that χ 2 ( ω) corresponds to Eq. (3.6) where T s,k ( ξ) is in fact T s,k ( ω, ξ) so the minimization over the pull parameters is performed for each given value of the model parameters. 5 The allowed regions are then determined by the set of model parameters for which the ∆χ 2 lies below a certain level, as indicated in each case.
For reference let us mention that the analysis performed within the SM results into χ 2 min,SM = 256.7 for a total of 192 data points and 27 pull parameters. Results for the magnetic moment scenario, assuming three independent parameters (one for each of the SM neutrino flavours). Since the neutrino magnetic moment is positive definite, the cuts have been taken according to a one-sided χ 2 distribution. In the central panels, we show the allowed confidence regions (at 1σ and 2σ, for 2 d.o.f. using one-sided intervals) while in the right-most panels we show the profile of the ∆χ 2 for each of the parameters individually. In each panel, the results are obtained after minimization over the undisplayed parameters.

Neutrino magnetic moment
For the analysis with neutrino magnetic moment we find χ 2 min,SM+µν = 256.5. In other words, the data do not show a statistically significant preference for a non-vanishing magnetic moment, which allows to set constraints on the parameter space for this scenario. The corresponding results are summarized in Fig. 2, assuming three independent parameters for the ν e , ν µ and ν τ magnetic moments. In each panel, the ∆χ 2 is minimized over the undisplayed parameters. From the one-dimensional projections we can read the upper bounds at 90% CL (for 1 d.o.f., one-sided intervals, i.e., ∆χ 2 = 1.64) Our results are fully compatible with those in Ref. [24], considering the different choice of oscillation parameters assumed. If we assume that all flavours carry the same value of the neutrino magnetic moment then the dependence on the oscillation parameters cancels out and we obtain the following bound (at 90 % CL for 1 d.o.f., one-sided) This coincides precisely with the result obtained by the collaboration in Ref. [24], and serves as validation for the χ 2 implementation and the choice of systematic uncertainties adopted in our analysis. We notice that the Borexino bounds are stronger than those obtained by other experiments under the same flavour assumptions and CL. This is the case for the (over a decade old) experimental results of GEMMA [9,45] (µ ν < 2.9 × 10 −11 µ B ) and TEX-ONO [8] (µ ν < 7.4 × 10 −11 µ B ), as well as the bounds derived from recent results such as CONUS [46] (µ ν < 7.5 × 10 −11 µ B ) and the combined analysis of Dresden-II and COHER-ENT [15] (µ ν < 1.8 × 10 −10 µ B ).

Non-Standard neutrino Interactions with electrons
Let us now discuss the case of NSI with electrons. In this scenario Borexino is sensitive to a total of 12 operators: 6 operators involving left-handed electrons (with coefficients ε L αβ ) and 6 operators involving right-handed electrons (with coefficients ε R αβ ), see Eq. (2.5). In all generality, a priori all operators should be considered at once in the analysis since all of them enter at the same order in the effective theory. However, including such a large number of parameters makes the problem numerically challenging, and the extraction of meaningful conclusions is also jeopardized by the multiple interference effects between the different coefficients. Thus, here we focus on four representative cases with a reduced set of operators, depending on the type of interaction that could lead to the terms in Eq. We start discussing the axial NSI case. For axial NSI, the effect cancels out in the matter potential, rendering neutrino oscillations insensitive to the new interactions. However, Borexino can still constrain these operators since they affect the interaction cross section in the detector. First, we notice that even allowing all 6 axial NSI coefficients to vary freely we get χ 2 min,SM+NSI A = 255.7, which is only one unit lower that the SM. So we proceed to derive the bounds on the full parameter space.
The results are presented in Fig. 3 where we show one-dimensional and two-dimensional projections of ∆χ 2 SM+NSI A after minimization over the parameters not shown in each panel. We notice that since the cross-section of flavour-conserving interactions receive contributions from both SM and BSM terms, which could lead to cancellations among them, the corresponding regions are generically larger than those involving only flavour-off-diagonal NSI, for which no SM term is present. From the figure we see that when all the parameters are included in the fit, Borexino is able to constrain NSI roughly only at the level  of |ε A αα | < O(2) for flavour-diagonal interactions, and |ε A α =β | < O(1) for the flavour-offdiagonal ones (see also Tab. 2). The only exception to this is the parameter ε A ee , for which tighter constraints are obtained, roughly at the level |ε A ee | < O(0.5) as it can be seen from the lowest row in the figure. The main reason behind this is that the ε A ee parameter enters the effective cross section in Eq. (2.5) multiplied by P ee ∼ 0.5, while the remaining NSI parameters are multiplied either by P eµ , P eτ ∼ 0.2 − 0.25 or by off-diagonal elements in the density matrix (which are at most ∼ 0.2).
Note also that the allowed regions for flavour-diagonal NSI are a bit asymmetric with respect to zero. In general, given the large parameter space it is difficult to identify a single reason for the specific shape of the allowed regions. Still, the origin of the asymmetry can be qualitatively understood if the generalized cross section matrix in Eq. (2.8) is written in terms of the axial and vector coefficients, defined (in analogy to Eq. (2.9)) as (4.4) In this case, the differential cross section for axial NSI, setting α = β, depends on the following combination of parameters: We notice that due to the interference between the SM and NSI parameters, a second SMlike solution appears in Eq. (4.5) with ε A αα −2 c Aα (i.e., for which C A αα = −c A,α ). It can actually be shown that for α = µ and α = τ such SM-like solution holds also after including the full y dependence of the cross section and so it provides a fit of similar quality to the SM solution. When only one NSI parameter is considered at a time this leads to the two disconnected allowed ranges seen for ε A µµ and ε A τ τ in Table 2. When all NSI parameters are included both allowed regions merge in a unique range which is asymmetric about zero. In the case of ε A ee , on the contrary, the y-dependent terms break such degeneracy and the second SM-like solution gets lifted.
From Eq. (4.5) we also see that for ε A αα −c Aα it is possible to reduce the cross section with respect to the SM case. This leads to a slightly better fit to the data (albeit mildly, at the level of ∆χ 2 = −1 as mentioned above) and results into the preferred solution to be at ε A ee 0 (since c Ae ∼ 0.5). Conversely, for ε A µµ and ε A τ τ the fit shows a mild preference for positive values since c Aµ = c Aτ ∼ −0.5.
Finally, it is worth stressing the impact of correlations in the fit, which are large for this scenario. For example, from the two-dimensional panels in Fig. 3 it is evident that there is a strong correlation between ε A µτ , ε A µµ , and ε A τ τ , which has a significant impact on the final limits obtained when all parameters are allowed to float freely in the fit. This is shown explicitly in Tab. 2, which compares the obtained bounds (at 90% CL) for the different NSI parameters in two different cases: turning on only one parameter at a time ("1 Parameter"), or including all operators simultaneously in the fit ("Marginalized"). As can be seen from the comparison between the two sets of results in this table, the impact due to correlations and degeneracies between different operators is significant and leads to considerably larger allowed regions if multiple NSI operators are included at once. For example, in the case of ε A µµ , the 90% CL range is reduced by a factor of ∼ 8 if the remaining NSI operators are set to zero.
Let us now turn to the analysis of vector NSI. In this case NSI affect the matter potential felt by neutrinos in propagation and a relevant question is how sensitive the analysis is to this effect, and whether the sensitivity to NSI is still dominated by the impact on the neutrino detection cross section. We have numerically checked that the results of the analysis are rather insensitive to NSI effects in propagation. This is partly so because in the energy window considered here the contribution from the higher-energy components (mainly 8 B, but also pep neutrinos) to the event rates is subdominant, and moreover it is affected by relatively large backgrounds uncertainties. As a result we find that the impact Allowed regions at 90% CL (∆χ 2 = 2.71)

Vector
Axial Vector of the matter potential on the fit is very limited and the sensitivity comes dominantly from NSI effects on the detection cross section. 6 We also find that, even allowing all 6 vector NSI coefficients to vary freely in the analysis, χ 2 min,SM+NSI V − χ 2 min,SM is only −0.2. Our allowed regions for vector NSI are presented in Fig. 4, where we see that the potential of Borexino to test NSI with electrons is again limited in this multi-parameter scenario: Borexino is only able to constrain NSI at the level of |ε V αα | < O(2 − 3) for flavourdiagonal interactions, and |ε V αβ | < O(1) for the off-diagonal parameters. Precise values for the allowed regions at 90% CL are given in Tab. 2.
The features describing the allowed regions in the two-dimensional projections of the ∆χ 2 are qualitatively similar to the results found for the axial case, albeit with some differences. In particular the regions for ε V µµ and ε V τ τ are slightly more symmetric than the axial case. Again, this can be qualitatively understood if the cross section in Eq. (2.8) is rewritten in terms of the axial and vector couplings which results in Eq. (4.5) but replacing c Aα ↔ c V α and ε A αα ↔ ε V αα . Thus, a similar interference as in the axial case takes place here; however since c V µ = c V τ ∼ −0.04 the quasi-degenerate SM-like solution lies closer to the SM so the region is more symmetric about zero.
For completeness, we also provide our results for NSI involving only left-handed and right-handed electrons in Tab. 3, where we again show the results obtained turning on only one operator at a time, and allowing all NSI operators simultaneously in the fit. We have numerically checked that, if only a single operator is included, we recover the results of the Borexino collaboration for both left-handed and right-handed NSI (Fig. 5 in Ref. [25]) with a very good accuracy. The only difference is the secondary minimum for ε L ee ∼ −1.3 which is outside of the range of parameters explored in Ref. [25].   when only one NSI operator is included at a time ("1 Parameter") or when the remaining NSI coefficients are allowed to float freely in the fit ("Marginalized"). Our results are compared to our derived bounds from GEMMA [45,50], TEXONO [50,51], COHERENT [15,16,52], and Dresden-II reactor experiment [14,15], as well as those derived from a global fit to oscillation data [50], all of them obtained with the same statistical criterion. We also show the bounds from CONNIE [13] (95% CL) and CONUS [54] (90% CL) which we have obtained applying an overall rescaling factor to their published bounds on a universally coupled vector (see text for details). For the B − L model, the line labelled "Borexino Rate Only" shows the approximate bound derived (at 95% CL, 1 d.o.f., for one-sided intervals) in Ref. [26] (see also Ref. [27]).

Light vector mediators
The potential of Borexino to probe simplified models with light vector mediators has been demonstrated in the literature for the B − L case [26,27], and the bound has been later recast to other U (1) models in Refs. [47,48]. However, Refs. [26,27] date from 2012 and 2015 and therefore did not include Phase-II data. Furthermore, their bound was only approximate: it was derived requiring that the new physics contribution to the event rates should not exceed the SM expectation. Here we report our limits for simplified models, for our data analysis which includes full spectral information for the Phase II data and a careful implementation of systematic uncertainties. Our results for vector mediators are presented in Fig. 5 for three anomaly-free models in which the vector couples to B − L, L e − L µ , or L e − L τ . As seen in the figure the results for the three models are very similar. This can be easily understood from Eqs. (2.3) and (2.16). In particular the dominant BSM contribution to the event rates arises from the interference term with the SM. So, the new physics will lead to a modification to the expected number of events that goes as where we have removed the sum over neutrino flavours, as c V e ∼ 1 c V µ(τ ) ∼ 0.04 and P ee > P eµ P eτ . Since for the three models considered q νe Z q e Z = 1, the derived bounds are about the same in the three cases.
For completeness, we also compare our Borexino results to previous bounds obtained in the literature. For the range of couplings considered here, the most relevant constraints come from CEνNS, neutrino scattering on electrons, and from neutrino oscillation data, while measurements of neutrino scattering on nuclei or electrons at high momentum transfer will be suppressed with the inverse of the momentum transfer squared (see, e.g., the discussion in Ref. [49]). In principle, additional bounds are obtained from fixed-target experiments or colliders (for example, the Z may be produced on-shell and decay to pairs of charged leptons). However, these bounds are a priori model-dependent: for example, if the new mediator decays predominantly to dark sector particles, some of them can be significantly relaxed. Thus, here we consider only those bounds from neutrino scattering and oscillation experiments, which always apply to the models considered (for a detailed compilation of constraints on light vector mediators we refer the interested reader to, e.g., Refs. [48,50]). Specifically, we consider the following set of constraints: • measurements of ν − e scattering at the reactor experiments GEMMA [45] and TEX-ONO [51], computed as in Ref. [50]. Here we note that our analysis of GEMMA Phase-II shows a poor fit to the data (with or without BSM contributions). Naively including this dataset in the ∆χ 2 analysis results into strong bounds on the BSM models which, however, are not statistically robust. Thus, the results shown here for GEMMA correspond to Phase-I data alone, 7 for which (under the SM hypothesis) we numerically find χ 2 /n d.o.f. ∼ 26/39, n d.o.f. being the number of bins.
• bounds from CEνNS measurements (on both CsI and Ar nuclei) at COHERENT [16,52] and at the Dresden-II reactor experiment [14], computed as in Ref. [15] (based on the analysis performed in Refs. [50,53]). Note that, although these experiments were designed to search for CEνNS (which dominates the sensitivity to B − L), they are also sensitive to ν-e elastic scattering and therefore can be used to set constraints on the L e − L µ and L e − L τ models as well (see Ref. [15] for details). For concreteness, here we show the results for Dresden-II using the quenching factor labelled as "Fef" in Ref. [15] (qualitatively similar results are obtained for the "YBe" quenching factor).
• constraints derived from the impact of new vector mediators on the matter potential in neutrino oscillations. These were obtained from a global fit to oscillation data in Ref. [50]. Note that these constraints only apply to flavoured U (1) models (where the neutrino flavours have different charges under the new interaction) and therefore are absent in the B − L case.
• bounds derived from searches for CEνNS at CONNIE [13], and from CEνNS and ES searches at CONUS [54]. In both cases the bounds were reported for a vector mediator which couples universally to all SM fermions. The bounds obtained from ES searches coincide for the B − L and the universal models, since q νe Z q e Z = 1 in the two cases. Regarding the bounds from CEνNS searches, the reported limits for the universal model are rescaled for the B − L case, assuming an average momentum transfer. 8 For CONNIE we quote the results obtained using the Lindhard quenching factor [55]. For CONUS, on the other hand, we take the results corresponding to a quenching factor 9 for k = 0.16 from Fig. 5 in Ref. [54]. Notice that the excluded regions derived by CONUS in Ref. [54] are given as 90% CL (see Ref. [12] for details), while for CONNIE [13] the excluded regions are shown at 95% CL.
• for the B − L model, we also show previous estimates of the Borexino constraint computed in Ref. [26,27]. These results are labelled as "Rate only". As outlined above, these are obtained using data available prior to 2015, requiring that the new physics contribution to the total event rates in the energy region corresponding to the 7 Be line does not exceed the SM expectation (within reported uncertainties, which at the time were at the level of 8%). These limits were reported as one-sided constraints, at 95% CL with seemingly 1 d.o.f.
As can be seen, our results for Borexino qualitatively agree with the approximate bounds from Ref. [26], but lead to an improved constraint by a factor ∼ 60% on the upper bound on Z for very light mediators (and by about 30% on g Z /M Z in the contactinteraction regime). Also, as seen in the figure, our analysis of Borexino Phase-II spectra leads to improvements on the bounds from scattering experiments in most of the parameter space, for the three models under consideration. On the other hand, the bounds from oscillation data prevail as the strongest constraints on the L e − L µ and L e − L τ models.

Light scalar and pseudoscalar mediators
Next let us discuss our results for a scalar (pseudoscalar) mediator coupled universally to the SM fermions, shown in the left (right) panel of Fig. 6. It is important to note that constraints from oscillation data do not apply to scalar (or pseudoscalar) mediators. For comparison, we also show the bounds derived from neutrino scattering measurements, computed following (or obtained from) the same references as in Sec. 4.3. Comparing the scattering constraints shown in Fig. 6 against the same bounds for vector mediators in Fig. 5, we find that, generically, the bounds imposed on the scalars (and more so for pseudoscalars) are somewhat weaker. This is expected because the dominant contribution for the vector interactions arise from the interference with the SM (which is quadratic on the new coupling constant) while the scalar (or pseudoscalar) contributions do not interfere with the SM and therefore the dependence on the new coupling constant is quartic. Furthermore 8 This approximation allows to factor out of the cross section the dependence on the BSM model parameters, through a proportionality constant that only depends on the weak charges of the nucleus in the SM and the BSM cases. We have checked that this approximation is excellent as we are able to correctly reproduce their allowed regions for the universal vector scenario. 9 As we will see this choice has no impact on the allowed regions for vector mediators since in this case the bound is driven by the ES contribution; however, in the scalar case some sensitivity arises from CEνNS searches as well and therefore there is a slight dependence on this choice.
as seen in Eqs. (2.12) and (2.14) the cross section in the scalar and pseudoscalar scenarios are suppressed as T e m e /E 2 ν and T 2 e /E 2 ν , respectively, with respect to the vector-mediated scattering cross section in Eq. (2.16).
From the figure we see that our bounds from Borexino Phase II data improve over previous constraints from ES reactor experiments (GEMMA and TEXONO) in both panels. However, for scalar interactions (left panel) the strongest limits come from the Dresden-II reactor experiment, and Borexino is only able to give a comparable limit for very light mediator masses M φ ∼ 0.5 MeV. This is so because for the universal model shown the scalar coupling to quarks would lead to additional contributions to CEνNS at the Dresden-II experiment (as well as for CONNIE and CONUS), which enhances their sensitivity in the high-mass region (see, e.g., Ref. [15]). Conversely, in the pseudoscalar scenario our bounds from Borexino dominate for mediator masses M ϕ 5 MeV, with TEXONO being the only other experiment which provides a comparable result. This is so because pseudoscalars interactions couple to spin and therefore the corresponding CEνNS cross section can be safely neglected. Thus, the Dresden-II contour shown in the right panel has been obtained considering only the new physics contribution to ES (which would pass their signal selection cuts), leading to a worse constraint compared to the one obtained in the scalar case. Finally, note that while in principle it should also be possible to derive bounds on the pseudoscalar scenario using ES searches at CONUS, the collaboration did not study this scenario and therefore no limits from CONUS are available in this case. From the results for the scalar case we expect, however, that these would be weaker than the bound obtained from Borexino, since the cross section in the pseudoscalar case would be suppressed to that in the scalar scenario by an additional factor of T e /m e (at CONUS, this factor can be estimated as O(keV/0.5 MeV) ∼ 2 × 10 −3 ).

Conclusions
In this work we have presented a detailed analysis of the spectral data of Borexino Phase II (corresponding to a total exposure of 1291.51 days × 71.3 tons) with the aim at exploiting its full potential to constrain BSM scenarios. Following closely the analysis of the collaboration, we have performed a detailed treatment of the background subtraction and careful accounting of all systematic uncertainties (particularly those associated with the energy reconstruction), including a total of 27 nuisance parameters in our fit. We have validated our constructed χ 2 by performing a fit to the solar neutrino fluxes which accurately reproduces the results of Borexino [22], as shown in Figs. 8 and 9.
We have used our analysis to derive constraints on neutrino magnetic moments, NSI, and several simplified models with light scalar, pseudoscalar, and light vector mediators. In particular, for NSI we consider at the same time flavour off-diagonal as well as flavourdiagonal interactions. This brings the issue of the interplay between neutrino oscillations and interactions which do not conserve neutrino flavour, which has been often overlooked in the literature. We correctly account for those in terms of a density matrix (see Eq. (2.4)), which we provide for the specific case of NSI in Appendix A. 10 −2 10 −1 10 0 10 1 Figure 6. Bounds from our analysis of Borexino Phase-II spectral data on a scalar (left panel) and pseudoscalar (right panel) mediator (at 90% CL for 1 d.o.f., using two-sided intervals, i.e., ∆χ 2 = 2.71) which couples universally to all fermions in the SM. Our results are compared to our derived bounds from GEMMA [45,50], TEXONO [50,51], COHERENT [15,16,52], and Dresden-II reactor experiment [14,15], obtained using the same statistical criterion. For comparison we also show the bounds from CONNIE [13] (95 %CL) and CONUS [54] (90% CL).
Our main results are presented in Sec. 4. For neutrino magnetic moments we have allowed different magnetic moments for the different neutrino flavours -albeit still imposing that they are flavour-diagonal -and obtained the bounds in Eq. (4.2). For flavourindependent magnetic moment we obtained the bound in Eq. (4.3). This is in perfect quantitative agreement with that obtained by the collaboration in Ref. [24] and serves as an additional validation of our χ 2 implementation.
Our analysis for NSI includes all operators with the same Lorentz structure: purely vector, axial vector, left-handed only, or right-handed only (6 parameters in each case, assumed to be real). Our results are shown in Figs. 3 and 4 as well as in Tables 2 and 3. We have numerically checked that, if only one NSI operator is included at a time, we recover the results in Ref. [24] for the left-handed only and right-handed only scenarios, to a very good approximation. However, we find that when all couplings are allowed simultaneously, strong correlations appear and the bounds are substantially relaxed, as explicitly shown in Tables 2 and 3.
We have also derived new constraints for simplified models with new light vector, scalar, and pseudoscalar mediators which we show in Figs. 5 and 6. For mediator masses M 0.1 MeV, the obtained experimental bound becomes independent of the mass of the mediator. Within this (effectively massless) regime, the bounds derived from Borexino data for the considered models read, at 90% CL (1 d.o.f., one-sided) Conversely, for sufficiently large mediator masses, the contact-interaction approximation is recovered and the event rates vary with a power of (|g|/M ) n with n = 4 for scalar mediators and with n ranging between 2 and 4 for vector mediators. In this limit we get the following bounds at 90% CL (1 d.o.f., one-sided) It is convenient to write U = OU 12 , where O = R 23 R 13 is real while U 12 is a complex rotation by angle θ 12 and phase δ CP . Then we can write: In order to further simplify the analysis, let us now assume that all the mass-squared differences involving the "heavy" states ν 3 can be considered as infinite: ∆m 2 3l → ∞ for l = 1, 2. In leading order, the matrixH takes the effective block-diagonal form: where H (2) is the 2 × 2 sub-matrix ofH corresponding to the first and second neutrino states. Consequently, the evolution matrix is: We are interested only in the elements S αe . From the definition of O we see that O e2 = 0.
Taking into account the block-diagonal form ofS, we obtain: In order to derive explicit expression, let us begin writing H (2) = H (2) vac + H mat with: 4E ν − cos 2θ 12 sin 2θ 12 e iδ CP sin 2θ 12 e −iδ CP cos 2θ 12 , The coefficients ε f D and ε f N are related to the original parameters ε f αβ by the relations: Finally, it is convenient to introduce the following effective quantities:

A.1 Neutrino conversion probabilities
Let us begin from the flavour conversion probabilities, P αe = |S αe | 2 . From Eq. (A.5) it is straightforward to derive: 11 S (2) * 21 where we have used the fact that the term containing a factor e −i∆ 33 L oscillates very fast, and therefore vanish once the finite energy resolution of the detector is taken into account. This expression can be generically written as: (A.14) From Eq. (A.11) we get:  16) in full analogy with Eq. (A.12). As before, this expression can generically be written as: General expressions for the A αβ , B αβ , C αβ and D αβ matrices can be easily obtained.

B.1 Energy Calibration
As mentioned in the text, we perform the fit to Borexino data in terms of the observed number of hits N h as estimator of the recoil energy of the electron scattered by the incoming neutrinos. Energy resolution is included as the Gaussian function in Eq. (3.4) which gives the probability that an event with electron recoil energy T e yields an observed number of hits N h in terms of two parameters, the mean value of N h for a given T e ,N h (T e ), and the width of the probability distribution σ h (T e ). To determine these two parameters as functions of T e we proceed in in two steps.
First we obtain the dependence of σ h N h from the calibration results of Borexino published in Ref. [56]. Concretely, by fitting the data from the upper panel of Fig. 22 of the published version (Fig. 21 in the arXiv version) with Gaussian distributions, we derive a relation betweenN h and σ h , which is depicted in the left panel of Fig. 7 and which reads: The calibration results were obtained with a γ-ray source. And the basic assumption is that Eq. (B.1) holds equally for electrons. Second we introduce that resolution function, which is now a function ofN h , in our evaluation of the expected event rates for the different solar fluxes and infer the dependence ofN h (T e ) which best reproduces the solar spectra depicted by Borexino in Fig. 7 of Ref. [22]. We assume a 3rd-degree polynomial relation and find the best agreement to be with The agreement we find between the Borexino solar spectra and our own calculations can be seen in the right panel of Fig. 7.

B.2 Implementation of systematics uncertainties, and assumed priors
In our analysis we introduce a total of 22 uncorrelated pulls to include background normalization uncertainties, signal and background systematic uncertainties related to the detector performance, and energy shift uncertainties. In addition 5 correlated pulls are used to introduce the SSM priors on the solar flux components. Here we discuss the implementation these different types of uncertainties separately.

Background normalizations
We consider 9 different background sources. Seven correspond to radioactive isotopes contaminating the scintillator, 14 C, 11 C, 10 C, 210 Po, 210 Bi, 85 Kr, and 6 He. Two additional additional backgrounds are due to pile-up of uncorrelated events, and residual external backgrounds (we refer in what follows to this last one as "ext"). We read the 18 (9 for each data sample) spectra and best fit normalization for these backgrounds from the results in Fig. 7 of Ref. [22]. Following the procedure described in that reference we allow their normalization to vary in the fit. We do so by introducing a set of 11 independent pulls of which 7 corresponding to the normalizations of 14 C, 10 C, 210 Bi, 85 Kr, 6 He, pile-up and ext which are assumed to be the same for both data samples, and 4 to the normalization of 11 C, 210 Po which are left to be different in the tagged and the subtracted samples. The normalizations corresponding to 10 C, 210 Bi, 6 He, 11 C tag , 210 Po tag , and 210 Po sub are left unconstrained and therefore for them no bias factor is included in the χ 2 . For the other normalizations we include priors following the information provided in Ref. [22].
The effect of these 11 pulls is included as described in Eq. where for convenience we label the r = 1, . . . , 11 pulls by a label c which specifies the background component they normalize.

Uncorrelated detector related signal and background uncertainties
We introduce eight pulls to parametrize the systematic uncertainties related to detector performance and background modelling. Their corresponding derivatives are: where σ r f and σ r c denote the assumed uncertainties for their effect each signal (S f s,i ) and background (B c s,i ) component, respectively. The values used in our fit are provided in the upper and lower sections of Tab. 4, respectively.  Table 4. Prior uncertainties assumed for each signal component σ r f (upper part of the table) and for each background component σ r c,sub = σ r c,tag (lower part), for the eight nuisance parameters used to account for systematic uncertainties in the fit. Here, ξ 1 − ξ 6 correspond to the six first sources of systematic uncertainties listed in Tab. IV of Ref. [22] where they only show their effect on the determination of the solar fluxes in the Borexino fit. ξ 2 corresponds to the choice of energy estimator which we find to be also to be relevant for the determination of the steeply falling 14 C background. ξ 3 corresponds to the pile-up modelling so we assign similar uncertainties for the effect of that pull on the pile-up background. ξ 8 corresponds to the combined uncertainty on the live time, scintillator density and fiducial volume uncertainty from the same table which we assume affect all signals and backgrounds. ξ 7 refers to an overall normalization error of the elastic scattering cross-section [36] (which only affects the signal).

Energy shift uncertainties
In addition to the systematic uncertainties describe above we also introduce a set of pulls to account for possible systematic shifts in the energy spectra for the different backgrounds that we have read from Fig. 7 of Ref. [22]. We came out to the need of this additional pulls by noticing that the Borexino data shown in that figure can be extracted in two different ways: directly from the upper panels in that figure, or using the residuals provided in the lower panels, and assuming that these have been obtained with respect to the best-fit lines shown in the upper panels. The results obtained in these two ways should be equivalent but we found that this was not the case.
From further inspection of the data obtained in these two ways, it was clear that they were partially shifted in energy with respect to one another. This is most evident for the signal and background components with steepest energy dependence which happen to occur in the lowest part of the spectrum where the statistics is higher, and therefore the effect in the quality of the fit most relevant. We find that the introduction of an add-hoc energy shift of about 0.25N h is able to reduce significantly the discrepancy observed between the two data sets, reconciling the two results. Thus, in our fit we use the data extracted from the upper panels in their Fig. 7, but we include additional pull terms in order to correct for a possible shift in energy. Since such energy shift has the largest impact in the low-energy part of the spectrum (and, most notably, on the pp solar neutrinos), we introduce three pulls: one for the 14 C and pile-up backgrounds, another one for the 11

Correlated uncertainties on the solar fluxes
The goal of our analysis of Borexino data is constraining the BSM scenarios presented in Sec. 2. For that we assume the solar fluxes to be as given by the Standard Solar Model (SSM). Technically this is enforced by introducing five pulls for the normalization of the five flux components of the solar neutrino signal subject to the priors of the SSM in both their central values, uncertainties and correlations. We introduce those by means of a covariance matrix Σ which we obtain from Ref. [57,58]. In our calculations, we use the high-metallicity (HZ) solar model for simplicity. 10 So in this case the corresponding derivatives are simply: As a first validation of our constructed χ 2 function we perform an analysis which intends to reproduce the fit of the collaboration in Ref. [22] to determine the solar neutrino fluxes. So we make a fit in which we do not include the covariance matrix for the solar fluxes previously described. Instead we leave the solar pp, 7 Be and pep neutrino signal fluxes completely free. We do introduce two priors for the CNO and 8 B the SSM (for concreteness we chose the high metallicity version of the model) as Borexino does. In the analysis we fit all background normalizations with the priors described above, and include the eight pulls for the uncorrelated systematics and the the three pulls for the energy shift. The result for our best fitted spectra (this is, our version of Fig. 7 in Ref. [22]) is shown in Fig. 8. The dependence of our ∆χ 2 on the normalization of the solar fluxes and the dominant backgrounds is shown in Fig. 9. The flux normalizations in the figure are normalized to the best fit values obtained by Borexino collaboration, hence a value of "1" means perfect agreement between our best fit fluxes and those of the collaboration. For comparison we also show for the solar fluxes as grey lines the inferred ∆χ 2 from the results quoted by the collaboration in Table II of Ref. [22] (assuming it is approximately Gaussian). The figures show that our constructed event rates, the best-fit normalization of all signal and background components, and the precision in the determination of the solar fluxes reproduce with very good accuracy those of the fit performed by the collaboration. Ref. [22] does not give enough information for a quantitative comparison of the ∆χ 2 for background normalizations, but we find that the precision we obtain for those is in good qualitative agreement with the results in Table III and Fig. 5 of Ref. [22].  Figure 9. Dependence of ∆χ 2 of fit to the Borexino phase-II spectra on the normalization of the solar fluxes and the dominant backgrounds (normalized to the corresponding best fit normalizations of the fit of the Borexino collaboration in Ref. [22]). For comparison we show as light blue curves the corresponding results of the determination of solar fluxes in Ref. [22]) (see text for details).