Final results of Borexino on CNO solar neutrinos

We report the first measurement of CNO solar neutrinos by Borexino that uses the Correlated Integrated Directionality (CID) method, exploiting the sub-dominant Cherenkov light in the liquid scintillator detector. The directional information of the solar origin of the neutrinos is preserved by the fast Cherenkov photons from the neutrino scattered electrons, and is used to discriminate between signal and background. The directional information is independent from the spectral information on which the previous CNO solar neutrino measurements by Borexino were based. While the CNO spectral analysis could only be applied on the Phase-III dataset, the directional analysis can use the complete Borexino data taking period from 2007 to 2021. The absence of CNO neutrinos has been rejected with>5{\sigma} credible level using the Bayesian statistics. The directional CNO measurement is obtained without an external constraint on the $^{210}$Bi contamination of the liquid scintillator, which was applied in the spectral analysis approach. The final and the most precise CNO measurement of Borexino is then obtained by combining the new CID-based CNO result with an improved spectral fit of the Phase-III dataset. Including the statistical and the systematic errors, the extracted CNO interaction rate is $R(\mathrm{CNO})=6.7^{+1.2}_{-0.8} \, \mathrm{cpd/100 \, tonnes}$. Taking into account the neutrino flavor conversion, the resulting CNO neutrino flux at Earth is $\Phi_\mathrm{CNO}=6.7 ^{+1.2}_{-0.8} \times 10^8 \, \mathrm{cm^{-2} s^{-1}}$, in agreement with the high metallicity Standard Solar Models. The results described in this work reinforce the role of the event directional information in large-scale liquid scintillator detectors and open up new avenues for the next-generation liquid scintillator or hybrid neutrino experiments.


Introduction
Solar neutrinos are produced in the core of the Sun by nuclear reactions in which hydrogen is transformed into helium.The dominant sequence of reactions is the so-called pp chain [1,2] which is responsible for most of the solar luminosity, while approximately 1 % of the solar energy is produced by the so-called Carbon-Nitrogen-Oxygen (CNO) cycle.Even though the CNO cycle plays only a marginal role in the solar fusion mechanisms, it is expected to take over the luminosity budget for main sequence stars more massive, older, and hotter than the Sun [3].Solar neutrinos have proven to be a powerful tool to study the solar core [4,5,6,7] and, at the same time, have been of paramount importance in shedding light on the neutrino oscillation phenomenon [8,9,10,11,12,13].
One important open question concerning solar physics regards the metallicity of the Sun, that is, the abundance of elements with Z > 2. In fact, different analyses of spectroscopic data yield significantly different metallicity results, that can be grouped in two classes: the so-called High-Metallicity (HZ) [14,15] and Low-Metallicity (LZ) [16,17,18] models.The solar neutrino fluxes, in particular that from the CNO cycle reactions, can address this issue.Indeed, the SSM predictions of the CNO neutrino flux depend on the solar metallicity directly, via the abundances of C and N in the solar core, and indirectly, via its effect on the solar opacity and temperature profile.
Borexino delivered the first direct experimental proof of the existence of the CNO cycle in the Sun with a significance of ∼7 σ, also providing a slight preference towards High-Metallicity models [4,19].This result was obtained with a multivariate analysis of the energy and radial distributions of selected events.To disentangle the CNO signal from the background, the multivariate fit requires an independent external constraint on the pep neutrino rate and on the 210 Bi rate; the latter is obtained by tagging 210 Bi- 210 Po coincidences in a temperature stabilized, layered scintillator fluid (see [4,19] for more details).For this reason, the CNO measurement has been performed only on approximately one third of the Borexino data, the so-called Phase-III.
In this paper, we present new results on CNO neutrinos obtained exploiting the "Correlated and Integrated Directionality" (CID) technique, which uses the directional information encoded in the Cherenkov light emitted alongside the scintillation, to separate the solar signal from non-solar backgrounds.Borexino demonstrated the viability of this technique using 7 Be solar neutrinos [20,21].Here we apply the CID technique to the CNO analysis, obtaining two important results: we show that we can extract the evidence of solar CNO neutrinos on the entire Borexino dataset following an alternative approach with respect to the standard multivariate analysis and, consequently, without the help of the 210 Bi constraint; we also show that by combining the information coming from the directionality with the standard multivariate analysis performed on Phase-III data we obtain an improved measurement of the CNO neutrino interaction rate.
The paper is structured as follows.Section 2 describes the Borexino detector and summarizes the event reconstruction techniques.
The CID analysis for the CNO neutrino measurement is illustrated in Sec. 3, outlining the methods, reporting the results, and detailing the main sources of systematic uncertainties.Finally, in Sec. 4 we show our best result on CNO neutrinos obtained combining the CID and the standard multivariate analysis.

The Borexino experiment
Borexino was a liquid scintillator (LS) neutrino detector [22] that ran until October 2021 with unprecedented radiopurity levels [23,5], a necessary feature of its solar neutrino measurements.The detector was located deep underground at the Laboratori Nazionali del Gran Sasso (LNGS) in Italy, with about 3800 m water equivalent rock shielding suppressing the cosmic muon flux by a factor of ∼10 6 .
The detector layout is schematically shown in Fig. 1.The Stainless Steel Sphere (SSS) with a 6.85 m radius supported 2212 8-inch photomultiplier tubes (PMTs) and contained 280 tonnes of pseudocumene (1,2,4-trimethylbenzene, PC) doped with 1.5% of PPO (2,5-diphenyloxazole) wavelenght shifter, confined in a nylon inner vessel of 4.25 m radius.The density of the scintillator was (0.878 ± 0.004) g cm −3 with the electron density of (3.307 ± 0.015) × 10 31 e − /100 tonnes.The PC-based buffer liquid in the region between the SSS and IV shielded the LS from external γ radiation and neutrons.The nylon Outer Vessel, that separated the buffer in two sub-volumes, prevented the inward diffusion of 222 Rn.The SSS itself is submerged in a domed, cylindrical tank filled with ∼1 kton of ultra-pure water, equipped with 208 PMTs.The water tank provided shielding against external backgrounds and also served as an active Cherenkov veto for residual cosmic muons passing through the detector.
Borexino detected solar neutrinos via their elastic scattering on electrons of the LS, a process sensitive, with different probability, to all neutrino flavors.Electrons, and charged particles in general, deposit their energy in the LS, excite its molecules and the resulting scintillation light is The effective total light yield was ∼500 photoelectrons per MeV of electron equivalent deposited energy, normalized to 2000 PMTs [23].The energy scale is intrinsically non-linear due to ionization quenching and the emission of Cherenkov radiation.The Geant4 based Monte Carlo (MC) software [24] simulates all relevant physics processes.It is tuned using the data obtained during calibration campaigns with radioactive sources [25].Distinct energy estimators have been defined, based on different ways of counting the number of detected photons [23,20].The position reconstruction of each event is performed by using the time-of-flight corrected detection time of photons on hit PMT [23].Particle identification is also possible in Borexino [23], in particular α/β discrimination [4], by exploiting different scintillation light emission time profiles.
The Borexino data-taking period is divided into three phases: Phase-I (May 2007-May 2010), Phase-II (December 2011-May 2016), and Phase-III (July 2016-October 2021).Phase-II started after the detector calibration [25] and an additional purification of the LS, that enabled a comprehensive measurement of the pp chain solar neutrinos [5].Phase-III is characterized by a thermally stable detector with greatly suppressed seasonal convective currents.This condition has made it possible to extract an upper limit constraint on the 210 Bi contamination in the LS, and thus, to provide the first direct observation of solar CNO neutrinos [4].

Correlated and Integrated Directionality for CNO
Cherenkov photons emitted by the electrons scattered in neutrino interactions retain information about the original direction of the incident neutrino.
Typically, in water Cherenkov neutrino detectors, this information is accessed through an event-by-event direction reconstruction, as demonstrated by the measurements of 8 B neutrinos, at energies larger than 3.5 MeV [26,7].Instead, the Borexino experiment has provided a proof-of-principle for the use of this Cherenkov hit information in a LS detector and at neutrino energies below 1 MeV through the so-called "Correlated and Integrated Directionality" (CID) technique.A detailed explanation of the method can be found in [21,20].
The CID method discriminates the signal originating in the Sun -due to solar neutrinos -from the background.Cherenkov light is sub-dominant in Borexino, but it is emitted almost instantaneously with respect to the slower scintillation light.Consequently, directional information is contained in the first hits of an event (after correcting for the time-of-flight of each photon).The CID analysis is based on the cos α observable: for a given PMT hit in an event, α is the aperture angle between the Sun and the hit PMT at the reconstructed position of the event (see also Fig. 3 in [20]).For background events, the cos α distribution is nearly uniform regardless which hit is considered.For solar neutrino events, the cos α distribution is flat for scintillation photons which are emitted isotropically, but has a characteristic non flat distribution peaked at cos α ∼ 0.7 for Cherenkov hits correlated with the position of the Sun.Since we cannot distinguish Cherenkov and scintillation photons, in our previous work [21] we have used only the 1 st and 2 nd hits of each event, which have the largest probability of being Cherenkov hits.In the new analysis presented in this paper we fully exploit the directional information contained in the first several hits.This choice is supported by Monte Carlo simulations and sensitivity studies as discussed in Sec.3.2.
The solar neutrino signal is obtained by fitting the cos α distributions of the selected first several hits, as a sum of signal and backgrounds contributions.The expected cos α distributions for signal and background are obtained from Monte Carlo simulations.As in Ref. [20], for each selected data event we simulate 200 MC events of solar neutrinos (represented by 7 Be or pep according to RoI, see below) and the same amount of the background events (represented by 210 Bi).These events are simulated with the same astronomical time as the data event and with the position smeared around the reconstructed vertex.From the fit we then obtain the total number of solar neutrinos N ν detected in the RoI.The fit also includes two nuisance parameters.The effective Cherenkov group velocity correction gv ch nuisance parameter accounts for small differences in the relative hit time distribution between scintillation and Cherenkov hits in data relative to the MC.The second parameter is the event position mis-reconstruction in the initial electron direction ∆r dir , an indirect effect of the Cherenkov hits, where the reconstructed position is slightly biased towards early hit PMTs of the corresponding event.Here ∆r dir is a free parameter of the fit, while gv ch is obtained independently and is constrained in the fit.
Compared to the previous proof-of-principle analysis [21,20], the current CID analysis has been improved in a variety of ways.The full detector live time can be used now thanks to a novel gv ch correction calibration.In the previous publication, gv ch has been obtained using 40 K γ calibration data (see [20]).In this work instead we calibrate  [2].The grey band shows the 7 Be-ν edge region used for the estimation of gv ch correction (RoI gvc ), while the CNO region used to measure the CNO-ν rate is shown in yellow band (RoI CNO ).gv ch by exploiting the 7 Be solar neutrino events which allows us to extend the analysis to the full Borexino dataset, as explained in Sec.3.1.Additionally, indirect Cherenkov information from the systematic influence on the vertex reconstruction and consequently on the cos α distribution was exploited by the inclusion in the analysis of later hits with negligible contribution of Cherenkov photons, see Sec. 3.2.Technical details of the CID fitting procedure can be found in Sec.3.3, while different systematic effects are discussed in Sec.3.4.
The final CID results regarding the CNO measurement are reported in Sec.3.5.

CID strategy for the CNO measurement with the full dataset
In the previous Borexino works [21,20], the calibration of the Cherenkov light group velocity gv ch has been performed using γ sources deployed during the Borexino calibration campaign in 2009 [20].The solar neutrino analysis was performed on the Phase-I dataset, that has been taken close in time to the source calibration of the detector [25].The gv ch found in this way was used to obtain the first measurement of 7 Be solar neutrinos with the CID method [21].For the CID measurement of CNO in this paper, the entire Borexino dataset is used (from 2007 until 2021).
Since the sub-nanosecond stability of the detector time response cannot be guaranteed for long periods, and no more calibrations have been performed after 2009, we developed a method to calibrate gv ch on the 7 Be shoulder data.This is done by using the same RoI as in [21,20] (here called RoI gvc electron equivalent energy range of 0.5 MeV ≲ T e ≲ 0.8 MeV) and performing the CID analysis where the 7 Be is constrained to the Standard Model predictions [2].The gv ch correction extracted in this way is then used in the CID analysis of the RoI CNO , in which the CNO contribution is maximized, and which is fully independent from RoI gvc .This step has been found to be justified according to MC studies, as the wavelength distribution of the detected Cherenkov photons produced by electrons from RoI gvc and RoI CNO is the same.
With this new strategy, the Cherenkov light gv ch can be calibrated on the same data-taking period as the one used for the CNO analysis.Two analyses have been performed in parallel for Phase-I (May 2007 to May 2010, 740.7 days) and Phase-II+III (December 2011 to October 2021, 2888.0 days).The gv ch correction obtained for Phase-I can be compared to the one previously obtained from the 40 K γ source [20].Additionally, the analyses of the two independent data-sets allows for the investigation of any variation of the detector response over time.The RoI gvc and the RoI CNO are shown for the Phase-II+III dataset in Fig. 2 for illustrative purposes.The results on gv ch are provided in Sec.3.5.1.
In the final analysis, the Three-Fold-Coincidence algorithm [4] is applied to the RoI CNO to suppress the cosmogenic 11 C background, preserving the exposure with a signal survival fraction of 55.77% ± 0.02% for Phase-I and 63.97% ± 0.02% for Phase-II+III.The radial (R FV ) and T e energy cuts of RoI CNO were optimized considering the expected number of solar neutrinos over the statistical uncertainty of the total number of events.The optimized cuts are R FV < 3.05 (2.95) m and 0.85 (0.85) MeV < T e < 1.3 (1.29) MeV) for the Phase-I (Phase-II+III).In addition, all other cuts including the muon veto and data quality cuts have been applied as in Refs.[4,19].The overall exposures for the CID CNO analysis are 740.7 days × 104.3 tonnes × 55.77% for Phase-I and 2888.0 days × 94.4 tonnes × 63.97% for Phase-II+III.The total exposure of Phase-II+III (477.81 years × tonnes) is about four times larger than that of Phase-I (118.04 years × tonnes).

N th -hit analysis approach
As mentioned above, the CID analysis is performed on the first several early hits of ToF corrected hit times from each event in the RoI.In this subsection we describe the optimization of the number of hits from each event to be used in the CID analysis.
The procedure is based on the comparison of the MC-produced cos α distributions of signal and background.
First, PMT hits of each individual event are sorted according to their ToF-corrected hit times and are labeled in this order as "N th -hit", with N = 1, 2, ... up to the total number of hits.Second, the cos α distributions are constructed for each N th -hit for both the signal and background MC.Third, for each N th -hit cos α distribution a number of 10,000 toy MC samples are simulated with the number of events observed as in the real data.Next, we perform a direct signal to background comparison, based on a standard χ 2 -test.Figure 3 shows the resulting ∆χ 2 for Phase-II+III in the RoI CNO averaged over the 10,000 toy datasets as a function of N th -hit.A larger average ∆χ 2 corresponds to a greater difference between the MC signal and background and thus a larger expected sensitivity for the CID fit, independent of the true signal to background ratio.While only the earliest ∼4 N th -hits have a relevant, direct contribution of Cherenkov hits to the cos α distribution of the neutrino signal, later N th -hits also contribute to the CID sensitivity due to the indirect Cherenkov influence on ∆r dir .A more in-depth explanation of the ∆r dir effect is shown in Fig. 9 in [20].
A possible impact of the gv ch and ∆r dir nuisance parameters on the N th -hit selection has been investigated and is presented in Fig. 3.The first hits of the events provide the largest ∆χ 2 values thanks to the direct Cherenkov light.A decrease of gv ch is decreasing the group velocity of Cherenkov photons and thus their contribution at early hits.The impact of ∆r dir can be seen for N th -hit > 4, where the contribution of direct Cherenkov hits becomes negligible relative to the scintillation hits, but the signal and background MC cos α histograms are still different from each other (∆χ 2 > 0).
In conclusion, the early hits selection for the CID analysis in both RoI gvc and RoI CNO is then performed from the first hit up to the N th -hit(max) = 15, 17 for Phase-I and Phase-II+III, respectively.This is an optimization where all direct and indirect Cherenkov information is used, while at the same time this selection keeps the contribution from delayed scintillation photons, undergoing various optical process during the propagation through the detector, relatively small.
The Cherenkov-to-scintillation photon ratio as a function of N th -hit has also been checked explicitly, as is shown in Fig. 4 for RoI CNO .As expected, it can be seen that the early N th -hits benefit from the largest Cherenkov-to-scintillation ratio of ∼ 13% for the first hit.
The overall total Cherenkov-to-scintillation ratio is small and found to be 0.475% in the MC.

CID fit procedure
The fitting strategy follows the procedure developed in our previous CID analysis [20].The data cos α distributions from the selected RoI, constructed for each N th -hit from the first up to the N th -hit(max), are fitted simultaneously with the MC produced, expected cos α distributions of the neutrino signal and background, where the signal cos α distribution depends on gv ch and ∆r dir .The nuisance parameter ∆r dir cannot be calibrated in Borexino and is left free to vary without a dedicated pull term.The number of cos α histogram bins used in the analyses is i = 60 for all energy regions and phases, as values of i < 30 reduce the expected CID sensitivity.

Fit in the RoI gvc
The CID analysis in RoI gvc used for the gv ch calibration is based on the χ 2 -test: where D n i and M n i are the numbers of cos α histogram entries at bin i for a given N th -hit n, for data and MC, respectively.The term N is the scaling factor between the MC and the data event statistics and the term "N 2 • M n i " in the denominator takes into account the finite statistics of MC.The explicit dependence of the fit on N ν , gv ch , and ∆r dir can be expressed by decomposing the MC contribution to the one from the signal S and the background B: The number of neutrino events N ν and ∆r dir are treated as nuisance parameters to produce the χ 2 (gv ch ) profile, where N ν is constrained by the SSM expectation.For this, the neutrino prior probability distribution P(N ν ) is given by the sum of the Gaussian probability distributions with mean and sigma from the high-metallicity (HZ) SSM and low-metallicity (LZ) SSM [2] predictions on the number of 7 Be+pep-ν in 7 Be-ν shoulder region, which is then convoluted with a uniform distribution of CNO-ν between zero and the HZ-SSM CNO expectation + 5σ.In this way, by leaving CNO reasonably free to vary, we avoid a potential correlation of the gv ch calibration and the subsequent measurement of the CNO-ν rate using this gv ch constraint.

Fit in the RoI CNO
The χ 2 -test for the measurement of number of solar neutrinos (N ν ) in RoI CNO is: using gv ch and ∆r dir as nuisance parameters.The gv ch parameter is now constrained by the previous calibration in RoI gvc through the pull term ∆χ 2 gv ch gv ch .

Systematic uncertainties
In this work we have performed a detailed evaluation of the systematic uncertainties.The quantitative evaluation of the systematic uncertainties on the CID analysis results are given in Sec.3.5.
It has been found that the choice of N th -hit(max) and the histogram binning do not introduce any systematic uncertainty.Backgrounds different than 210 Bi also contribute in the analysed energy intervals, but have been found to be indistinguishable in the CID analysis and do not contribute to the systematic uncertainty budget.Even if the external γ events are not uniformly distributed in the FV, due to their attenuation in the LS, the difference in the cos α distribution between these events and uniform background events is found to be safely negligible, given the statistics of the data.
The following effects have an impact on the final results: PMT selection Some PMTs feature an intrinsic misbehaviour of their hit time distribution, if compared to all other PMTs.They are identified using an enriched sample of 11 C events.In addition, a small number of PMTs feature inconsistency between the data and MC in the relative contribution to the first hits.The systematic effect has been evaluated by varying the selection of usable PMT.
Relative PMT hit time correction It has been found by analyzing an enriched sample of 11 C events, that the data PMTs have a small relative time offset between each other with an average value of 0.3 ns.This time offset has been measured with an uncertainty of up to ±0.1 ns.It is reasonable to correct this relative offset in data, as it does not exists in MC.This makes it also necessary to propagate the uncertainty of the PMT time correction through the entire analysis chain, which introduces a systematic uncertainty.
Influence of low number of signal events As described above, signal and background MC are produced on an event-by-event basis.This could introduce an additional systematic uncertainty through the particular choice of the event positions and neutrino directions used for the production of the signal MC.This systematic uncertainty is estimated by producing a large number of the signal MC cos α distributions corresponding to a random selection of the expected number of signal events from data in each phase.The data is then analyzed again with these reduced signal MC cos α histograms.This effect contributes to the systematic uncertainty only in the RoI CNO for Phase-I.

CNO-ν and pep-ν cos α distributions
The CNO-ν and pep-ν events show a significantly different energy distribution in the selected RoI.The expected Cherenkov to scintillation hits ratio for pep-ν events (0.475%) is higher than for CNO events (0.469%) due to their different energy distribution in RoI CNO .
Moreover, the angular distribution of recoiled electrons by CNO-ν and pep-ν is also different due to its dependence on energy distributions of the neutrino and recoiled electron.The final fit on the number of solar neutrinos is performed with the pep-ν MC and the systematic uncertainty is estimated by performing the CID analysis again with the CNO-ν MC.The absolute difference between the two analyses is used conservatively as the systematic uncertainty.This systematic is negligible for the gv ch calibration, as CNO+pep neutrinos are sub-dominant to the 7 Be neutrinos.

Constraint on non-CNO neutrinos
The CID analysis is only sensitive to the measurement of the total number of solar neutrinos N ν and cannot differentiate between different solar neutrino species.
Therefore, a measurement of the number of CNO neutrinos depends on the subtraction of the non-CNO neutrinos from N ν obtained in the RoI CNO .The number of pep neutrinos is constrained by the SSM predictions [2], while 8 B is constrained using the high precision flux measurement of Super-Kamiokande [6].

Exposure
In this category we cover systematic uncertainties on the determination of the fiducial mass and on the fraction of solar neutrinos in the RoI.These uncertainties are estimated using toy-MC studies, based on the results of the calibration campaign [25] on the performance of the position reconstruction and on the uncertainty on the energy scale described in [4], respectively.The uncertainty on the fiducial mass includes also the uncertainty on the scintillator density.Additionally, in the gv ch calibration using the constraint on the expected number of all solar neutrino events, we consider also the uncertainty on α/β discrimination applied to suppress 210 Po α decays in the RoI gvc .

Effective gv ch calibration on the 7 Be edge
The effective calibration of the Cherenkov light as a results of the CID analysis on the 7 Be edge, using Eq. 1, is gv ch = (0.140 ± 0.029) ns m −1 for Phase-I and gv ch = (0.089 ± 0.019) ns m −1 for Phase-II+III, including the systematic uncertainties summarized in Table I.The compatibility between the data and the MC model, illustrated in Fig. 5, is good with χ 2 /ndf = 874.9/897,p value = 0.70 for Phase-I and χ 2 /ndf = 1036.2/1017,p value = 0.33 for Phase-II+III.The χ 2 /ndf and p values have also been investigated for the individual N th -hits cos α histograms, with different binning choices to investigate the fit performance.For all cases, the best fit MC model cos α distribution is always in agreement with the data.Figure 5 shows an illustration of the best fit results (red) relative to a pure background hypothesis (blue).
For early hits, direct Cherenkov light causes the peak around cos α ∼ 0.7 and the influence of ∆r dir induces the negative slope for cos α < 0. For later N th -hits the Cherenkov peak washes away, but the indirect impact of the Cherenkov hits on the position reconstruction bias ∆r dir makes it still possible to distinguish between the neutrino signal and the background.The non-flat background cos α distribution originates from the live PMT non-isotropic distribution relative to the position distribution of the Sun around Borexino.
These gv ch values for Phase-I and Phase-II+III differ by less than 1.5σ and both are in agreement with the previous calibration performed at the end of Phase-I using a 40 K γ source: gv ch = (0.108 ± 0.039) ns m −1 (Fig. 13 in [20]).

CNO measurement with CID
This section describes the results of the measurement of CNO solar neutrinos with CID using the full Borexino live Figure 5.The CID data (black) and the best fit results (red) for the measurement of the gv ch parameter.While the analysis is done separately for Phase-I and Phase-II+III, here the sum of Phase-I+II+III is shown for illustration purposes.There are in total 78632 events in the RoI gvc .The best fit for the constrained number of neutrino events is N ν = 50063, while the best fit values for the parameter of interest are gv ch = 0.140 ns m −1 for Phase-I and gv ch = 0.089 ns m −1 for Phase-II+III.For comparison, the background MC (blue) scaled to the same total number of events is shown.(a) The sum of the first to fourth N th -hits cos α histograms shows the Cherenkov peak.(b) The sum of the fifth to the N th -hit(max) cos α histograms shows the effect the ∆r dir parameter on the later hits.time from May 2007 to October 2021.The gv ch values presented in Sec.3.5.1 are used as independent pull terms in the Eq. 3 for the fit in RoI CNO of their respective phases.This takes into account the potential systematic differences of the detector response between Phase-I and Phase-II+III.The resulting number of solar neutrino events N ν in the RoI CNO can be converted into the number of CNO neutrinos detected in the same energy region after constraining the contributions from pep and 8 B neutrinos, but without any a-priori knowledge of the backgrounds.This number of CNO events can be further transformed into the measurement of the CNO-ν interaction rate in Borexino and the CNO    betwwen the data and the MC model is good with χ 2 /ndf = 884.8/897,p value = 0.61 for Phase-I and χ 2 /ndf = 1000.7/1017,p value = 0.64 for Phase-II+III.The MC model is able to reproduce the data cos α distribution, which has also been investigated for the individual N th -hits cos α histograms.Figure 6 illustrates the best fit results (red) relative to a pure background hypothesis (blue), in which the CID cos α histograms of data (black) are shown for the sum of Phase-I + Phase-II+III, as well as for the sum of the early first to fourth N th -hits (top) and the sum of the later N th -hits from the fifth to N th − hit(max) (bottom).The actual fit is performed on Phase-I and Phase-II+III independently.The same observations made for Fig. 5 hold also true for Fig. 6, where the early hits show the Cherenkov peak and the later hits show the impact of ∆r dir .

Fit response bias correction
Toy-MC analyses found that the fit of the number of solar neutrinos in RoI CNO shows a small systematic shift between the injected number of neutrinos and the best fit number of neutrinos, due to the correlation between the nuisance parameters (gv ch , ∆r dir ) and the relatively low total number of neutrino events.This fit response bias is induced by the two nuisance parameters as they only impact the shape of the neutrino signal MC cos α distribution but not that of background.We note that this effect was found to be negligible in RoI gvc due to the relative large number of signal events and the large signal to background ratio.
The value of the fit response bias in RoI CNO is estimated Figure 6.Illustration of the CID data (black) and the best fit results (red) summed for the Phase-I + Phase-II+III with a total of 8964 events in the RoI CNO .The best fit of the total number of neutrino events is N ν = 3519 without any systematic correction.For comparison, the background MC (blue) scaled to the same total number of events is shown.(a) The sum of the first to fourth N th -hits cos α histograms shows the Cherenkov peak.(b) The sum of the fifth to the N th -hit(max) cos α histograms shows the effect the ∆r dir parameter on these later hits.using the Bayesian posterior distribution of N ν [27], which is produced through a toy-MC rejection sampling, described in summary below.The prior distribution for the number of neutrino events is chosen to be uniform between zero and the number of selected data events (2990 for Phase-I and 5974 for Phase-II+III), the prior distribution of ∆r dir is also uniform, and the prior distribution of gv ch is given by the measurement at the 7 Be-ν edge RoI gvc P gv ch ∝ exp − 1 2 ∆χ 2 (gv ch ) .The pseudo-data inputs N sim ν , gv sim ch , ∆r sim dir are sampled from the MC signal and background cos α distributions following these model parameter prior distributions.The analysis is then performed in the same way as for the real data and results in best fit values of the pseudo-data N fit ν , gv fit ch , ∆r fit dir .The real data result now defines a multivariate Gaussian distribution P accept (N ν , gv ch , ∆r dir ) with a mean value given by its best fit values and with a standard deviation given by the systematic uncertainty of the PMT time corrections.The sampled true values of the triplet N sim ν , gv sim ch , ∆r sim dir are then saved only with a probability of P accept (N fit ν , gv fit ch , ∆r fit dir ), given by the best fit result of the pseudo-data, otherwise they are rejected.The resulting distributions of the true values for N ν , gv ch , ∆r dir then correspond to their Bayesian posterior distributions.
The fit response bias is illustrated in Fig. 7 for Phase II+III.The likelihood distribution P(N ν ) ∝ exp − 1 2 ∆χ 2 (N ν given by the χ 2 fit of data with Eq. 3 and averaged over the 1000 fits with different PMT time offsets is shown in blue.The black distribution is given by the simulation of 20k pseudo-data analyses, selected through the rejection sampling MC described above.The red distribution is produced by shifting P(N ν ) by a value of ∆N ν = −109 ± 4 and this distribution is well in agreement with the black rejection sampled distribution.It is therefore used as the posterior distribution of the CID analyses.We note that for the Phase-I the situation is similar and the shift is found to be −50 ± 4 events.

Inclusion of systematics
The final result of the CID analysis for the number of solar neutrinos is given by the Bayesian posterior distribution of N ν , marginalized over the nuisance parameters and convoluted with the systematic uncertainties.The relevant systematic uncertainties are shown in Table II and assumed to be normally distributed.

CID results on CNO
The interpretation of the CID results requires the correct treatment of the physical boundaries of the analysis, i.e. 0 ≤ N ν ≤ 2990 (5974) for Phase-I (Phase-II+III), respectively.This is done in a Bayesian interpretation, based on the posterior distribution P (N ν ) shown in Fig. 8. Next, the distribution of the number of CNO-ν events is estimated by constraining the expected number of pep and 8 B neutrino events (N pep+ 8 B ) where the constraint on the number of pep neutrinos uses the SSM predictions [2] and 8 B is constrained using the high precision flux measurement of Super-Kamiokande [6] including model uncertainties, the difference between HZ-SSM and LZ-SSM predictions, as well as the Borexino FV and energy systematic uncertainties from Table II.This is done through the convolution of the N ν posterior distributions from Fig. 8 with the predicted P(N pep+ 8 B ) probability distribution: P(N CNO ) = P(N ν ) * P(−N pep+ 8 B ).The resulting P(N CNO ) posterior distributions are shown in Fig. 9.
The CID measurement for the number of CNO-ν events is then N CNO = 270  uncertainty corresponds to the equal-tail 68% CI within the physical boundaries, including all systematics.
It has been observed that Phase-I and Phase-II+III do not show prohibitively different behavior for the full CID analysis-chain and the MC model is well able to reproduce the data cos α histograms for both phase selections and for each selected energy region.It is then reasonable to combine the conditionally independent results of Phase-I and Phase-II+III, through the convolution of both posterior distributions P (N CNO ) I+II+III = P (N CNO ) I * P (N CNO ) II+III .The probability that exactly zero CNO-ν events contribute to the measured data CID cos α distribution is P(N CNO = 0) = 1.35 × 10 −3 for Phase-I, P(N CNO = 0) = 5.87 × 10 −5 for Phase-II+III, and P(N CNO = 0) = 7.93 × 10 −8 for the combined result.This corresponds to a one-sided exclusion of the zero-CNO hypothesis at about 5.3σ credible level for the combination of Phase-I and Phase-II+III.
The CNO-ν rate probability density function is calculated from the measured posterior distribution of CNO-ν events, using the exposure of the respective phases.The effective exposure is given by the product of the fiducial mass, the detector live time, the TFC-exposure, the trigger efficiency, and the fraction of CNO-ν events within the selected energy region.The final CID result for the CNO-ν rate, using the full dataset of Phase-I + Phase-II+III, is R CID CNO = 7.2 ± 2.5 (stat) ± 0.4 sys +1.1 −0.8 (nuisance) cpd/100 tonnes = 7.2 +2.8 −2.7 cpd/100 tonnes The quoted uncertainties now also show the systematic uncertainties from Table II separately from the influence of the nuisance parameters gv ch and ∆r dir .The quoted statistical uncertainty corresponds to a hypothetical, perfect calibration of these CID nuisance parameters.The results are summarized in Table III.

Combined CID and multivariate analysis
In this Section we combine the CID analysis with the standard multivariate fit of Phase-III to improve the result on CNO neutrinos.This is done by including the posterior distributions of solar neutrinos from the CID analysis, shown in Fig. 8, in the multivariate analysis likelihood.
By statistically subtracting the sub-dominant 8 B neutrinos contribution and converting N CNO+pep to the corresponding interaction rate, it is possible to use these posterior distributions as external likelihood terms in the minimization routine.Following this procedure, two multiplicative pull terms constraining the number of CNO and pep neutrino events are used: the first one is related to the Phase-I (L P−I CID ), while the second one refers to Phase-II+III datasets (L P−II+III CID ).The overall combined likelihood used for this analysis becomes: where the first three terms correspond to an improved version of the standard multivariate analysis described in [19].This approach couples the one-dimensional Poisson likelihood for the TFC-Tagged dataset with a two-dimensional one (energy and radius) for the TFC-subtracted, to enhance the separation between signal and backgrounds.We have improved the binning optimization and used an updated version of Monte Carlo.
The pep neutrinos interaction rate is constrained with 1.4% precision to the 2.74±0.04cpd/100 tonnes value, by combining the Standard Solar Model predictions [2], the most current flavor oscillation parameters set [29] and the solar neutrino data [30,31].This constraint is applied with the Gaussian pull term L pep .An upper limit on the 210 Bi rate of (10.8 ± 1.0 cpd/100 tonnes) is applied with the half Gaussian term L210 Bi .This upper limit is obtained from the rate of the 210 Bi daughter 210 Po (see [4,19] for more details).

Results
As in [19], the energy RoI for the multivariate analysis is 0.32 MeV < T e < 2.64 MeV for electron recoil kinetic energy.The reconstructed energy spectrum scale is quantified in the N h estimator, representing total number of detected hits for a given event (see Sec. 2).The dataset is the same one analyzed in [19], in which the exposure amounts to 1431.6 days × 71.3 tonnes.
Along with CNO solar neutrinos, the free parameters of the fit are divided into three categories: internal ( 85 Kr and 210 Po) and external ( 208 Tl, 214 Bi, and 40 K) backgrounds, cosmogenic backgrounds ( 11 C, 6 He, and 10 C), and solar neutrinos ( 7 Be).Since 8 B solar neutrinos exhibit a flat and marginal contribution, the corresponding interaction rate is fixed at high-metallicity expectations from Solar Standard Model.As discussed in Eq. 4, the interaction rates of pep neutrinos and 210 Bi background are constrained with likelihood pull terms, and CID results reported in Sec.3.5 are accounted for as additional external constraints.The result of the fit for the energy and radial projections is shown in Fig. 11.For both projections, the sum of the individual components from the fit (magenta) is superimposed on the data (grey points).CNO neutrinos, 210 Bi and pep neutrinos contributions are displayed in solid red, dashed blue and dotted green lines, respectively, while the other spectral components ( 7 Be and 8 B neutrinos, other backgrounds) are shown in grey.The analysis has been performed using N h as energy estimator and the conversion to keV energy scale was performed only for the plotting purposes.The radial fit components, that are the uniform and the external backgrounds contributions, are shown in solid blue and dashed grey lines respectively.
The multivariate fit returns an interaction rate of CNO neutrinos of 6.7 +1.2 −0.7 cpd/100 tonnes (statistical error only).The agreement between the model and data is quantified with a p value of 0.2.
To account for sources of systematic uncertainty, the same Monte Carlo method described in [4,19] has been adopted.In a nutshell, hundred thousands Monte Carlo pseudo-experiments were generated, including relevant effects able to introduce a systematic error, such as the energy response non linearity and non uniformity, the time variation of the scintillator light yield and the different theoretical models for the 210 Bi spectral shape.The analysis is performed on these pseudo datasets assuming the standard response to study how this impacts the final result on CNO, yielding a total systematic uncertainty of +0.31 −0.24 cpd/100 tonnes.Other sources of systematic error are included in the estimation of the upper limit on 210 Bi contamination, as discussed more in detail in [4,19].The negative log-likelihood profile as a function of the CNO rate is reported in Fig. 12.The solid and dashed black lines show the results with and without systematic uncertainty, respectively.The result without CID constraint reported in [19] is included (blue line), for comparison.The improvement is clear especially for the upper value of the CNO rate.The CNO interaction rate is extracted from the 68% quantile of the likelihood profile convoluted with the resulting systematic uncertainty, as R MV+CID (CNO) = 6.7 +1.2 −0.8 cpd/100 tonnes.The significance to the no-CNO hypothesis reaches about 8σ C.L., while the resulting CNO flux at Earth is Φ(CNO) = 6.7 +1.2 −0.8 × 10 8 cm −2 s −1 .Following the same procedure used in [19], we use this result together with the 8 B flux obtained from the global analysis of all solar data to determine the abundance of C + N with respect to H in the Sun with an improved precision, for which we find N CN = 5.81 +1.22  −0.94 × 10 −4 .This error includes both the statistical uncertainty due to the CNO measurement, and the systematic errors due to the additional contribution of the SSM inputs, to the 8 B flux measurement, and to the 13 N/ 15 O fluxes ratio.Similar to what was inferred from our previous publication, this result is in agreement with the High Metallicity measurements [14,15], and features a 2σ tension with Low Metallicity ones [16,17,18].Similarly, if we combine the new CNO result with the other Borexino results on 7 Be and 8 B in a frequentist hypothesis test based on a likelihood ratio test statistics we find that, assuming the HZ-SSM to be true, our data disfavours LZ-SSM at 3.2σ level.The blue line shows the same profile obtained in the previously published analysis without CID constraint [19].The blue, violet, and gray vertical bands show 68% confidence intervals (CI) for the low-metallicity SSM B16-AGSS09met (3.52 ± 0.52 cpd/100 tonnes) and the high-metallicity SSM B16-GS98 (4.92 ± 0.78 cpd/100 tonnes) predictions [2,28], and the new Borexino result including systematic uncertainty, respectively.

Conclusions
In this work we have presented the results on CNO solar neutrinos obtained using the "Correlated and Integrated Directionality" (CID) technique.
We have shown that the CID technique can be used to extract the CNO signal without any a priori assumptions on the backgrounds, in particular that of 210 Bi.The Phase-I (May 2007 to May 2010, 740.7 days) and Phase-II+III (December 2011 to October 2021, 2888.0 days) datasets have been analyzed independently to investigate possible variations of the detector response over time.By adopting the Bayesian statistics, we have combined the conditionally independent results of Phase-I and Phase-II+III: the resulting CNO rate obtained with CID only is 7.2 +2.8 −2.7 cpd/100 tonnes.The no-CNO hypothesis including the pep constraint only is rejected at 5.3σ level.This result, albeit less precise than the one published by Borexino using the standard multivariate analysis, is the first obtained without the application of a 210 Bi constraint.
We have also obtained an improved CNO solar neutrino result by combining the standard multivariate analysis with the CID technique.The CID technique helps in separating the solar signal from non solar backgrounds, improving the significance and precision of the CNO measurement with respect to the result previously published by Borexino.The resulting CNO interaction rate is 6.7 +1.2 −0.8 cpd/100 tonnes and the significance against the absence of a CNO signal, considered as the null hypothesis, is about 8σ.The C+N abundance with respect to H is calculated from this result following the procedure adopted in [19] and is found to be N CN = 5.81 +1.22  −0.94 × 10 −4 , compatible with the SSM-HZ metallicity measurements.
In conclusion, we have shown that the directional information of the Cherenkov radiation can be effectively combined with the spectral information coming from scintillation, for solar neutrino studies.This combined detection approach provides a measurement that is more powerful than the individual methods on their own.The sensitivity of the CID method could be significantly improved in future liquid scintillator-based detectors by optimizing the Cherenkov-to-scintillation ratio and by performing dedicated calibrations campaigns.

8 Figure 2 .
Figure2.Illustration of the two RoIs used in the analysis on the energy spectrum of the Phase-II+III data in a fiducial volume of 2.95 m radius.Monte Carlo PDFs of different solar neutrino components are scaled to high-metallicity SSM prediction[2].The grey band shows the7 Be-ν edge region used for the estimation of gv ch correction (RoI gvc ), while the CNO region used to measure the CNO-ν rate is shown in yellow band (RoI CNO ).

Figure 3 .
Figure 3. ∆χ 2 between the Phase-II+III RoI CNO neutrino signal and background MC cos α distributions for different selections of nuisance parameters.∆r dir = 2.7 cm corresponds to the nominal value observed in the neutrino MC.

Figure 4 .
Figure 4. Cherenkov-to-scintillation PMT hit ratio as a function of the time-of-flight sorted N th -hits for the pep neutrino Monte Carlo of Phase II+III in the RoI CNO (0.85 MeV -1.3 MeV).

Figure 7 .
Figure 7. Illustration of the fit response bias for Phase II+III in the RoI CNO .The data fit result, i.e. the likelihood P(ν) given by the ∆χ 2 profile of Eq. 3 and averaged over 1000 fits with different PMT time offsets is shown in blue.The posterior distribution of 20k pseudo-data analyses, selected through rejection sampling, is shown in black.The red line corresponds to the posterior distribution that includes only the systematics from the PMT time alignment correction.

Figure 8 .
Figure 8.The CID measured posterior probability distributions for the number of solar neutrinos N ν in the RoI CNO for Phase-I (blue) and Phase-II+III (red).All systematic effects are included.

Figure 9 .
Figure 9.The CID measured posterior probability for the number of CNO-ν events after constraining pep and 8 B neutrinos for Phase-I (blue) and Phase-II+III (red).All systematic effects are included.

Figure 11 .
Figure 11.Multivariate fit results for the TFC-subtracted dataset, projected over the energy (top panel) and the radius (bottom panel) dimensions.For both projections, the sum of the individual components from the fit (magenta) is superimposed on the data (grey points).CNO neutrinos,210 Bi and pep neutrinos contributions are displayed in solid red, dashed blue and dotted green lines, respectively, while the other spectral components ( 7 Be and 8 B neutrinos, other backgrounds) are shown in grey.The analysis has been performed using N h as energy estimator and the conversion to keV energy scale was performed only for the plotting purposes.The radial fit components, that are the uniform and the external backgrounds contributions, are shown in solid blue and dashed grey lines respectively.

Figure 12 .
Figure 12.CNO-ν rate negative log-likelihood (−2∆ ln L) profile obtained from the 2-dimensional multivariate spectral fit combined with the CID analysis constraint, with and without folding in the systematic uncertainties (black dashed and solid lines respectively).The blue line shows the same profile obtained in the previously published analysis without CID constraint[19].The blue, violet, and gray vertical bands show 68% confidence intervals (CI) for the low-metallicity SSM B16-AGSS09met (3.52 ± 0.52 cpd/100 tonnes) and the high-metallicity SSM B16-GS98 (4.92 ± 0.78 cpd/100 tonnes) predictions[2,28], and the new Borexino result including systematic uncertainty, respectively.

TABLE I .
Systematic uncertainties of the gv ch measurement in the RoI gvc , relative to the best fit value.

TABLE II .
Systematic uncertainties on the number of solar neutrino events N ν in RoI CNO , relative to the best fit value.The uncertainty from pep+ 8 B-ν constraint is relevant only for N CNO .The last two rows are relevant only for the CNO-ν rate (R CNO ) calculation.
1.1 σ) is 1.7 times less likely to be true, given the results of the Borexino CID analysis.

TABLE III .
CID CNO-ν results with systematic uncertainties.