Borexino results on neutrinos from the Sun and Earth

Borexino is a 280-ton liquid scintillator detector located at the Laboratori Nazionali del Gran Sasso in Italy. Since the start of its data taking in May 2007, it has provided several measurements of low-energy neutrinos from various sources. At the base of its success, lie unprecedented levels of radio-purity and extensive thermal stabilisation, both resulting from a years-long effort of the collaboration. Solar neutrinos, emitted in the hydrogen-to-helium fusion in the solar core, are important for the understanding of our star, as well as neutrino properties. Borexino is the only experiment that has performed a complete spectroscopy of the \emph{pp} chain solar neutrinos (with the exception of the \emph{hep} neutrinos contributing to the total flux at $10^{-5}$ level), through the detection of \emph{pp}, $^7$Be, \emph{pep}, and $^8$B solar neutrinos and has experimentally confirmed the existence of the CNO fusion cycle in the Sun. Borexino has also detected geoneutrinos, antineutrinos from the decays of long-lived radioactive elements inside the Earth, that can be exploited as a new and unique tool to study our planet. This paper reviews the most recent Borexino results on solar and geoneutrinos, from highlighting the key elements of the analyses up to the discussion and interpretation of the results for neutrino, solar, and geophysics.


Introduction
Neutrinos have a special position within the Standard Model (SM) of elementary particles. Having no electric charge and no colour, they interact with matter only via weak force. Thus, the probability of their interaction is very low, and this fact has a two-fold consequence. On a positive side, neutrinos bring us unperturbed information from otherwise inaccessible locations-including the solar core and the Earth's interior. On the other hand, this very same characteristic makes neutrino detection an experimental challenge. Large-volume detectors, constructed from specially developed materials with radioactivity levels of many orders of magnitude below the ambient values, must be placed in underground laboratories to shield them from cosmic radiation. Such measures are necessary in order to acquire high statistics of neutrino events with acceptable signal-to-background ratio. An enormous years-long effort is behind each neutrino experiment.
Borexino, a 280 ton liquid scintillator detector located at the Laboratori Nazionali del Gran Sasso (LNGS) in Italy, became an "ideal scenario" for future projects regarding the achieved radio-purity. Many years dedicated to the selection of construction materials and their treatment, involving surface cleaning or purification of the liquid scintillator, preceded the start of the data taking in May 2007. Section 2 of this paper is dedicated to the description of the basic features of the experiment and its main phases, defined according to the milestones achieved in terms of the radio-purity of the liquid scintillator and detector's thermal stability. The described arguments range from the detector setup (Section 2.1), through the event reconstruction (Section 2.2), calibration and Monte Carlo simulation (Section 2.3), background contamination (Section 2.4), up to the interactions used to detect neutrinos and antineutrinos in Borexino (Section 2.5).
This paper reviews the most recent Borexino results on neutrinos from the Sun (Section 3) and Earth (Section 4). Each section introduces the field and importance of the measurement of the respective neutrinos, summarizes the details of the analysis and describes the latest results, including their interpretation and discussion.
Solar neutrino production and propagation from the solar core to our detector are discussed in Section 3.1. In particular, the fusion of Hydrogen to Helium proceeding via the principal pp chain and the subdominant CNO cycle are covered in Section 3.1.1. Section 3.1.2 instead is dedicated to the Standard Solar Model (SSM) that predicts solar neutrino fluxes as a function of the not-yet-well-understood abundances of elements heavier than Helium ("solar metallicity problem"). The neutrino-energy-dependent process that converts part of the solar neutrino flux from a pure electron flavour to a mixture of all flavours is discussed in Section 3.2.1. Common features of all Borexino solar neutrino analyses are synthesised in Section 3.2. The comprehensive spectroscopy of the pp chain solar neutrinos based on [1] is contained in Section 3.3: the strategy of this particular analysis (Section 3.3.1), the results on pp, pp, 7 Be, and 8 B neutrinos (Section 3.3.2) and their implications for solar and neutrino physics (Section 3.3.3). These include: the confirmation of the origin of solar energy and of the Sun's thermal stability over the scale of 100 000 years, a mild preference towards the high-metallicity SSM with respect to the low-metallicity counterpart, and observation of the energy dependent electron neutrino survival probability for different neutrino species, excluding purely vacuum-dominated flavour conversion. The seasonal variation of the solar neutrino flux, expected due to the eccentricity of the Earth's orbit around the Sun, has been observed for 7 Be neutrinos [2], as briefly discussed in Section 3.3.4. About 1% of the solar energy is expected to be produced in the CNO cycle, a fact confirmed by the recent Borexino discovery [3]. Section 3.4 is dedicated to this observation of solar neutrinos from the CNO fusion with the analysis strategy summarised in Sections 3.4.1 and 3.4.2. The expected sensitivity ( [4] and Section 3.4.3) is in agreement with the results (Section 3.4.4) obtained from both a spectral fit and a counting analysis. Precise measurement of solar neutrinos can lead to constraining of the parameter space of suggested non-standard neutrino interactions (Section 3.5.1 based on [5]) or placing strong limits on the effective neutrino magnetic moment (Section 3.5.2 based on [6]).
The part of this paper dedicated to geoneutrinos is introduced by providing essential information about the Earth's structure and heat budget (Section 4.1.1), Bulk Silicate Earth models (Section 4.1.2) that predict the composition of the primitive Earth and consequently also the abundances of 238 U and 232 Th, which produce geoneutrinos, a new tool to probe our planet (Section 4.1.3). The latest Borexino geoneutrino analysis [7] is discussed in Section 4.2, including the expected levels of geoneutrino, and other antineutrino signals, as well as non-antineutrino backgrounds, event selection cuts, spectral fit, and sources of systematic errors in Sections 4.2.1 to 4.2.6, respectively. The results and their interpretation in terms of measured geoneutrino signal at LNGS, extracted geoneutrino signal from the Earth's mantle and the corresponding radiogenic heat, as well as imposed limits on the power of a hypothetical georeactor at different locations inside the Earth are discussed in Sections 4.3.1 to 4.3.4, respectively.
The final concluding remarks about the Borexino solar and geoneutrino analyses are given in Section 5.

Detector setup
The Borexino detector has an onion-like structure with the radio-purity of materials increasing towards the centre. A schematic of the detector is shown in Figure 1. The main neutrino target is 280 ton of liquid scinitillator (LS) located in the core of the detector. Pseudocumene (PC, 1,2,4-trimethylbenzene) is used as a solvent with 1.5 grams per litre of PPO (2,5-diphenyloxazole) as a solute. The scintillator density is (0.878 ± 0.004) g cm −3 [7], where the error considers the changes due to the temperature variations during the whole data-taking period. The scintillator is contained in a thin spherical nylon inner vessel (IV) with a radius of 4.25 m. The LS is surrounded by a non-scintillating buffer liquid (inner buffer) made up of PC doped with dimethylphthalate as a quencher. The shape of the IV changes with time, because of a leak of the LS from the IV to the buffer region which started around April 2008 [9,7]. The leak was identified by reconstructing a large number of events outside the IV. The IV shape is reconstructed based on contaminants on its surface selected between 0.8-0.9 MeV, dominated by external backgrounds such as 40 K, 214 Bi, and 208 Tl (see Section 2.4). The density of the inner buffer is almost the same as for the scintillator material. This buffer region is held by a nylon outer vessel (OV) with a radius of 5.50 m, followed by a second outer buffer region, which in turn is surrounded by a stainless steel sphere (SSS) with a radius of 6.85 m, which holds 2218 8-inch photomultiplier tubes (PMTs), facing inwards. The 2.6 m thick buffer region shields the inner volume against external radioactivity from the PMTs and the SSS. Moreover, the OV serves as a shielding against inward-diffusing Radon. The inner components contained inside the SSS are called the inner detector (ID). The SSS is enclosed in a cylindrical tank filled with high-purity water, additionally endowed with 208 external PMTs, which define the outer detector (OD). This water tank serves as an extra shielding against external gammas and neutrons, and as an active Cherenkov veto for residual cosmic muons passing through the detector.
Charged particles interacting with molecules of the LS produce scintillation light, the amount of which is roughly proportional to the deposited energy. The exact amount of the emitted light depends on the particle type. In addition, the energy scale is to some extent intrinsically non-linear due to the ionisation  Figure 1: Schematic view of the Borexino detector. From inside to outside: the liquid scintillator contained in the nylon inner vessel, nylon outer vessel, stainless steel sphere, and the water tank for the Cherenkov muon veto. From [1].
quenching [11] and emission of the Cherenkov light. Particles causing high ionisation density experience high levels of quenching. Due to this, the visible energy of αs' is quenched in the LS with a factor of about 10 with respect to electrons of the same energy.
The PMTs in Borexino convert the light to photoelectrons (p.e.), defined as the electrons removed from the photocathode of the PMT through incident photons. In Borexino, the effective light yield is about 500 p.e. per 1 MeV of electron equivalent. The PMTs have random dark-noise coincidences with a rate of about 1 hit per 1 µs in the whole detector. The average number of working PMT channels in the ID varies in time and is 1576 and 1238 for Phase-II and Phase-III, respectively [10,3].
The ID PMTs are coupled to an analog front-end (FE board, FEB) that amplifies the signal. Further processing differs for the two data acquisition (DAQ) systems [8]: the main DAQ and a semi-independent flash analog-to-digital converter (FADC) sub-system commissioned in November 2009. The two systems run independently and have different triggers. The main DAQ treats every PMT information separately, has a threshold of about 50 keV, and a read-out window of 16 µs. The FADC system instead processes an integrated signal of 24 FEB outputs, has an energy threshold of 1 MeV, and 1.28 µs read out window.
In the main DAQ, the FEB signal is fed to a digital circuit (Laben board, LB). A fast amplified timing signal can trigger a threshold discriminator, set to about 1% of the p.e., defining a PMT hit. The FEB also integrates the PMT current and provides the second input to the digital LB. This second signal provides the charge of the hit, integrated in 80 ns, that is proportional to the number of photons hitting the PMT during the integration time. The photons eventually falling on the same PMT in the time interval from 80 to 180 ns are not detected.

Event reconstruction
In this Section, we describe the event reconstruction algorithms applied to the triggers of the main DAQ system. Firstly, the clustering algorithm identifies an accumulation of hits in the 16 µs DAQ gate that correspond to a single physical event. Typically, there is just one cluster present, but multi-cluster triggers do exist. Higher level event reconstruction algorithms such as position, energy, and particle identification are then applied on each identified cluster. The main variables used in the solar and geoneutrino analyses presented in this paper are obtained from the data of the main DAQ system. The FADC system is optimised for multi-MeV events [12,13] and was also successfully used to improve the muon tagging efficiency and to identify noise events [7].

Position reconstruction
The position reconstruction is based on the time-of-flight (TOF) technique, that subtracts from each measured hit time t i a position-dependent time-of-flight T TOF i from the point r 0 of particle interaction at time t 0 to the PMT at position r i , where the i th hit was detected: Here, n LS = 1.68 is the effective refraction index of the LS determined using calibration data [14] and c 0 the speed of light in vacuum. The algorithm determines the most likely vertex (t 0 , r 0 ) of the interaction, using the arrival times t i of the detected hits on each PMT and the position vectors r i of the corresponding PMTs. The likelihood maximization uses the probability density functions (PDFs) of the hit detection as a function of the time elapsed from the emission of scintillation light due to the interaction of an electron. The shape of the PDFs changes according to the charge of each hit [9]. The position resolution is about 10 cm at 1 MeV at the centre of the detector [9]. For other positions with larger radii, the resolution decreases on average by a few centimetres.
Energy reconstruction The visible energy in Borexino is different for different particle types. The detected light is proportional to the deposited energy, up to the leading order. There are intrinsic non-linearities of the energy scale due to the particle dependent ionisation quenching [11] and a small contribution from Cherenkov radiation [10]. The deposited energy is parameterized via the following energy estimators: • N h : number of PMT hits, including multiple hits on a single PMT.
• N p : number of triggered PMTs, ignoring multiple hits on the same PMT.
• N dt 1(2) p : a variant of N p , restricting the considered time interval to 230(400) ns after the cluster start time.
• N pe : charge expressed in number of photoelectrons collected in all PMT hits.
The energy estimators are normalised to 2000 working PMTs, since not all PMTs are active during the datataking. In addition, the energy estimators can also be geometrically normalised, considering the relative weight of each PMT to be proportional to the solid angle with respect to the reconstructed position of the event. This normalisation takes into account the fact that the amount of light seen by each PMT depends on the distance of this PMT to the event. An electron with kinetic energy of 1 MeV produces approximately 500 photoelectrons in the Borexino detector. This results in 5%/ E (MeV) energy resolution.
α/β discrimination The discrimination of α and β particles is based on the different types of interactions of these particles. αs' have a high ionisation quenching compared to βs', leading to different hit time profiles, as shown for 214 Bi-βs' and 214 Po-αs' in Figure 2(a). In Borexino, α/β discrimination is performed on an event-by-event basis. The algorithm is tuned based on 222 Rn-correlated 214 Bi(β − )-214 Po(α) fast coincidences introduced into the detector by a small air leak during the water extraction cycles performed  214 Bi-βs' (in black) after TOF subtraction. From [9]. (b) The MLP variable, shown for the same events, characterised by a high separation power between α and β/γ like events. From [7].
between June 2010 and August 2011. The α/β discrimination variable currently used in Borexino is the so-called MultiLayer Perceptron (MLP), while other variables have also been used in the past, as described in [9]. The MLP variable is based on neural networks and, in Borexino, it consists of one input layer, one hidden layer, and an output layer. It can distinguish between two classes of events by training on the acquired data. Thus, it can exploit not only the time profiles but also the pulse shape variables such as mean times, variances, skewness, and kurtosis associated with the given training data sample [7]. The α particles in Borexino have an MLP value around 0, while the β/γ particles have a value around 1, as shown in Figure 2(b). β + /β − discrimination The hit time profiles of electrons and positrons look very similar in the liquid scintillator, making their separation very challenging. Therefore, the discrimination is done on a statistical basis and not on an event-by-event basis. Positrons emitted in the LS build ortho-positronium in 50% of the cases, as discussed in [15]. This formation leads to a delay of the e + e − → γγ annihilation process, with a typical lifetime of ∼3 ns. The lifetime of para-positronium is about 125 ps in vacuum, making its contribution indistinguishable from the prompt light emission caused by the positron [9,15]. The pulse shape discrimination of β − s' and β + s' in Borexino is based on the likelihood of the position reconstruction, normalised to the number of PMTs N p . This parameter is called PS-L PR . The position reconstruction is based on the emission profiles for electrons, as discussed in the Position reconstruction paragraph. While the spatial position reconstruction for β − s' and β + s' is precise within the position reconstruction resolution, the likelihood value is worse for positrons than for electrons. This causes the difference in the PS-L PR variable, making pulse shape discrimination possible. Figure 3 shows the β − and β + hit time distributions (Figure 3(a)) and the PS-L PR variable distributions (Figure 3(b)) in the region from 400 to 650 N h , corresponding to approximately 1.0 to 1.8 MeV. The PS-L PR variable for β − s' is taken from Monte Carlo simulations, while for β + s' it is taken from a pure sample of the 11 C(β + ) cosmogenic background, selected via the three-fold coincidence algorithm (see Sections 2.4 and 3.2).

Calibration and Monte Carlo simulation
In order to understand the whole Borexino detector and validate the physics model adopted for the description of the light emission, propagation, and detection by PMTs, a dedicated Geant4-based Monte Carlo (MC) simulation code has been developed by the Borexino collaboration [16]. The MC code has been tuned based on calibration of the detector with radioactive sources and various laboratory From [9]. (a) PS-L PR variable in the energy window from 400 to 650 N h . From [10].
measurements [14]. The Borexino calibration campaigns were performed in November 2008, January 2009, June-July 2009, and July 2010 using different types of radioactive sources as extensively discussed in [14] and [16].
Calibrations The goal of the calibration campaign was to (1) determine the accuracy and resolution of the position reconstruction, (2) measure the absolute energy scale and resolution, (3) estimate the energy response and non-uniformity depending on the energy and position of an event, (4) tune the MC simulation framework. The different calibration sources were used to study the responses of αs', βs', γs', and neutrons, covering an energy range of 0.122-10 MeV. These sources were deployed in 295 locations. The energy scale was determined through the usage of different monochromatic γ sources ranging from 0.122 to 2.615 MeV located at the centre of the detector and at few positions along the vertical detector axis. To study the uniformity of the trigger efficiency, some of these γ sources were also deployed at different positions and at larger distances from the centre. 222 Rn was used as a α/β-source, while 241 Am- 9 Be was used as a neutron-source (see Sections 3.2 and 4.2). The calibration sources used to calibrate the event position reconstruction were 222 Rn and 241 Am-9 Be, which have been placed in 182 and 29 positions in the scintillator, respectively [14]. The external calibration was done using 228 Th whose daughter nuclide 208 Tl is a strong gamma source. This helped in studying the exact determination of the IV shape, the external γ background near the IV, and the asymmetries in the energy response near the IV. The source positions have been measured with a charge coupled device (CCD) camera system 1 [14]. In total, 7 CCD cameras were mounted on the SSS. The differences (nominal-to-reconstructed) in the (x, y, z) coordinates have been determined with a precision better than 2 cm [9]. In addition to these dedicated calibration campaigns, there are constant offline checks of the detector's stability and regular online PMTs' charge and timing calibration [9], in order to monitor the quality of the acquired data.

Monte Carlo Simulations
The Borexino Monte-Carlo (MC) simulation [16] is able to simulate all the processes after the interaction of a particle in the detector, including the knowledge of detector effects. It is based on the GEANT4 software (v4.10.5). The code is able to generate the physical event and track the light production, propagation and detection. Furthermore, the electronics and the trigger responses are fully simulated. The framework monitors the detector evolution in time, based on the electronics status, trigger settings, and active PMTs. The outputs from the MC simulation and the real data are treated in exactly the same way. The energy response and position for the sources placed at the detector centre have been reproduced with a precision better than 0.8% and 1%, respectively [16]. The MC simulations are needed  for every Borexino analysis and are especially relevant for the evaluation of systematic uncertainties (see  Sections 3 and 4.) The overall agreement of data and MC in the energy region below 3 MeV is within the order of 1% while above 3 MeV it is of the order of 1.9% [13].

Background levels
In Borexino, the backgrounds can be classified into internal, external, and cosmogenic [9,1,17]. The internal background isotopes, namely, 238 U and 232 Th chains, 14 C, 85 Kr, 210 Pb, 210 Bi, 210 Po, and 208 Tl are contaminants of the liquid scintillator itself. The external background components originate from materials outside of the LS and are typically represented by γs' that are able to reach the LS volume. The cosmogenic background consists of muons and consecutive events created by muon spallation in the detector and the surrounding rock. The backgrounds are summarised in Tables 1 and 2 and are discussed below.

Internal background
238 U chain (τ = 6.4×10 9 years): a primordial long-lived radioactive isotope with 99.3% abundance in natural Uranium. The chain contains eight α and six β decays. The chain contains the fast 214 Bi(β − , Q = 3.272 MeV)-214 Po(α, Q = 7.686 MeV) decay sequence with τ = 238 µs, allowing a coincidence tagging assuming secular equilibrium. This chain is highly suppressed by the water extraction campaign as a consequence of the LS purification, leading to an upper limit of < 9.5 × 10 −20 g/g (95% C.L.) on the whole chain. 232 Th chain (τ = 20.2×10 9 years): a primordial long-lived radioactive isotope with 100% abundance in natural Thorium. The chain has six α and four β decays. Assuming secular equilibrium, it is possible to determine its content through the fast 212 Bi(β − , Q = 2.252 MeV)-212 Po(α, Q = 8.955 MeV) decay sequence with τ = 433 ns. The content of this chain after the LS purification campaign reached an upper limit of < 5.7 × 10 −19 g/g (95% C.L.). From this chain, 208 Tl (Q = 4.999 MeV, τ = 4.4 min), which simultaneously emits an electron and gamma rays, has high relevance for the 8 B analysis (see Section 3.2). Since there is also 208 Tl background originating from the IV contamination, the component coming from the 232 Th contamination of the LS is specifically called Bulk 208 Tl. 14 C (β − decay, Q = 0.156 MeV, τ = 8270 years): this isotope is a natural component of the organic liquid scintillator and is chemically identical to the stable isotope 12 C. Therefore, it cannot be removed through purification. It dominates the low-energies, relevant for the pp-ν analysis, and dictates the trigger rate which is about 25 Hz at 50 keV threshold. Its rate is stable in time and determined as (40 ± 2) Bq/100 ton. In Borexino, the relative ratio of 14 C to 12 C is ≈ 10 −18 g/g [18]. 85 Kr (β − decay, Q = 0.687 MeV, τ = 15.4 years): this isotope occurs in the air due to nuclear explosions. Its decay rate can be determined by following the procedure described in [9], which exploits the 85 Kr-85m Rb fast delayed coincidence decay with a branching ratio of 0.43%. It is a major background of 7 Be solar neutrinos (Section 3.2) and has been suppressed by a factor of about 4.6 after the LS purification campaign. 210 Pb (β − decay, Q = 0.064 MeV, τ = 32.2 years) and 210 Bi (β − decay, Q = 1.162 MeV, τ = 7.2 days): 210 Pb is an isotope which is contained in the LS. Its very low Q value is below the Borexino analysis threshold. It has a long lifetime, so it can be considered stable in time during few-years long analysis periods. The 210 Pb in the LS is assumed to be in secular equilibrium with 210 Bi, its short-lived daughter nuclei. The content of 210 Bi has been reduced by a factor of about 2.3 after the water extraction period. It is an important background for 7 Be solar neutrinos (Section 3.3) and a major challenge for CNO solar neutrinos (Section 3.4). While 210 Pb is most probably present also on the surface of the IV, there is no evidence that it would be leaching to the inside of the scintillator, causing an additional source of 210 Bi and its daughter, 210 210 Po events, that can be selected on an event-by-event basis using the MLP variable, are an important "standard candle" to follow the stability of the detector response over time.
External background External background is represented by the particles created outside of the scintillator but reaching the fiducial volume of the analysis. Typically, only γ-rays represent external background and other particles, as αs' and βs' produced in external materials, are absorbed before being able to enter the central parts of the LS. The contamination levels of the detector's construction materials, such as PMTs, SSS, or IV, are extensively discussed in [21]. The external background can be divided into three categories: (1) γs' from 40 [17]. The exception is the 11 C rate that is taken from the Phase-II results on solar neutrinos, measured in 71.3 ton [1]. The n-capture time in Borexino was measured using the 241 Am- 9 Be neutron calibration source [20].
analysis [13], which uses peripheral areas of the IV for detection. When 208 Tl decays, it can be located within the nylon membrane (surface 208 Tl) or in the fluid in a close proximity to the IV (emanated 208 Tl). The latter component can be caused, for example, by 220 Rn, a volatile progenitor that can diffuse out of the IV.
[(α, n), X] → γ: this background is represented by high energy gammas produced in captures of radiogenic neutrons. The latter are produced by (α, n) interactions triggered by αs' from decays of 238 U/ 235 U and 232 Th chains occurring in the SSS and the PMTs' glass. The MC simulation shows that these neutrons are mainly captured on the SSS iron and on the protons and 12 C in the 80 cm buffer layer adjacent to the SSS. The energy of the gammas from these neutron captures extends upto 10 MeV. This background is considered in the energy ranges from 3.2 to 5.7 MeV and from 5.7 to 16.0 MeV, targeted for the 8 B neutrino analysis (see Section 3.2).
Cosmogenic background Cosmogenic backgrounds in Borexino can be divided into three main categories: cosmic muons, cosmogenic neutrons, and cosmogenic radioisotopes [17]. Cosmic muons are created due to the interaction of high energy primary cosmic rays with the nuclei in the atmosphere. Cosmogenic fast neutrons can arise from the spallation of muons passing either through the OD and/or the ID, or the surrounding rocks and can penetrate through the detector materials, due to their high energy and no charge. The cosmogenic radioisotopes are created due to the spallation of cosmic muons on detector materials. In contrast to cosmogenic neutrons, the charged ions of radioactive isotopes have low penetration ability and act as backgrounds only when produced inside the liquid scintillator. The three categories of cosmogenic backgrounds in Borexino are discussed below. The detection and measurement of cosmogenic background are extensively discussed in [20,17,7].
Cosmic muons: The primary cosmic muon flux arriving at the Earth's surface is about 6.5×10 5 m −2 h −1 and is attenuated by a factor of ∼10 6 at LNGS due to the mountains above and this corresponds to a measured flux of (3.432±0.003) · 10 −4 m −2 s −1 [22]. The mean energy of muons at LNGS is about 280 GeV, compared to about 1 GeV at the Earth's surface, since the lower energy muons incident at the surface are absorbed, and the flux steeply falls as a function of, abou energy. Muons in Borexino are classified into three types: internal, external and special muons, which are extensively discussed in [7]. The internal muons, about 4300 per day, are the ones crossing both the OD and the ID. They are identified using three flags: either by a special electronics flag called the Muon Trigger Flag, which means that the Water Cherenkov OD triggered; or through the software reconstruction algorithm called the Muon Cluster Flag, that identifies clusters of hits among those detected by the OD; or via the Inner Detector Flag, that uses different cluster variables for the reconstruction of muon pulse-shape information in the ID. The dead time applied after these muons differ for different analyses, depending on their relevant cosmogenic backgrounds. External muons cross only the OD and do not form clusters in the ID. They are detected by either the Muon Trigger Flag or the Muon Cluster Flag with an overall rate similar to that of internal muons. A 2 ms dead time, nearly 8 times the neutron capture time, is applied after all external muons to eliminate fast neutrons crossing the LS after these muons. Special muon flags are designed to tag a very small category of muons that cross the detector, typically in times when the detector was in a special state. These special categories of muons include also noise events and are particularly important in the geoneutrino analysis (see Section 4.2), where the signal rate is extremely low. Therefore, depending on the analysis, special muons can also be treated as internal muons [7]. In addition to the different muon tagging methods mentioned above, the FADC-sub system allows for an accurate pulse shape discrimination of muons. It plays a key role in analyses where the muon tagging is extremely important such as the geoneutrino analysis [7] and the hep solar neutrino analysis [13]. The combined muon tagging efficiency of the main DAQ and the FADC sub-system in Borexino is 99.9969% [7]. The measurement of muons using 10 years of Borexino data, and their seasonal and annual modulations are discussed in detail in [22].
Cosmogenic neutrons: Cosmic muons in Borexino can lead to the creation of cosmogenic neutrons, due to different spallation processes on carbon nuclei. Neutrons in the Borexino LS are captured with a lifetime of (254.5 ± 1.8) µs (measured using 241 Am- 9 Be neutron calibration source [20]) and create a 2.2 MeV γ when captured on a proton, or a 4.95 MeV γ when captured on a 12 C nucleus. The 2.2 MeV γ is not relevant for the solar neutrino analysis, but the 4.95 MeV γ is important for the 8 B solar neutrino analysis [13]. A 1.28 ms gate is opened after each internal muon to guarantee high detection efficiency of cosmogenic neutrons. The above-mentioned 2 ms dead time (∼8 times the n-capture time) is usually enough to suppress these fast neutrons arising from the passage of muons. Fast neutrons from muons crossing the water tank and from undetected muons crossing the surrounding rocks are relevant backgrounds for the geoneutrino analysis. The neutrons from the water tank are estimated through the analysis of the acquired data, while a dedicated Monte-Carlo simulation is required to estimate the contribution from the surrounding rocks. This is discussed in detail in [7].
Cosmogenic radioisotopes: The spallation of cosmic muons on 12 C nuclei leads to the creation of 11 C isotope (β + decay, Q = 0.960 MeV, τ = 29.4 min), an important background for pep solar neutrinos. Due to its long lifetime of 29.4 min, 11 C cannot be suppressed with a simple time veto. It is tagged through the three-fold coincidence (TFC) algorithm, which exploits the fact that 11 C is mostly produced in time coincidence with neutrons, and further divides the data into TFC-enriched and TFC-depleted spectra for the solar neutrino analysis, which is discussed in more detail in Section 3.2. In addition, since the decay mode of this isotope is via β + , it can be distinguished, on a statistical level, from the solar neutrino signal (β − ) through pulse-shape discrimination, as discussed in Section 2.2. The other cosmogenic isotopes relevant for the solar and geoneutrino analyses are listed in detail in Table 2. All these backgrounds are relevant for the 8 B solar neutrino analysis. They are suppressed using a long dead time of 6.5 s after a muon signal, with the exception of 10 C and 11 Be. 10 C requires a longer time veto and an additional space veto, while the 11 Be background is treated using a multivariate fit and its residual rate is found to be compatible with zero [13]. 9 Li represents an important non-antineutrino background for geoneutrinos due to its (β − + n) decay mode, which imitates the geoneutrino signal (Section 2.5). They are suppressed using sophisticated time and spatial vetoes [7]. Other isotopes such as 8 He and 12 B are are also relevant for the geoneutrino analysis and are discussed in [7], but their contribution is negligible when compared to the 9 Li isotope.

Neutrino and antineutrino detection
The detection principles of neutrinos and antineutrinos are significantly different in Borexino and are discussed in this Section.
Neutrino detection Neutrinos ν of all flavours are detected via the neutrino-electron elastic scattering process: ν e,µ,τ + e → ν e,µ,τ + e, in which neutrinos interact with the electrons present in the LS, which has a density of (3.307 ± 0.015) × 10 31 electrons per 100 ton [3]. In this process, a fraction of the neutrino energy is transferred to the electron, which is finally responsible for the generation of scintillation light in the detector. The electron recoil spectrum, continuous even in the case of mono-energetic neutrinos, extends up to a maximum energy T max e given by: where m e is the electron mass and E ν the neutrino energy. The elastic scattering process has no threshold. The cross section for ν e s' is in the order of 10 −45 to 10 −43 cm 2 for solar neutrino energies, i.e. below 20 MeV [23]. It is about 5 times larger with respect to the ν µ,τ − e scattering process. This is due to the fact that the latter interact only through neutral current (NC) interactions, while ν e s' additionally interact via charged current (CC) interactions. In Borexino, however, the scattering process induced by ν µ s' and ν τ s' cannot be distinguished from each other with the current amount of data. The cross sections used for the solar neutrino analysis consider leading order radiative corrections and are taken from [23], with improved measurements taken from [24].

Antineutrino detection
Electron antineutrinosν e are detected via the Inverse Beta Decay (IBD) reaction: in which an electron flavour antineutrino is captured on a free proton (Hydrogen nucleus), producing a neutron and a positron. Hereby, the LS has a density of (6.007 ± 0.001) × 10 30 protons per 100 ton [7]. The IBD has a kinematic threshold of 1.806 MeV due to the larger mass of the neutron compared to the proton. The positron first deposits its kinetic energy and then annihilates, producing two gammas with E γ = 0.511 MeV each. These two processes cannot be distinguished and lead to the creation of a prompt event. The energy of the antineutrino is mostly transferred to the positron and thus the visible energy of a prompt event E e + vis can be directly connected with the energy of the incident antineutrino Eν e : After its creation, the neutron is thermalised and then it is captured after a typical time of (254.5 ± 1.8) µs [20]. The capture is accompanied by a de-excitation gamma. In the majority of cases, the neutron is captured on a proton and the energy of the gamma is 2.2 MeV. With a probability of 1.1% [7], a 4.95 MeV γ is emitted after the neutron capture on 12 C. Each gamma of this energy range deposits its energy in the scintillator predominantly by multiple Compton scatterings. Several Compton electrons are then detected as a single delayed event. At MeV energies, the IBD cross section is in the order of 10 −42 cm 2 [25], which is about 100 times larger compared to neutrino-electron elastic scattering. Antineutrinos are able to interact also via elastic scattering, but it is much more convenient to detect them via the IBD process, providing a golden channel to identify the rare interactions and significantly suppress backgrounds, exploiting the fast prompt-delayed coincidence signal.

Solar neutrinos
The Sun is a strong natural source of neutrinos, and the emitted flux of solar neutrinos is of the order of 10 10 cm −2 s −1 , with their energy spectrum extending up to about 15 MeV. They are produced in the electron flavor (ν e ) along the nuclear fusion processes that occur in the core of the Sun. Differently from photons, also produced in these interactions, neutrinos are able to travel directly from the production site to the Earth, without being deflected or absorbed. Therefore, solar neutrinos are a direct probe to the Sun's interior. Indeed, they are being extensively used to understand the fundamental processes powering our star since decades [26-30, 9, 31]. Historically, solar neutrino measurements led to the experimental evidence of neutrino flavor transformation [32,33]. Even today, they are at the base of the most precise determination of the θ 12 mixing angle [34]. More recently, they are being used in searches for physics beyond the Standard Model [35,5] and are among the goals of future experiments [36][37][38][39]. The vast majority (about 99%) of the energy produced in the Sun comes from a series of reactions fusing Hydrogen to Helium, called pp chain. The associated neutrino flux is generated by various sub-processes, and includes the so-called pp, pep, 7 Be, 8 B, and hep neutrinos. The remaining small fraction of solar energy is produced in the so-called CNO cycle, in which the Hydrogen-to-Helium fusion is catalysed by the presence of Carbon, Nitrogen, and Oxygen. More details about the production mechanism and propagation of solar neutrinos from the Sun to the Earth are given in Section 3.1. In particular, Section 3.1.1 reports about the pp chain and CNO cycle. Section 3.1.2 discusses the so-called Standard Solar Model (SSM) that predicts the fluxes of different species of neutrinos, that depend on the so-called metallicity, i.e., the abundance of elements heavier than Helium. Solar neutrinos arrive on the Earth as a mixture of all flavours. The process of the flavour conversion of solar neutrinos, maximised for neutrinos with energy greater than ∼2 MeV by the presence of the dense solar matter via the Mikheyev-Smirnov-Wolfenstein (MSW) effect [40,41] is briefly discussed in Section 3.2.1. Borexino is the only experiment that has performed a complete spectroscopy of solar neutrinos. Section 3.2 presents the basic principles of Borexino solar neutrino analysis and underlines the features common to various analysis aimed to extract rates of different solar neutrino species. The following sections then discuss the particularities of different analyses, their results, and physics implications. Section 3.3 describes specifically the measurement of pp chain [1], while Section 3.4 is devoted to the discovery of CNO neutrinos [3]. Finally, Section 3.5 gives a brief overview of Beyond-the-Standard-Model physics probed by Borexino. Searches for spectral deformation of electrons scattered off 7 Be solar neutrinos due to eventual Non-Standard Interactions (NSI) are described in Section 3.5.1, while the tight limits set on Neutrino Magnetic Moment (NMM) are discussed in Section 3.5.2.

Hydrogen to Helium fusion: pp chain and CNO cycle
Solar neutrinos are emitted during the fusion of protons to Helium nuclei taking place in the solar core: The dominant fusion process is the pp chain, while a subdominant fraction of solar energy is produced in the so-called CNO cycle, in which the fusion is catalysed by the presence of heavier elements, namely    Table 3 after constraining the solar luminosity as presented in [43].
Carbon, Oxygen, and Nitrogen. Thus, the CNO contribution to the Sun's fusion processes depends directly on the core's metallicity, i.e. the mass abundance of elements heavier than Helium. In addition, the net contribution of the CNO cycle is strongly dependent on the core's temperature. The Sun is a relatively small and cold star in the Universe and the CNO contribution to the luminosity budget is around 1%. For heavier stars, having a mass greater than ∼1.3 solar masses, the CNO cycle is instead believed to be the dominant process which burns Hydrogen into Helium. It is therefore considered as the main nuclear fusion process occurring in the Universe. The precise CNO contribution in the Sun, however, is unknown, since the Sun's metallicity is not known with precision (see Section 3.1.2). The top part of Figure 4 shows the schemes of the pp chain and of the CNO cycle, while its bottom part shows the energy spectrum of solar neutrinos. The fluxes, and thus the normalisation of the energy spectra, are predicted by the SSM [42] as discussed in the following Section. The pp chain consists of three branches, indicated as pp-I, pp-II, and pp-III, each terminated by the production of 4 He. The flux of solar neutrinos is dominated by the pp neutrinos (order of 10 10 s −1 cm −2 ) having a continuous energy spectrum with 0.420 MeV endpoint. In the pp chain, also mono-energetic 7 Be (10% branching at 0.384 MeV (excited state) and 90% branching at 0.862 MeV (ground state)) and pep (1.44 MeV) neutrinos are produced, as well  Table 3: Solar neutrino fluxes predicted by Standard Solar Models B16-GS98 (High Metallicity, HZ-SSM) and B16-AGSS09met (Low Metallicity, LZ-SSM) [42] in units of cm −2 s −1 with the exponential factors given in the last column. The relative difference between the two model predictions is also quoted. In the Borexino analysis, we call the CNO spectrum the weighted sum of all three components.

Standard Solar Model and the metallicity problem
The so-called Standard Solar Model (SSM) is a spherically symmetric quasi-static model of a star in hydrostatic equilibrium with one solar mass M , including several differential equations derived from basic physical principles. SSM assumes that the Sun was initially chemically homogeneous and that the mass loss is negligible during the whole 4.57 Gyr of its existence [42]. The calibration is done by adjusting the mixing length parameter and the initial Helium and metal 2 mass fractions in order to satisfy the constraints imposed by the present-day solar luminosity L , radius R , and surface metal-to-hydrogen abundance ratio (so-called metallicity, (Z/X) ) [42]. SSM assumes that solar energy is generated by the pp chain and the CNO cycle.
Among the outputs of the SSM are the neutrino fluxes, summarised in Table 3 for the new generation of SSM called B16 [42], that includes updates on nuclear reaction rates, more consistent treatment of the equation of state, and a novel treatment of opacity uncertainties. The prediction is given separately for the two canonical sets of solar abundances. So-called low metallicity (LZ or AGSS09met-LZ) scenario [44,45] represents the most recent revision of solar abundances based on development of three-dimensional hydrodynamical models of the solar atmosphere, of techniques to study line formation, and the improvements of atomic properties such as transition strengths. The other solar metal composition scenario, so-called high metallicity (HZ or GS98-HZ) scenario [46], is based on older one-dimensional modeling of the solar atmosphere and predicts higher metal abundances. In this regard, emerges the so-called metallicity puzzle. The fact is that newer LZ-SSM spoil the earlier agreement between the observed sounds speed profile (helioseismology data) and the corresponding SSM predictions. The origin of this discrepancy is currently not understood. However, the SSM prediction of neutrino fluxes depend on the metallicity: the metallicity influences the opacity of the Sun and consequently also the temperature in the core and the fusion rates. There is a sizeable difference of 8.9% and 17.6% between the HZ and LZ SSM predictions of 7 Be and 8 B fluxes, respectively ( Table 3). The largest difference between the fluxes predicted by the LZ and HZ SSM results for the CNO cycle and amounts to about 32%. The metallicity indeed directly influences the efficiency of the CNO cycle, since the "metals", Carbon, Nitrogen, and Oxygen are the elements which catalyze the process. The CNO neutrino flux is then directly dependent upon the core's metallicity, which keeps memory of the Sun's elemental composition at the time of formation. Since the metal abundance in the core is believed to not be influenced by the surface, CNO neutrinos preserve the core's information in its pristine conditions. Thus, neutrinos produced in the CNO cycle are a unique probe to the Sun's primordial composition. In summary, precise measurements of the solar neutrino fluxes can provide important boundary conditions for the future development of the SSMs and our understanding of the stars in general.

Neutrino flavour conversion and matter effects
In the standard three-flavour neutrino framework, the electron neutrino survival probability, i.e., the probability to measure solar neutrino in the same flavour as it was produced, can be cast in the form [47]: where P 2ν (ν e → ν e ) corresponds to the 2ν probability for θ 13 = 0, which depends on the solar mixing angle θ 12 and the mass splitting ∆m 2 12 only. Should the oscillation happen in vacuum, this survival probability could be approximated by [1]: that does not depend on energy and has an approximate value of 0.54. In reality, solar neutrinos are crossing the dense solar matter. The electron flavour neutrinos experience an extra potential due to the charge-current interaction with the electrons present in the Sun. This affects the neutrino oscillation probability that changes with respect to a pure vacuum oscillation scenario. This effect is called Mikheyev-Smirnov-Wolfenstein (MSW) effect [40,41]. Thus, the survival probability P MSW ee depends not only on the oscillation parameters, but also on the neutrino-energy-dependent potential. Assuming an adiabatic decrease of the electron density with radius, it can be expressed as follows [48,1]: (ν e → ν e ) = P MSW ee 0.5 cos 4 θ 13 1 + cos 2θ M 12 cos 2θ 12 , where cos 2θ M 12 = cos 2θ 12 − β (cos 2θ 12 − β) 2 + sin 2 2θ 12 (11) and where θ M is the mixing angle in matter, E ν the neutrino energy, n e the electron density in matter, and G F the Fermi coupling constant. This MSW survival probability stays at the vacuum value at low energies typical for pp neutrinos (vacuum dominated region), while for high energy 8 B neutrinos it decreases to about 0.32 (matter dominated region). The situation is further complicated by the fact, that different neutrino species have their production regions at different radii [49], and thus, propagate through regions of different electron densities. Calculation of the survival probabilities considering non-adiabatic corrections and  Table 4: Details of the solar neutrino analysis of the Borexino Phase II(+I) and Phase III data. The first column specifies the solar neutrino species extracted by the fit of the data using the exposure mentioned in the second column. The third column shows the variables used in the fit: global variable E, standing for various energy estimators (Section 2.2) or R, reconstructed radius. The parameters in square brackets identify the variables used in the multivariate part of the fit in the LER: pulse shape variable used to constrain residual 11 C(e + ) and radial distribution constraining external background. The fourth column gives the energy range of the global variable, while the last column highlights the species whose rates were constrained in the fit. The source of information used to extract these independent constraints is given in brackets, while more details are given in text. UL = Upper Limit; Cosm. = cosmogenic background; indep. data = independent data; LER/HER = low/high energy region; LPoF = Low Polonium Field.
averaging over the production region for each solar neutrino species has been presented in [50]. Finally, the exact form of the transition region between the vacuum and matter dominated regions might be sensitive to different models of the non-standard neutrino interactions and thus, is a point of interest for searches for physics beyond the Standard Model [47]. Solar neutrinos are detected via the elastic scattering of electrons, as it was discussed in Section 2.5. Therefore, even for mono-energetic neutrinos, as 7 Be or pep neutrinos, the spectrum of scattered electrons is a continuous one, characterised only by the Compton-like edge, corresponding to a maximal energy of the scattered electrons (Equation 4). It is impossible to distinguish the electrons scattered off by solar neutrinos from the background components, on an event-by-event basis 3 . Therefore, the analysis proceeds in two steps: (1) the event selection, with a set of cuts to maximize the signal-to-background ratio, and (2) the extraction of the neutrino and residual background rates through a combined fit of the distributions of global quantities of the events surviving the cuts. Table 4 summarizes some basic details of the solar neutrino analysis discussed in this paper such as the extracted solar species, analyzed periods, exposure, fit variable, energy range, and constraints used in the fit. Measurement of the pp chain neutrinos is discussed in Section 3.3, based on [1,10,13]. The interaction rates of the pp, 7 Be, and pep neutrinos were obtained through a spectral fit of the Phase-II data in the so-called LER (Low Energy Region) below 3 MeV. The measurement of 8 B neutrinos was performed on the combined Phase-I+II data, through a fit of the radial distribution. This strategy avoids any assumption on the spectral shape, i.e., on the shape of the survival probability in the transition region, defining the so called "upturn", the energy interval sensitive to possible NSI (Section 3.2.1). This analysis is performed in the HER (High Energy Region) above 3.2 MeV and below 16 MeV. This energy interval is divided into two sub-parts namely, HER-I (below 5.7 MeV) and HER-II (above 5.7 MeV), each characterised by different backgrounds. The experimental confirmation of the existence of the CNO fusion in the Sun [3] was performed using the Phase-III data in the energy interval similar to LER, with an increased energy threshold, due to the worsened resolution associated with the loss of PMTs. The CNO analysis is discussed in more detail in Section 3.4.
The main event selection criteria are conceptually similar for all solar neutrino analyses and are conceived to reject cosmic muons surviving the mountain shield, reduce the cosmogenic background, select an optimal spatial region of the scintillator (the fiducial volume, FV), and remove eventual noise events.
The muon detection combines the information of the external Cherenkov veto with that of the inner detector, including the pulse shape analysis. Cosmogenic background is reduced by applying a veto following each muon. A 2 ms veto suppresses neutron captures from external muons, which cross only the water buffer. A vast majority of captures occur on protons, emitting 2.2 MeVγs', while around 1% of the captures occur on 12 C nuclei, emitting 4.95 MeV γs'. For the internal muons crossing the scintillator, the veto time differs. In the LER, a time cut of 300 ms is used to suppress the saturation effects of electronics. In the HER, a much longer veto of 6.5 s is applied, in order to suppress 12 B, 8 He, 9 C, 9 Li, 8 B, 6 He, and 8 Li decays. In addition, the presence of untagged 11 Be (Q = 11.5 MeV, β − , τ = 19.9 s) is estimated through a multivariate approach and found to be compatible with zero. In HER, an additional spherical cut of 0.8 m radius is also applied for 120 s around the capture position of cosmogenic neutrons to remove 10 C(Q = 3.6 MeV, β + , τ = 27.8 s) background.
Specifically in the LER analysis, the cosmogenic 11 C(Q = 0.96 MeV, β + , τ = 29.4 min) background requires a dedicated treatment. Due to its relatively long life-time, it cannot be removed by a simple veto. The so-called three-fold coincidence (TFC) algorithm [9,10] uses the fact that 11 C, created by the spallation of muons on 12 C, is mostly produced together with neutrons: Based on the characteristics and mutual configuration of the muon and cosmogenic neutrons, the TFC algorithm identifies the space-time regions with increased probability to create 11 C and/or assigns to each event a probability to be 11 C. Based on this, all events are divided into the TFC-subtracted and TFC-enriched categories. The TFC-subtracted spectrum preserves about 64% of the total exposure, while the 11 C suppression preserves more than 90%. In both, LER and HER, 214 Bi and 214 Po events are removed with about 90% efficiency, through the space-time correlation of their fast β + α delayed coincidence.
The low-and high-energy analyses use different fiducial volumes (FV). In the LER, the FV represents the central region of 71.3 ton, selected to maximally suppress external γs' from 40 K, 214 Bi, and 208 Tl, originating from the materials surrounding the scintillator. It is contained within the radius R < 2.8 m and the vertical coordinate -1.8 m < z < 2.2 m. The HER is above the energy of the aforementioned external background. The analysis in HER-I requires only a z < 2.5 m cut to suppress background events related to a small pinhole in the inner vessel that causes liquid scintillator to leak into the buffer region. The total selected mass in this case is 227.8 ton. In contrast, the analysis in HER-II uses the entire scintillator volume of 266 ton, since the above-mentioned external background does not affect this energy window.
Backgrounds which survive the selection cuts are treated in a Poissonian binned likelihood fit to disentangle them from solar neutrino signals. In the HER-I and HER-II, a fit of the radial distribution of events is performed to separate the 8 B neutrino signal (uniformly distributed in the scintillator) from the external background. The LER analysis follows a multivariate approach, which simultaneously fits the two (TFC-subtracted and TFC-tagged) energy spectra, the radial, and in Phase-II, also the pulse-shape estimator distributions. The radial fit helps to constrain the external background. The pulse shape distribution in Phase-II is used to constrain residual 11 C(e + ) in the TFC-subtracted spectrum. The likelihood function is constructed by the multiplication of the likelihoods corresponding to the listed distributions: where the last term was used only in Phase-II. The free parameters of the fit are the rates of solar neutrinos and background components from zero threshold. When needed, some of these can be constrained by multiplicative pull terms (Gaussian or semi-Gaussian) to break the correlation among species having similar spectral shapes. This is, of course, only possible, when there is an independent way to evaluate these rates. This will be discussed in the following Section 3.3 for Phase II and Section 3.4 for Phase-III analyses. A summary of these conditions is shown in the last column of Table 4. The probability distribution functions (PDFs) used in the fit for signal and backgrounds are typically obtained by the Geant-4 based MC simulation (Section 2.3) 4 . This MC-based method has the advantage that the detector response is automatically taken into account by the simulation. The main disadvantage is that an extensive analysis campaign is needed in order to evaluate the systematic uncertainties related to the imprecision of the MC-based PDFs. In addition, this approach cannot adjust for eventual changes that might appear in the detector. An alternative is the so-called analytical fit of the energy spectra used in the Phase-II analysis [10]. Here, the detector response is represented by analytical functions, with parameters such as the effective light yield, non-linearity of the energy scale, resolution and non-uniformity parameters. In the fit, some of these parameters are kept free and some are constrained or fixed. Thus, this approach is more flexible to adjust for changes in the detector and does not require to evaluate systematic uncertainty towards the parameters left free in the fit. The disadvantage of this approach is the usage of several fit parameters, which makes this approach generally more prone to correlations.

Spectroscopy of pp-chain neutrinos
This section is dedicated to the latest Borexino analysis of pp chain solar neutrinos [1,13,10]. After the discussion of the main principles of this analysis in the previous section, the particularities of the analysis strategy of pp chain solar neutrinos is presented in Section 3.3.1 and the results in Section 3.3.2. The discussion of their implications is presented in Section 3.3.3. Observation of the seasonal modulation of 7 Be neutrino rate -a direct evidence of its solar origin, is briefly presented in Section 3.3.4.

Analysis strategy
As mentioned in the previous section, the LER and HER follow different analysis strategies. These are explained below.

LER Analysis
The goal of this analysis is to measure pp, 7 Be, and pep neutrino interaction rates with the same multivariate fit in the energy interval from 0.19 to 2.93 MeV. Due to the different energy interval of these three neutrino species, the extraction of each species faces different challenges. In order to make it easier to follow this section, we exemplify the spectral fit in Figure 5.
For low energy pp neutrinos, there is a high correlation with the irreducible 14 C background. Nevertheless, the 14 C rate can be independently determined from the fit of the energy spectrum of the second clusters [51], i.e., events randomly falling in the last part of the 16 µs long data acquisition window triggered by the preceding event (first cluster). The obtained 14 C rate is (40 ± 2) Bq/100 ton and this value was used to constrain it in the fit.
The other challenge for pp solar neutrinos are pile-up events, i.e., events happening so close in time to each other, that they are reconstructed as single events. In this case, events occurring far away from each other in space, even out of the FV, can be erroneously reconstructed in the FV. Pile-up is dominated by the overlap of 14 C+ 14 C, but non-negligible contributions arise also from the pileup between external background with either 14 C or 210 Po. Pile-up is treated using the following two methods described in [51,16]. In one case, synthetic pile-up spectrum is constructed, starting from real data or MC, that is used as an additional spectral component, fully constrained in shape and rate during the fit. In the other case, all spectral components are convoluted with a randomly acquired spectrum, i.e. with events acquired with a solicited, external trigger. Thus, in this approach, all energy PDFs are slightly deformed (with the dominant change only on 14 C spectrum) and no additional component is added in the spectral fit.
The Compton-like edge of the 7 Be neutrino signal is a clearly visible feature in the spectrum, see Figure 5. This spectral feature makes the fit relatively easy, even if there is a correlation with 85 Kr and 210 Bi backgrounds.
The measurement of pep neutrinos is complicated by the presence of 11 C background, that is treated by the TFC technique discussed in Section 3.2. It is the measurement of pep neutrinos that required the multivariate fit approach. In addition, there is a strong correlation between pep, 210 Bi, and CNO spectral shapes. Thus, in order to break this degeneracy, the CNO rate is constrained in the fit to the SSM prediction including MSW-LMA oscillations. The analysis is repeated for both HZ-SSM and LZ-SSM, with the expected CNO rate of (4.92 ± 0.55) cpd/100 ton and (3.52 ± 0.37) cpd/100 ton, respectively. In case of different results, these are quoted separately.
The 8 B neutrino rate does not affect the LER analysis and is always fixed to its SSM prediction, while the hep neutrino rate is simply neglected due to its fully negligible expected rate.

HER Analysis
The strategy to extract the 8 B solar neutrino rate is based on radial fits in the HER-I and HER-II, that are shown in Figure 6. In both fits, the dominant uniform contribution is from 8 B solar neutrinos. Some residual uniform background due to muons, cosmogenic isotopes, and 214 Bi decays surviving the cuts is very small. This is estimated following the procedure in [19], and constrained in the fit. Only in HER-I there is an additional uniform background from bulk 208 Tl(Q = 5 MeV, β − , γ), which comes from the residual 232 Th contamination of the liquid scintillator and is constrained in the fit to the value based on the 212 Bi-212 Po (β + α) fast delayed coincidences. External 208 Tl contamination contributes to the HER-I with two distinct components: one from contamination directly on the inner vessel surface, and another from decays of nuclei that have recoiled from the inner vessel into the liquid scintillator or originated from the volatile progenitor of 208 Tl, 220 Rn, which has emanated from the nylon. The rates of both components are left free to vary in the radial fit. Finally, HER-I and HER-II are also polluted by γ-rays following the capture of radiogenic neutrons produced via (α, n) or spontaneous fission reactions of 238 U, 235 U, and 232 Th in the Stainless Steel Sphere (SSS) and PMTs. This rate is also a free parameter in the fit.   and Phase-III [3] solar neutrino analyses. The rates and fluxes are integral values without any threshold; the first error is statistical, the second systematic. The rate-to-flux conversion assumes neutrino flavour conversion [50] with the neutrino oscillation parameters from [47]. The last column shows the fluxes as predicted by the HZ-and LZ-SSM [42] (see Table 3). The result for 7 Be neutrinos contains the ground and excited state lines as shown in Figure 4. The fluxes of pp, 7 Be, pep, CNO, 8 B, and hep neutrinos are normalised to 10 10 , 10 9 , 10 8 , 10 8 , 10 6 , and 10 3 (see Exp factor definitions in Table 3), respectively.
The neutron captures are the only background in the HER-II analysis, as there are no naturally long-lived isotopes above 5 MeV.

Results on pp, pep, 7 Be, and 8 B neutrinos
The analysis strategies for the LER and HER were explained in the previous Section. The results from these analyses are discussed below.   [10]. Further details are discussed in the text.

LER Analysis
The result of the multivariate MC fit in the LER is shown in Figure 5, which illustrates the four fits, namely the TFC-subtracted, TFC-tagged, pulse shape, and radial distributions, as defined in Equation 14. The fit results in terms of interaction rates of solar neutrinos in counts per day per 100 ton (cpd/100 ton) are given in Table 5 including systematic errors. The fact that the fit is repeated with the CNO rate constrained to HZ-and LZ-SSM predictions influences only the resulting pep neutrino rate and is thus given separately with the label HZ and LZ. In both cases, the absence of the pep reaction in the Sun is rejected with > 5σ significance, which makes this measurement the discovery of solar pep neutrinos.
In spite of the remarkable understanding of the detector response throughout the scintillator volume and in a large energy range (Section 2.3), an extensive study of the possible sources of systematic errors has been performed. The results of these studies are summarised in Table 6, which lists the various contributions to the systematic error individually for the pp, 7 Be, and pep measurements. The main contribution to the systematic error comes from the fit model, that is, possible residual inaccuracies in the modelling of the detector response (energy scale, uniformity of the energy response, pulse-shape discrimination shape) and uncertainties in the theoretical energy spectra used in the fit. The second source of systematics is related to the fit method, i.e. eventual differences between the MC-based and analytical fit approach. Further systematic effects arise from the choice of the energy estimator, the details of the implementation of the pile-up, different fit energy ranges and binning, the inclusion of an independent constraint on 85 Kr, and the estimation of the FV. This last uncertainty is determined with calibration data, using sources deployed in known positions throughout the detector volume.

HER Analysis
The radial fits in HER-I and HER-I ranges to obtain the 8 B neutrino interaction rate are shown Figure 6. The resulting 8 B interaction rates are given in Table 5. The most important systematic uncertainties arise from the determination of the target mass, that is complicated by the presence of the small leak in the inner vessel. The evolution of the scintillator mass is monitored on a week-by-week basis, by studying the inner vessel shape, which is obtained from the spatial distribution of its surface contamination. Additional sources of systematic error include the energy scale uncertainty, and the application of the z-cut in HER-I. The uncertainties from the live-time determination and the knowledge of the scintillator density are almost negligible.
To complete the pp chain analysis, a search for hep neutrinos has been performed. Its flux expectation is two orders of magnitudes smaller than that of 8 Tables 4 and 5).

Implications for solar and neutrino physics
The measured interaction rates of solar neutrinos, as discussed in the previous Section, can be used to test our understanding of both the Sun and the basic neutrino properties. Assuming that the physics of neutrino interactions and oscillation are known, the measured rates can be converted to neutrino fluxes, to be then compared individually with the HZ and LZ SSM predictions, which is important for constraining the solar metallicity. In addition, one can quantify the relative intensity of the two primary terminations of the pp chain (pp-I and pp-II in Figure 4) and evaluate the solar neutrino luminosity. On the other hand, assuming the SSM predictions for the neutrino fluxes, one can evaluate the electron neutrino survival probability P ee for different energies and compare them with the standard prediction of the 3-flavour neutrino oscillations including the MSW effect. The following paragraphs discuss these points in more detail.
pp-chain solar neutrino fluxes Considering solar neutrino oscillations, the expected neutrino interaction rate in Borexino R ν is [9]: where N e is the number of target electrons (see Section 2.5), Φ e is the solar neutrino flux, dλ/dE ν is the differential energy spectrum of solar neutrinos, P ee is the electron neutrino survival probability (see Section 3.2.1 and [50]), and dσ e,µ,τ /dT e are the differential cross sections for the scattering reaction discussed in Section 2.5. The spectrum of solar neutrinos, normalised according to SSM predicted fluxes (last column of Table 3) is shown in the bottom part of Figure 4. The cross-sections of elastic scattering for different neutrino flavours were discussed in Section 2.5. We remind that Borexino has no sensitivity to distinguish between the shapes of recoiled electron spectra from ν e and ν µ,τ . However, for the same neutrino energy, the σ e is about 4-5 times larger than σ µ,τ . Thus, in order to convert the measured interaction rate to flux, it is important to know the relative proportion of the flavours in the measured flux and therefore, P ee . This conversion is relatively simple for mono-energetic neutrinos, such as 7 Be and pep. For solar neutrinos with a continuous energy spectrum and analysis performed in a restricted energy interval of scattered electrons, the situation is more complicated. One has to take into account the energy dependent detector response and assume energy dependence of P ee . This procedure is described in Appendix of [13]. The solar neutrino fluxes converted following this procedure are given in the third column of Table 3. In This conversion assumes that all interacting neutrinos are of electron flavour, which has a higher probability to interact with respect to other neutrino flavours. Thus, in order to comply with the measured rate, a smaller flux is sufficient. This number does not depend on P ee and is therefore very useful to compare the results among different experiments. The Borexino result is compatible with the high-precision result of Super-Kamiokande 2.345 ± 0.014 (stat.)±+0.036 (syst.)×10 6 cm −2 s −1 [31,13].
Metallicity As discussed in Section 3.1.2, the SSM prediction for solar neutrino fluxes depends on the assumption of solar metallicity. For pp chain neutrinos, the difference between the HZ-and LZ-SSM predictions is largest for 7 Be and 8 B neutrinos, 8.9% and 17.6%, respectively, as it is reported in Table 3. When comparing the measured 7 Be and 8 B neutrino fluxes with these theoretical SSM predictions, as shown in Figure 7, it is possible to evaluate the agreement between the data and the HZ-and LZ-SSM predictions. Note that the errors in the Borexino measurements are in both cases smaller than the theoretical uncertainties.
The Borexino results are compatible with the temperature profiles predicted by both HZ-and LZ-SSMs. However, the 7 Be and 8 B solar-neutrino fluxes measured by Borexino provide an interesting hint in favour of the HZ-SSM prediction. A frequentist hypothesis test based on a likelihood-ratio test statistics (HZ versus LZ) was performed by computing the probability distribution functions with a MC approach. Assuming HZ to be true, the LZ is disfavoured at 96.6% C.L. Additionally, a Bayesian analysis gives a Bayes factor of 4.9 showing a mild preference towards HZ-SSM. However, this hint weakens when the Borexino data are combined with data of all the other solar neutrino experiments + KamLAND reactor antineutrino data.
pp-chain branches From the measured solar neutrino fluxes, it is possible to evaluate the ratio R I/I I between the 3 He-4 He and the 3 He-3 He fusion rates, which quantifies the relative intensity of the two primary terminations of the pp chain, a critical probe of the solar fusion. Neglecting the 8 B neutrino contribution, this ratio can be expressed as: The result obtained with Borexino is R I/I I = 0.1780 +0.027 −0.023 . This value is consistent with both the HZ-and LZ-SSM predictions, 0.180 ± 0.011 and 0.161 ± 0.010, respectively.

Solar luminosity and thermal stability
The neutrino fluxes determined experimentally can be used to derive the total power generated by nuclear reactions in the Sun's core [43]. Using the measured Borexino fluxes from Table 3, the obtained luminosity L = (3.89 +0. 35 −0.42 ) × 10 33 erg s −1 is in agreement with the luminosity calculated using the photon output [52,53], L = (3.846 ± 0.015) × 10 33 erg s −1 . This is a robust and direct evidence of the nuclear origin of the solar power. While neutrinos provide a real time picture of the solar core, it takes around 10 5 years for the photons to reach the solar photosphere, from where they are free to escape. The comparison of the two luminosities then also proves that the Sun has been in thermodynamic equilibrium over this timescale.
Electron neutrino survival probability The measured interaction rates of solar neutrinos can be used to extract the electron neutrino survival probability at different energies. This can be done using already discussed Equation 15, assuming standard neutrino interactions and, in this case, SSM fluxes. Figure 8 shows the extracted P ee as a function of the neutrino energy for each measured solar neutrino species. The obtained neutrino survival probabilities are P ee (pp, 0.267 MeV) = 0.57 ± 0.09, P ee ( 7 Be, 0.862 MeV) = 0.53 ± 0.05, P ee (pep, 1.44 MeV) = 0.43 ± 0.11, P ee ( 8 B HER , 8.1 MeV) = 0.37 ± 0.08, P ee ( 8 B HER−I , 7.4 MeV) = 0.39 ± 0.09, and P ee ( 8 B HER−I I , 9.7 MeV) = 0.35 ± 0.09. For continuous neutrino spectra, i.e. for pp and 8 B, the P ee is quoted for the average energy of neutrinos that produce scattered electrons in the given energy range. The quoted errors include the uncertainties on the SSM solar-neutrino flux predictions.
Borexino is the only experiment that can simultaneously test neutrino flavour conversion both in the vacuum and in the matter-dominated regime, providing the most precise measurement of the P ee in the LER. In HER, where the flavour conversion is dominated by matter effects in the Sun, the Borexino results are in agreement with the high-precision measurements of Super-Kamiokande [54,55] and SNO [32,33]. The Vacuum-LMA 5 prediction is shown as a grey band in Figure 8, calculated with Equation 9 using θ 12 and θ 13 values based on measurements of KamLAND [56] and Daya Bay [57], respectively and obtained without Borexino data. The pink band instead shows the MSW-LMA solution [50] with the oscillation parameters indicated in [58]. Borexino data disfavours the vacuum-LMA hypothesis at 98.2% C.L. and are in excellent agreement with the expectations from the MSW-LMA paradigm.

7 Be flux seasonal modulation
The flux of solar neutrinos at Earth is not constant in time, but is instead modulated due to Earth's movement around the Sun on an elliptical trajectory. Given the variation of the solid angle covered by Earth during the year, a variation in the solar neutrino flux is predicted to exist. The net flux variation between the maximum and the minimum is estimated to be around 6.7%.
Borexino was able to measure the annual modulation of solar neutrinos with high significance [9,2], confirming the solar origin of the measured 7 Be signal. In order to maximize the signal-to-background ratio, events were selected in the 7 Be shoulder energy range. In the latest analysis [2], the energy window has been tuned to be [0.215,0.715] MeV and the events have been selected inside a large FV of 98.6 ton. The chosen energy region is also rich in 210 Po α-decays, which are efficiently suppressed by means of the MLP cut (see Section 2.2). In order to extract the modulation signal, three different analytical approaches were used: an analytical fit to event rate, a Lomb-Scargle periodogram, and an Empirical Mode Decomposition analysis. Figure 9 shows the β − event rate in the energy region of interest, in ∼30-day time bins. All methods yield compatible results, confirming the observation of solar neutrino flux modulation. The fit values for the modulation periodicity and its amplitude are well consistent with the expectations. Borexino was able to reject the hypothesis of no modulation with a confidence level of 99.99%.

Detection of CNO neutrinos
This Section reports the details about the analysis performed on Borexino data to extract the CNO signal. Section 3.4.1 describes the challenges of this analysis and the main strategy towards the first observation of the CNO neutrino signal on Phase-III data [3]. The achieved 5σ significance of this results is compatible with the expected sensitivity [4] described in Section 3.4.3. Section 3.4.2 describes the method used to evaluate the rate of the 210 Bi background contaminating the scintillator. This was used as a constraint in the spectral fit as well as in the counting analysis, leading to the final results discussed in Section 3.4.4.

Analysis strategy
The energy of recoiling electrons after CNO neutrino interactions, shows a continuous distribution with an endpoint at 1.517 MeV. The main sensitivity to CNO neutrinos comes from the energy region above the 7 Be shoulder, i.e. from 0.8 to 1.0 MeV [4]. This energy region is populated by other backgrounds, which limit the sensitivity towards the CNO neutrino events. Exactly as in the pp chain LER neutrino analysis (Section 3.3), the cosmogenic 11 C(β + ) background is treated via the TFC algorithm (Section 3.2). In the TFC-subtracted spectrum, the most important backgrounds are the β − particles emitted by the 210 Bi contaminant and the pep-neutrino recoil electrons. These two species have spectral features very similar to that of CNO neutrinos and their absolute rate must then be constrained by means of independent inputs.
In the Phase-II LER analysis, it was not possible to set an independent constraint on 210 Bi. The rate of pep neutrinos was constrained indirectly, by constraining the ratio of pp and pep neutrino fluxes to the SSM predictions of 47.7 ± 0.8 (HZ-SSM) and 47.5 ± 0.8 (LZ-SSM), while the absolute pp and pep rates were free fit parameters. This corresponds to an effective constraint of 10% precision for the pep neutrino rate, a value dominated by the precision with which Borexino can measure pp neutrinos. The pp-pep ratio is known very well from nuclear physics due to the fact that both reactions have the same nuclear matrix element. By performing a ∆χ 2 -profiling, an upper limit on CNO interaction rate of 8.1 cpd/100 ton (95%C.L.) was obtained.
In the Phase-III analysis, the energy threshold was increased above the end-point of pp neutrinos. This was due to a worsened resolution at low energies as a consequence of dying PMTs. The rate of pep neutrinos was then constrained directly to a value of (2.74 ± 0.04) cpd/100 ton, corresponding to a precision of 1.4% [43,4]. This value results from a combination of robust theoretical assumptions and a global fit to the solar neutrino data, excluding Borexino Phase-III. The other main background for the CNO measurement, consisting of 210 Bi decays, was also constrained in the Phase-III analysis. 210 Bi has a short lifetime of 7.2 days, and its overall rate in the detector can be assumed to be in secular equilibrium with its parent nucleus 210 Pb. The decay chain of 210 Pb is summarised in Equation 2. 210 Pb decays are well below the analysis threshold and can therefore be considered invisible. The 210 Po, on the other hand, is an unstable isotope, which produces mono-energetic α particles. α-decays are efficiently detected in Borexino on an event-by-event basis, by means of the MLP selection (see Section 2.2). Provided that the secular equilibrium in Equation 2 is maintained, the measured 210 Po corresponds to the 210 Bi rate.

Low Polonium Field
As discussed previously, the rate of 210 Bi decays can be constrained via its link with the 210 Po decay rate, with the assumption that this latter term is only supported by in-equilibrium 210 Pb decay chain. Data collected by Borexino since its start, however, indicate that an out-of-equilibrium component of 210 Po is present in the detector. The source of this component is likely the surface of the Inner Vessel, from which 210 Po is detached into the scintillator. The mean free path of 210 Po atoms is calculated to be very small in stable conditions. However, the presence of convective motions in the Borexino scintillator allow 210 Po to spread throughout the scintillator volume. Under this conditions, the measured value of 210 Po decay rate would be much higher than the 210 Bi decay rate, spoiling any possible constraint.
To limit convective motions in the scintillator volume, the Borexino Collaboration pursued a longlasting effort, culminated in the detector thermal insulation in 2015 and the subsequent installation of active temperature controls (Section 2.1). This way, 210 Po mixing has been strongly suppressed since 2016, leading to the formation of a very clean region around the centre of the detector, called the Low Polonium Field (LPoF). The in-equilibrium 210 Po decay rate in the LPoF region can be then measured. However, there might be still some residual contribution of convective 210 Po in this region and the measured 210 Po rate can therefore only be translated into an upper limit for the 210 Bi rate.
The 210 Po events in the LPoF have been chosen in the energy range 0.30-0.54 MeV and selected by means of the MLP (see Section 2.2). A paraboloid equation in 2D, assuming rotational symmetry along the x − y plane, has been used to fit the data. The minimum 210 Po rate (R( 210 Po min )) has been then obtained through the fit function: Here, ρ 2 = x 2 + y 2 , z 0 is the minimum position of the LPoF along the z-axis, a and b are shape parameters along the respective axes, ε E and ε MLP are the efficiency of the energy and MLP cuts, respectively, applied to select 210 Po events, and R β is the residual rate of β events after the selection of 210 Po events. Since the LPoF slowly moves along the z-axis due to residual convective motions, data in the LPoF needs to be aligned along the z-direction before performing the fit on the full dataset. This has been done by "blindly" aligning the data in the LPoF every month (or every two months) using the centre z 0 obtained by fitting the data of the previous month. Monthly fits have been performed in large volumes of 70 or 100 ton. After the blind alignment using the centres of every month, the final fit has been performed on the aligned dataset in around 20 ton (∼5000 events) using either a simple paraboloid in Equation 17 with four free parameters (R( 210 Po min ), a, b, z 0 ) or with more free parameters, depending on the method. The simple paraboloid fit can be performed either as a likelihood fit with ROOT [59] or with the MultiNest Bayesian tool [60][61][62].
The assumption of the rotational symmetry has been also verified. In addition to the 2D paraboloid fit, 3D ellipsoidal fits have been also performed with MultiNest without assuming rotational symmetry along the x − y plane, resulting in statistically compatible results. In order to account for the complexity of the LPoF along the z-axis, a cubic spline function was implemented along the z-axis in equation (17) and the fit was performed with MultiNest. Despite its better fit on the LPoF data, the method was statistically compatible with the simple paraboloid fit and both the methods were finally used for the 210 Bi upper limit. In addition to the statistical uncertainty of the fit, the other sources of uncertainties considered in this analysis are: • Systematic uncertainties from the fit: mass of the fit region, and binning of the data histogram.
• Homogeneity of β-events: Since the 210 Bi upper limit is estimated from the LPoF, which is only 20 ton, it is necessary to study the homogeneity of β-events in the entire FV of the CNO analysis and in the energy region of 210 Bi. The radial homogeneity has been studied by dividing the FV into 25 iso-volumetric shells. The angular homogeneity was studied by extending Fourier decomposition over a sphere surface, by projecting the spatial co-ordinates of the selected events on a sphere.
The final 210 Bi upper limit obtained through the estimation of the minimum 210 Po rate in the LPoF of Borexino, including both statistical and systematic contribution, is: R( 210 Bi) ≤ (11.5 ± 1.3) cpd/100 ton (1σ).

Sensitivity
The strength of the constraint on 210 Bi rate directly affects the precision of the CNO measurement. In order to evaluate quantitatively the precision of the CNO measurement of Borexino as a function of the 210 Bi rate constraint, a MC-based sensitivity study has been performed [4]. To assess the expected discovery potential to CNO neutrinos, a frequentist hypothesis test has been performed. Within this framework, two hypotheses have been considered: the null hypothesis H 0 , meaning that no CNO is assumed to exist, and the alternative hypothesis H 1 , that includes the presence of CNO. By indicating with L(H 0 ) and L(H 1 ) the two respective maximum likelihood values, the following likelihood ratio q 0 can be used as a test statistic [63]: Two sets of toy data have been produced, one with CNO injected one without CNO. Each dataset has been fit twice, i.e. with the H 0 and the H 1 hypotheses. The distribution of the test statistics q for the dataset without CNO injected, is called q 0 . The median of the q distribution for the dataset with CNO injected is called q med . Then, the p value of the q 0 distribution with respect to q med defines the discovery potential. Figure 10 reports the CNO median discovery significance assuming the HZ-SSM hypotheses and under different assumptions of pep and 210 Bi rate constraints, and assuming an exposure of 1000 days × 71.3 ton (93% of the Phase-III exposure). Neutrino interaction rates have been chosen according to HZ-SSM predictions. Background rates have been extracted from a MV fit on Phase-III data with the exception of the 210 Bi rate that has been set to 10 cpd/100 ton, a value similar to the upper limit estimated in the Phase-III data (see previous Section). This small difference does not influence the sensitivity, that is dominated by the precision of the constraint. The strong dependence of the sensitivity to CNO neutrinos on the strength of the external constraints on pep and 210 Bi rates is evident in Figure 10. By assuming the uncertainties similar to the one obtained in Phase-III data (Sections 3.4.1 and 3.4.2), a 5σ significance on the CNO neutrino signal can be reached by Borexino.

Spectral analysis
To extract the CNO neutrino interaction rate from data, a multivariate analysis has been performed, by simultaneously fitting the energy spectra between 0.320 MeV and 2.640 MeV and the radial distribution of events (Table 4), after all the selection procedure as already introduced in Section 3.2. The Phase-III dataset has been considered for the analysis. Apart from the constrained pep-neutrino and 210 Bi rates, and the fixed 8 B rate, all other species rates (including CNO neutrinos) have been left free to vary in the fit. The reference PDFs used to fit both the energy and radial distributions have been built by means of a complete Monte Carlo simulation. Considering only statistical uncertainty, the best fit returns a CNO interaction rate of R = 7.2 +2.9 −1.7 cpd/100 ton (68% confidence interval) with Borexino [3]. The effect of the fit configuration (such as the fit ranges) has been found negligible in the analysis. Since the fit heavily relies on the simulated MC PDFs for signal and backgrounds, a mismatch between data and simulations could potentially affect the CNO result and introduce a bias. To take this effect into account, as a possible source of systematics, a toy MC-based study has been performed. By generating several millions of pseudo-datasets, deforming signals and backgrounds every time and fitting with the same non-deformed PDFs, the impact on the CNO measurement has been evaluated. As possible sources of deformation, the following contributions have been considered: • Detector energy response, in terms of the scintillator energy scale (0.23%), non-uniformity (0.28%) and non-linearity (0.4%). The size of the deformation has been chosen based on the allowed values from calibration data and the 'standard candles' namely, 210 Po and 11 C.
• Deformation of the spectral shape of the cosmogenic 11 C isotope, induced by noise cuts not fully reproduced by MC (2.3%).
The final systematic contribution has been found to be +0.6 −0.5 cpd/100 ton, evaluated by comparing the CNO output from toy MC with and without injecting systematic distortions.  . From [4].
Counting analysis As a cross-check of the spectral analysis, a counting analysis has been performed on the same data sample. Data events have been counted inside an energy region where the CNO signal-tobackground ratio is maximised, corresponding to [0.780, 0.885] MeV. The pep-neutrino and 210 Bi numerical rate constraints, used in the spectral analysis, have been used in the counting analysis as well. The rate of 210 Bi has been symmetrically constrained and signals and backgrounds have been described by means of analytical functions. The CNO interaction rate has been extracted by subtracting all the background contributions, evaluated within a certain uncertainty, and then propagating those uncertainties. The final measured rate is (5.6±1.6) cpd/100 ton, confirming the presence of CNO at 3.5σ level. The quoted uncertainty considers both statistical and systematic terms. The uncertainty related to the energy response, in particular, is the dominant contribution in the overall balance. Figure 11 summarizes the Borexino result on the measured CNO neutrino interaction rate. The result of the spectral fit is reported in terms of log-likelihood profiles, with statistical uncertainty (dotted black line) and also with folded systematic contributions (solid black line). The best fit value is R = 7.2 +3.0 −1.7 cpd/100 ton (68% confidence interval), including systematic uncertainties. The inferred flux of CNO neutrinos at Earth is Φ = 7.0 +3.0 −2.0 × 10 8 cm −2 s −1 (68% confidence interval). The probability density function obtained from the counting analysis is also reported (solid red line). From the profiling of the log-likelihood, an exclusion of no-CNO hypothesis is achieved with 5.1σ significance. A further hypothesis test with 13.8 million Figure 11: Summary of CNO counting and spectral analysis results. Left, counting analysis bar chart. Right, CNO-neutrino rate negative log-likelihood (ln L) profiles, together with PDF from counting analysis and models predictions. From [3].

Final result
pseudo-datasets excludes the no-CNO hypothesis with 5.0σ at 99% confidence level.
The CNO neutrino measurement performed by Borexino is compatible with both HZ-SSM (0.5σ) and LZ-SSM (1.3σ) metallicity scenarios. A combined hypothesis test, including 7 Be and 8 B solar neutrino fluxes measured by Borexino, shows a preference for the HZ-SSM hypothesis at 2.1σ level.

Search for BSM physics
This section deals with searches beyond the Standard Model of particle physics with Borexino. Section 3.5.1 is dedicated to the search for Non-Standard neutrino interactions while Section 3.5.2 describes the procedure for the determination of the neutrino magnetic moment.

Non-standard neutrino interactions
Events selected in Borexino data for the solar neutrino analysis can be used to search for interactions not predicted by the Standard Model of Particle Physics. This class of phenomena is usually referred to as Non-Standard Interactions (NSI). The effect of NSI, regarding the solar neutrino flux, is to modify the electron neutrinos survival probability P ee , as well as the couplings between neutrinos and the scattered electrons. By analyzing Borexino Phase-II data [5], possible presence of NSI in the ν e e and ν τ e couplings have been studied. In particular, the analysis has looked for flavor-diagonal NC interactions. By considering the predictions from SSM, including oscillation effects, as well as both HZ-SSM and LZ-SSM metallicity scenarios, deviations in the measured spectrum have been searched for. Figure 12 reports the allowed regions for NSI extracted from Borexino data, considering ν e e and ν τ e couplings. Both regions in Figure 12 are compatible with null values of the ε parameters, which indicate the absence of NSI. In any case, the limits set by Borexino represent a significant reduction in the allowed parameter space, with respect to other experiments.
The same Phase-II dataset and analysis approach have been used to measure the value of the square sine of the Weinberg angle sin 2 θ W , by considering it as a free parameter in the analysis. From the likelihood profile, the best fit value results in: Allowed region for the NSI parameters ε R τ and ε L τ . Contours are reported for both HZ-SSM and LZ-SSM metallicity scenarios. From [5].
which is in agreement with results obtained from other neutrino-electron scattering experiments.

Neutrino magnetic moment
The flux of solar neutrinos detected by Borexino can be also used to derive indications of a neutrino anomalous magnetic moment. The presence of this condition would have the effect of modifying the neutrino-electron cross section and, consequently, the visible energy spectrum of solar neutrinos. The analysis of Borexino Phase-II data, including an additional component in the fit originating from an anomalous neutrino magnetic moment, has been performed in [6]. Figure 13 shows the likelihood profile used to estimate the limit on the neutrino magnetic moment.
The sum of the solar neutrino fluxes have been constrained using inputs from radiochemical gallium experiments. By including possible systematics effects, the derived limit for the neutrino magnetic moment is µ eff ν < 2.8 ×10 −11 µ B 6 . at 90% confidence level. The result is also free from uncertainties associated with predictions from the SSM neutrino flux. The limit on the effective value can be translated into a limit for magnetic moments of individual flavors, according to the measured values of flavor oscillation parameters. Limits on the value of the neutrino magnetic moment can be obtained from neutrino-electron elastic scattering experiments, either using solar [67,68] or reactor neutrinos [69,70]. The value obtained by Borexino is numerically the most stringent constraint on the neutrino magnetic moment.

Geoneutrinos
Geoneutrinos are mostly antineutrinos emitted as byproducts of radioactive decays inside the Earth. The fundamental quests of Earth sciences involve answering long-standing questions about the thermal, geodynamical, and geological evolutions of the Earth. Although geoneutrinos cannot answer all these questions, their measurement can serve as a unique tool in understanding the abundance of radioactive elements inside the Earth, especially in the inaccessible mantle. This can thereby help to determine the contribution of radiogenic heat to the total heat flux measured on Earth -a key parameter in many aspects of geosciences. This part of the review is divided into the following sections: Section 4.1 describes the structure, composition, and the heat flow of the Earth and it also presents geoneutrinos as a new tool for geoscience. The different components required for the geoneutrino analysis in Borexino are described in Section 4.2. The final results and geological interpretations of Borexino's recent geoneutrino measurement [7] are discussed in Section 4.3.

The Earth and geoneutrinos
This section is dedicated to the detailed structure of the Earth and its heat budget (Section 4.1.1) and the different models predicting the composition of the silicate part of the primordial Earth and thus, the Earth's radiogenic heat (Section 4.1.2). Geoneutrinos and their role in geosciences, in particular in the determination of the Earth's radiogenic heat, are discussed in Section 4.1.3.

The Earth's structure and heat budget
The Earth was formed by the accumulation of matter from the solar nebula over time [71,72]. All bodies with sufficient mass undergo the process of differentiation, i.e., the transformation of a homogeneous body into a layered structure. Primitive homogeneous Earth went through the first differentiation through segregation of metals, when the core separated from the silicate primitive mantle or bulk silicate earth (BSE). The BSE further differentiated into the present mantle and crust. The metallic core has Fe-Ni chemical composition and is expected to reach temperatures up to about 6000 K in its central parts. The inner core (∼1220 km radius) is solid due to high pressure, while the 2263 km thick outer core is liquid. The outer core has an approximate 10% admixture of lighter elements and plays a key role in the geodynamo process, which generates the Earth's magnetic field. The core-mantle boundary (CMB), a seismic discontinuity divides the core from the mantle. The mantle reaches temperature of about 3700 K at its deeper part. It is solid but viscous on long time scales, leading to mantle convection processes. This drives the movement of tectonic plates at the speed of few cm per year. At a depth of 400-700 km, the mantle is characterised by a transition zone, where a weak seismic-velocity heterogeneity is measured. The upper portion of the mantle contains the viscous asthenosphere on which the lithospheric tectonic plates are floating. These comprise the uppermost, rigid part of the mantle (i.e. the continental lithospheric mantle (CLM)) and the two types of crust: oceanic crust (OC) and continental crust (CC). The CLM is a portion of the mantle beneath the CC at a typical depth of ∼175 km [73]. The CC has a thickness of (34 ± 4) km [74] and is the most differentiated and heterogeneous layer, due to its complex history. It consists of igneous, metamorphic, and sedimentary rocks. The OC with (8 ± 3) km [74] thickness is created along the mid-oceanic ridges, where the basaltic magma differentiates from the partially melting mantle, up-welling towards the ocean floor.
The heat flow from the Earth's surface to the space results from a large temperature gradient across the Earth. The current best estimate for the total surface heat flux on Earth is H tot = (47 ± 2) TW [75]. Neglecting the small contribution (<0.5 TW) from tidal dissipation and gravitational potential energy released by the differentiation of crust from the mantle, the H tot is typically expected to originate from two main processes: (i) secular cooling H SC of the Earth, i.e. cooling from the time of the Earth's formation when gravitational binding energy was released due to matter accretion, and (ii) radiogenic heat H rad . The radiogenic heat of the present Earth arises mainly from the decays of isotopes with lifetimes comparable to, or longer, than the Earth's age (4.543 · 10 9 years): 232 Th (τ = 2.02 · 10 10 years), 238 U (τ = 6.447 · 10 9 years), 235 U (τ = 1.016 · 10 9 years), and 40 K (τ = 1.801 · 10 9 years) [76]. All these isotopes are labeled as heat-producing elements (HPEs). The relative contribution of radiogenic heat to the H tot is crucial in understanding the thermal conditions that existed during the formation of the Earth and the energy that is now available to drive the dynamical processes such as the mantle and outer-core convection. The convective urey ratio (UR CV ) quantifies the ratio of internal heat generation in the mantle over the mantle heat flux and is expressed as follows [77]: where H CC rad is the radiogenic heat produced in the continental crust, that is relatively well known to be 6.8 +1.4 −1.1 TW [74]. By adding the contribution from the oceanic crust and continental lithospheric mantle, the radiogenic heat from the lithosphere results to be 8.1 +1.9 −1.4 TW [7]. The mantle radiogenic heat is poorly constrained and ranges between 1.2 and 39.8 TW [7]. This is discussed further in Section 4.2.1. No radiogenic heat is expected to be produced in the metallic core. Preventing dramatically high temperatures during the initial stages of Earth formation, the present-day UR CV must be in the range between 0.12 to 0.49 [78].

Bulk Silicate Earth models
The Bulk Silicate Earth (BSE) models define the original chemical composition of the primitive mantle, including the abundances of HPE's and thus, the respective radiogenic heat. The elemental composition of BSE is obtained assuming a common origin for celestial bodies in the solar system. It is supported, for example, by the strong correlation observed between the relative (to Silicon) isotope abundances in the solar photosphere and in the CI chondrites, a special group of rare stony meteorites belonging to the carbonaceous chondrites [79]. Such correlations can be then assumed for the material from which the Earth was created. The BSE models agree with each other in the prediction of major elemental abundances (e.g. O, Si, Mg, Fe) within 10%. Uranium and Thorium abundances are assumed based on relative abundances in chondrites, and dramatically differ between different models. There are three main classes of BSE models: the cosmochemical, geochemical, and geodynamical models, as defined in [80,81]. The cosmochemical (CC) model [73] is characterised by a relatively low amount of U and Th, producing a total H rad = (11 ± 2) TW. This model assumes that the Earth is composed of enstatite chondrites. The geochemical (GC) model predicts intermediate HPE abundances for primordial Earth. It adopts the relative abundances of refractory lithophile elements as in CI chondrites, while the absolute abundances are constrained by terrestrial samples [82,83]. The geodynamical (GD) model shows relatively high U and Th abundances. It is based on the energy dynamics of mantle convection and the observed surface heat loss [84]. Additionally, an extreme fully radiogenic (FR) model can be obtained following the approach described in [85], where the terrestrial heat H tot of 47 TW is assumed to be entirely from radiogenic heat production H rad . Apart from these four main classes of models, there are also individual geological models used for comparison in the geological interpretations of Borexino's measurement discussed in Section 4.3. The predictions of these models are discussed in detail in [7].
Th/U ratio A global assessment of the Th/U mass ratio of the primitive mantle could provide information about the early evolution of the Earth and its differentiation. The most precise estimate of the planetary Th/U mass ratio has been refined to a value of M Th /M U = (3.876 ± 0.016) [86] and is relevant for the geoneutrino analysis as it will be discussed in the later sections. Significant deviations from this average value can be found locally, in the heterogeneous continental crust surrounding the detectors. The area surrounding the Borexino detector is characterised by a Th/U mass ratio ranging from ∼0.8 (carbonatic rocks) to ∼3.7 (terrigenous sediments) [87].
K/U ratio Potassium is the only semivolatile HPE. The bulk mass of Potassium predicted by different Earth models varies widely. Due to different possible scenarios, the K/U ratio predicted by different BSE models differs from 9700 to 16000 [81]. According to these ratios, the mantle radiogenic heat from 40 K varies between 2.6 and 4.3 TW. This results in an average contribution of 18% to the total mantle radiogenic power and is used in the determination of radiogenic heat from Borexino's geoneutrino measurement in Section 4.3.3.

Probing the Earth with geoneutrinos
Geoneutrinos are (anti)neutrinos emitted during the decay of the long-lived HPEs, discussed in the previous subsection, that can be summarised by the following equations: 238 U → 206 Pb + 8α + 8e − + 6ν e + 51.7 MeV (22) 235 U → 207 Pb + 7α + 4e − + 4ν e + 46.4 MeV (23) 232 Th → 208 Pb + 6α + 4e − + 4ν e + 42.7 MeV (24) 40 K → 40 Ca + e − +ν e + 1.31 MeV (89.3%) (25) 40 K + e − → 40 Ar + ν e + 1.505 MeV (10.7%) (26) In the β decays of 238/235 U and 232 Th chains and that of 40 K, all geoneutrinos are antineutrinos. Neutrinos are emitted only in the electron-capture decays of 40 K, having 10.7% branching ratio. In the radioactive decays of HPEs, the amount of released geoneutrinos and radiogenic heat are in a well-known ratio [equations (22) to (26)]. Thus, a direct measurement of the geoneutrino flux provides useful information about the composition of the Earth's interior. Consequently, it also provides an insight into the radiogenic heat contribution to the measured Earth's surface heat flux. Even though their small interaction cross section limits our ability to detect them, it makes them a unique probe of inaccessible innermost parts of the Earth. The first ideas to detect geoneutrinos started in the 1960s [88][89][90], further developed in 1984 [91], and the potential of liquid scintillator detectors to measure geoneutrinos was suggested in the 1990s [92,93]. It took several years to prove the methods feasible and the first investigation was performed by KamLAND in 2005 [94]. So far, only KamLAND [95,56] and Borexino [96][97][98]7] have been able to measure geoneutrinos.
The determination of the radiogenic component of Earth's internal heat budget has proven to be a difficult task, since an exhaustive theory is required to satisfy geochemical, cosmochemical, geophysical and thermal constraints, which are often based on indirect arguments. Therefore, direct geoneutrino measurements can be of great help. Geoneutrinos have also the potential to determine the mantle radiogenic heat, the key unknown parameter. This can be done by constraining the relatively-well known lithospheric contribution, as it will be shown in Section 4.3.2. Since the oceanic crust is depleted of HPEs and has a very small lithospheric contribution that can be easily determined, it would make the ocean floor an ideal environment for geoneutrino detection. Geoneutrino measurements can also contribute to the discussion about possible additional heat sources, whproposed by some authors. For example, stringent limits can be set on the power of a hypothetical Uranium natural georeactor suggested in [99][100][101][102] and discussed in Section 4.2.2 and 4.3.4.

Geoneutrino analysis
The various components required for the geoneutrino analysis with Borexino are discussed in this section. These include: the expected geoneutrino signal at LNGS (Section 4.2.1), the evaluation of other antineutrino signals (Section 4.2.2) and non-antineutrino backgrounds (Section 4.2.3), the data selection cuts (Section 4.2.4), the spectral fit used to extract the signal (Section 4.2.5), and the evaluation of various systematic uncertainties (Section 4.2.6).

Expected geoneutrino signal
The Earth shines in a flux of antineutrinos with a luminosity L ∼ 10 25 s −1 . Due to the IBD threshold discussed in Section 2.5, only 238 U and 232 Th geoneutrinos can be measured by LS detectors. For a detector placed on the continental crust, the expected 238 U and 232 Th geoneutrino flux is of the order of 10 6 cm −2 s −1 and is typically dominated by the crustal contribution. The differential flux of geoneutrinos emitted from isotope i = ( 238 U, 232 Th) and expected at LNGS location r is calculated using the following expression: where ε ν (i) is the specific antineutrino production rate for isotope i per 1 kg of naturally occurring element and Eν is geoneutrino energy. The geoneutrino energy spectra dn(i;Eν) dEν is normalised to one. The electron-flavour survival probability P ee after the propagation of geoneutrinos from a geological reservoir located at r to the detector is calculated considering oscillations in vacuum. The matter effect is estimated to be of the order of 1% [103], i.e., much less than other uncertainties involved in the geoneutrino signal prediction. The average survival probability P ee is 0.55. Despite the dependence of P ee w.r.t. energy, the geoneutrino energy spectrum is unchanged since the geoneutrinos undergo several oscillations before reaching the detector. The ρ( r ) is the density of the voxel emitting geoneutrinos and it is taken from geophysical models of lithosphere [74] and mantle [104]. The abundances a(i; r ) of isotope i are expressed per mass unit of rock. The integration is done over the whole volume of the Earth, considering geological constraints of the main HPE reservoirs.
The antineutrino signal can be conveniently expressed in Terrestrial Neutrino Units (TNU). 1 TNU corresponds to 1 antineutrino event detected via IBD (Section 2.5) over 1 year by a detector with 100% detection efficiency containing 10 32 free target protons (roughly corresponding to 1 kton of LS). To convert the differential geoneutrino flux dΦ(i;Eν, r) dEν to geoneutrino signal S(i) expressed in TNU, it is necessary to perform an integration over the geoneutrino energy spectra, considering the energy dependence of the IBD cross section [105], and 10 32 target protons in 1 year measuring time. The geoneutrino energy spectrum that can be observed by Borexino extends from 1.8 MeV (IBD threshold) to 3.7 MeV. Note that for a reference oscillated flux of 10 6 cm −2 s −1 , the geoneutrino signals from 238 U and 232 Th are S(U) = 12.8 TNU and S(Th) = 4.04 TNU, respectively. Considering the specific antineutrino production rates ε ν (U, Th) 7 one can calculate the signal ratio R S for a homogeneous reservoir characterised by a fixed a(Th)/a(U) ratio (where a denotes the abundance): This signal ratio thus depends on the composition of the reservoir. Adopting the CI chondrites a(Th)/a(U) = 3.9 for the bulk Earth, the signal ratio R s = 0.27 can be extracted. Geophysical and geochemical observations of the lihtosphere constrain the a(Th)/a(U) = 4.3 [7], implying a signal ratio of 0.29. As a consequence of maintaining the global chondritic ratio of 3.9 for the bulk Earth, the inferred mantle ratio a(Th)/a(U) results to be 3.7, which corresponds to a R S = 0.26. This is used for the extraction of mantle signal explained in Section 4.3.2.
The expected geoneutrino signal in Borexino S(U+Th) can be expressed as the sum of three components and this is shown in detail in Figure 14.
• S LOC (U+Th), the local crust (LOC) signal produced from the 492 km × 444 km crustal area surrounding LNGS, • S FFL (U+Th), the signal from the far field lithosphere (FFL), which includes the continental lithospheric mantle (CLM), i.e. the brittle portion of the mantle underlying the CC, and the remaining crust obtained after the removal of the LOC.
• S mantle (U+Th), the signal from the mantle.
The sum of S LOC (U+Th) and S FFL (U+Th) is the expected signal from the bulk lithosphere (U + Th) and is evaluated to be 25.9 +4.9 −4.1 TNU. This constraint is used for the extraction of mantle signal at LNGS, described in Section 4.3.2. The mantle geoneutrino signal and the mantle radiogenic heat are highly debated topics in Earth sciences. The distribution of HPEs in the mantle is not known and the assumed different scenarios influence the strength of the mantle signal. There is a high scenario, which assumes homogeneous HPE distribution, an intermediate scenario, that assumes HPEs distributed in two layers, a lower enriched mantle (EM) and an upper depleted mantle (DM) separated at 2180 km of depth, and a low scenario, where the HPEs are placed just above the CMB. Assuming these scenarios, the predictions of the mantle geoneutrino signal S mantle (U+Th) vary between 0.9 and 33.0 TNU for different BSE models. Using the relatively wellknown contribution from the bulk lithosphere, and the various predictions for the mantle signal, the total geoneutrino signal at LNGS S(U+Th) varies between 28.

Other antineutrino signals
The antineutrino backgrounds relevant for the geoneutrino analysis (Section 4.2) are described below. The main source of background in geoneutrino detection is the production of electron antineutrinos by nuclear power plants, the strongest man-made antineutrino source. In addition, atmospheric neutrinos can have a small impact. This has been treated as a systematic source of uncertainty in the recent measurement [7]. Another source of antineutrinos can be a hypothetical georeactor inside the Earth, searched for in a separate analysis (Section 4.3.4).
Reactor antineutrinos Many nuclei, produced in the fission process of nuclear fuel, decay through β-processes, with the consequent emission of electron antineutrinos, the so-called reactor antineutrinos. Their energy spectrum extends up to 10 MeV, well beyond the end point of the geoneutrino spectrum  [76]. An accurate determination of the expected signal and spectrum of reactor antineutrinos requires a wide set of information, spanning from the characteristics of nuclear cores to neutrino properties and is discussed in detail in [7]. Various experimental results from reactor antineutrino experiments show that the IBD positron energy spectrum deviates significantly from the spectral predictions of Mueller et al. 2011 [106] in the energy range between 4 -6 MeV, i.e. the so-called "5 MeV excess" and also an overall deficit with respect to the prediction. These features are treated as a systematic uncertainty (Section 4.2.6) using the Daya Bay high precision measurement [107]. The expected signal from reactor antineutrinos can be expressed in TNU, assuming 100% detection efficiency for a detector containing 10 32 target protons and operating continuously for 1 year. The uncertainties related to reactor antineutrino production, propagation, and detection processes are estimated using a Monte Carlo-based approach discussed in [103]. The expected reactor antineutrino signal with and without the Atmospheric neutrinos Atmospheric neutrinos can act as a potential background source for the geoneutrino measurement and this was studied for the first time in [7]. Atmospheric neutrinos originate in sequential decays of π ± and K ± mesons and µ ± leptons produced in cosmic rays' interactions with atmospheric nuclei. The flux of atmospheric neutrinos contains both neutrinos and antineutrinos, and the muon flavour is roughly twice abundant with respect to the electron flavour. The process of neutrino oscillations then alters the flavour composition of the neutrino flux passing through the detector. Atmospheric neutrinos interact in many ways with the nuclei constituting the Borexino scintillator. The most copious isotopes in the Borexino scintillator are 1 H (6.00 × 10 31 /kton), 12 C (4.46 × 10 31 /kton), and 13 C (5.00 × 10 29 /kton). Besides the IBD reaction itself, there are many reactions with 12 C and 13 C atoms that may, in some cases, mimic the IBD interactions. They have the form of ν + A → ν(l) + n + · · · + A , where A is the target nucleus, A is the nuclear remnant, l is charged lepton produced in CC processes, n is the neutron, and dots are for other produced particles like nucleons (including additional neutrons) and mesons (mostly π and K mesons). A dedicated simulation code was developed to precisely calculate this background in Borexino [7]. The expected number of IBD-like interactions due to atmospheric neutrinos in the energy window used for the analysis was evaluated via MC to be 6.7 ± 3.7 events in the analysed period.
Georeactor A possible existence of a georeactor, i.e. a natural nuclear fission reactor in the Earth's interior, was first suggested by Herndon in 1993 [99]. Since then, several authors have discussed its possible existence and its characteristics. Different models suggest the existence of natural nuclear reactors at different depths [100][101][102]. New upper bounds on the power of a potential georeactor are discussed in Section 4.3.4. In order to be able to set such limits for different hypothetical locations of the georeactor, the expected antineutrino spectra of a 1 TW point-like georeactor were calculated. It was assumed to operate continuously during the entire analysis period with the constant power fractions of fuel components, as suggested in [108] ( 235 U : 238 U ≈ 0.76 : 0.23). The energy released and the antineutrino spectra per fission were calculated similar to reactor antineutrinos, using only the flux parametrisation of [106]. Three different depths were considered for the georeactor location: 1) GR1: the Earth's centre (d = R Earth ), 2) GR2: the CMB at just below the LNGS site (d = 2900 km), and 3) GR3: the CMB on the opposite hemisphere (d = 2R Earth − 2900 = 9842 km). In the calculation of the survival probability, the matter effect is taken into account by assuming a constant Earth density (density variation in between the crust and the core causes about 3% signal variation). Unfortunately, Borexino's energy resolution does not allow to distinguish the spectral differences due to oscillations. The errors were calculated via the method similar to that used for reactor antineutrinos. The expected georeactor signals for a 1 TW georeactor located at GR1, GR2, and GR3 are 43.1± 1.3 TNU, 8.9 ± 0.3 TNU, and 3.7 ± 0.1 TNU, respectively.

Non-antineutrino backgrounds
There are various non-antineutrino backgrounds that can mimic an IBD and should therefore be estimated for the geoneutrino analysis. The three important non-antineutrino backgrounds that are constrained in the final spectral fit (Section 4.2.5) -cosmogenic 9 Li, accidental, and (α, n) backgrounds -are described in more detail below. In addition to these, there are various other minor backgrounds that are also estimated. These include: backgrounds due to untagged muons, fast neutrons from the surrounding rocks and the water tank, (γ, n) reactions, fission reactions in PMTs, and 214 Bi-214 Po coincidences due to 222 Rn contamination during the water extraction period. Table 8 summarizes the expected number of events from all non-antineutrino backgrounds passing the IBD selection cuts discussed in Section 4.2.4.
Cosmogenic 9 Li Cosmic muons reaching Borexino can produce many spallation products that can in turn produce (β − + n) pairs as shown in Table 2. Such decays are indistinguishable from the IBD events, since the prompt is represented by the β − and the delayed is due the neutron, exactly as in the IBD. The hadronic background for geoneutrino analysis is dominated by the 9 Li isotope. It has a Q-value of 13.6 MeV, a lifetime of 257.2 ms [17] and decays according to the following process: 9 Li → e − +ν e + n + 2α.
The 9 Li background can be estimated by studying the events that follow immediately after a cosmic muon, using the same IBD selection cuts described in Section 4.2.4. The residual background after the time and space muon vetoes can then be estimated using the measured lifetime of the 9 Li prompts and their distance to the muon, respectively. In order to account for other cosmogenic isotopes such as 8 He and 9 B, which can also be present in the muon daughter sample, the measured lifetime from Borexino data (260 ± 21) ms [7] is used to evaluate this background instead of the 9 Li lifetime.  Accidental background Accidental coincidences happen due to the random coincidences of promptlike and delayed-like events in Borexino. This is an important background for the geoneutrino analysis and depends on the selection cuts applied. They are mostly due to external backgrounds, since their reconstructed positions are mostly near the IV. They are usually searched for in an off-time window of 2 s to 20 s after the prompt and scaled back to the IBD time window of 1.28 ms (Section 4.2.4). Such a long time is used to maximize the number of selected accidental coincidences and thus to reduce the error with which the shape and rate of this background can be constrained.
(α, n) coincidences Decays along the chains of 238 U, 235 U, and 232 Th can produce α-particles, which in turn can give rise to (α, n) interactions that can mimic IBD coincidences. In Borexino, this is only due to the 210 Po isotope; while it is a part of the 238 U chain, it is found fully out of the secular equilibrium and its contamination is several orders of magnitudes higher than the rest of the chain. The α particle can then interact with other isotopes, mostly 13 C in the LS: The cross section of this interaction is 200 mb [109] and the produced neutron can have energies up to 7.3 MeV, almost indistinguishable from the IBD delayed. There are three different possibilities for a promptlike interaction, which are described in detail in [7]. This background is estimated using the amount of 210 Po events (tagged on an event-by-event basis via pulse shape discrimination techniques, see Section 2.2) in the entire analysed exposure, the probability of an (α, n) interaction to produce IBD-like coincidences passing the selection cuts, and the neutrino-yield of the reaction in pseudocumene.

Selection cuts
The various criteria used to select IBD candidates for the geoneutrino analysis in Borexino are described below. The detection efficiency of geoneutrinos for these selection cuts was estimated through MC to be (87.0 ± 1.5)%.

Muon vetoes
Muon vetoes are applied after both internal and external muons in Borexino, in order to reduce the cosmogenic background. Generally, a 2 ms veto is applied after external muons, similarly to the solar neutrino analysis, in order to remove fast neutrons. In the geoneutrino analyses performed until 2015, a conservative 2 s veto (8 times the lifetime of 9 Li) was applied after internal muons to eliminate muon daughters. In the improved analysis performed in 2019 [7], 4 different kinds of more complex space-time vetoes were applied, depending on the type of muon. This reduced the exposure loss from 10-11% in the previous analyses to 1.2% in the latest analysis.
Energy cuts Energy cuts are applied on both prompt and delayed, after the consideration of the IBD threshold and the n-capture energy peaks described in Section 2.5. The energy spectrum of prompt starts at ∼1 MeV, which corresponds to the two 0.511 MeV gammas. The threshold energy of the prompt is chosen as 408 p.e., corresponding to about 0.8 MeV. There is no upper limit on the energy of the prompt. The energy of the delayed can either be due to the neutron capture on 1 H (2.2 MeV gamma) or on 12 C with 1.1% probability (4.95 MeV gamma). However, at large radii, gammas can partially deposit their energy in the buffer, decreasing the visible energy. Consequently, the gamma-peak develops a low-energy tail and even the peak position can be shifted to lower values. Therefore, in the latest geoneutrino analysis, the energy ranges 700-1300 p.e and 1300-3000 p.e. were used for IBD selection, to account for the n-captures on 1 H and 12 C, respectively. However, the lower energy threshold was increased to 860 p.e. for the water extraction period between 2010 and 2011, because of increased 222 Rn contamination.

Time space correlation
The cuts applied on the time and distance between the prompt and the delayed are important background suppression cuts. The IBD events can either have both the prompt and delayed clusters in the same 16 µs DAQ window (double cluster events) or have them in separate DAQ windows as single cluster point-like events. A time coincidence window of 20-1280 µs is used for IBD selection of single cluster events. The lower threshold guarantees that the delayed event can trigger after the DAQ deadtime of 3-4 µs after the prompt event and the higher threshold corresponds to 5 times the measured neutron capture time of (254.5 ± 1.8) µs [20]. In the latest geoneutrino analysis, the time window 2.5-12.5 µs was also considered, in order to include double cluster events. The lower threshold was tuned on reactor antineutrino MC and guarantees that even for the highest energy prompt, the light from the prompt cannot alter the delayed event after 2.5 µs. The higher threshold was chosen by studying the variable position of the prompt event inside the DAQ window in the entire analysed period. In addition, it also allows the delayed event to have a cluster duration up to 2.5 µs before the end of the DAQ window.
The reconstructed distance between the prompt and the delayed is larger than the distance between their points of production due to different reasons [7]. The important factors to be considered before choosing the spatial cut are the MC efficiency and the accidental background. The efficiency drops quickly below 1 m and the accidental background steadily increases with increasing distance. This spatial cut was optimised to be 1.3 m, using MC-based sensitivity studies.
Pulse-shape discrimination Pulse shape discrimination is applied to the delayed in order to distinguish it from other α-like coincidences. The MLP parameter mentioned in Section 2.2 was used for the latest geoneutrino analysis. A cut > 0.8 was applied on the delayed to distinguish it from (α+γ) decays of 214 Po 8 , arising from 222 Rn contamination during the water extraction period.
Dynamical Fiducial Volume cut The vessel shape of Borexino changes in time due to the presence of a small leak and it can be reconstructed, as discussed in Section 2.1. The Dynamical Fiducial Volume (DFV) cut is applied in the geoneutrino analysis based on the distance of the prompt from the IV. The cut was optimised to be 10 cm, using MC-based sensitivity studies in the latest geoneutrino analysis. This distance is sufficient enough to suppress background events near the IV, and to account for the uncertainty due to the IV shape reconstruction. The enlarged fiducial volume in the recent analysis resulted in a 15 % increase in exposure when compared to the previous analyses.

Multiplicity cut
The multiplicity cut is applied to suppress neutron-neutron or buffer muon-neutron pairs that can enter the IBD selection due to undetected muons. This cut requires that there are no high energy events (> 400 p.e.) present in a ±2 ms time window (about 8 times neutron capture time) around the prompt or delayed. The corresponding exposure loss due to accidental coincidences of IBD candidates with >400 p.e. events within 2 ms is of the order of 0.01%, which is negligible for this analysis.

Spectral fit
After the selection of IBD candidates, the geoneutrino signal is extracted from the unbinned-likelihood spectral fit of the charges of all prompts: where N p pe is the vector of individual prompt charges N p i pe , and index i runs from 1 to N IBD , i.e. the total number of IBD candidates. The symbol θ indicates the set of the variables with respect to which the function is maximised, namely the number of events of signal and backgrounds. The shapes of all spectral components used in the fit are taken from the Probability Distribution Density Functions (PDFs) constructed using MC, with the exception of the accidental background which can be measured with sufficient precision using off-time coincidences in data. The result of the fit is the number of events due to each spectral component.
The spectral fit is generally performed in the range 408-4000 p.e. (≈0. . The number of geoneutrinos is always kept free. One way of doing it is by having one free fit parameter for geoneutrinos, when using the PDF in which the 232 Th and 238 U contributions are summed and weighted according to the chondritic mass ratio of 3.9, corresponding to a signal ratio of 0.27 (Section 4.2.1). Alternatively, 232 Th and 238 U contributions can be fit as two independent contributions. Additional combinations are also possible. The number of reactor antineutrino events is typically kept free. This is an important cross-check of the ability of Borexino to measure electron antineutrinos and can be verified by comparing the fit results to the well-known prediction of reactor antineutrinos in Section 4.2.2. The three non-antineutrino backgrounds (Section 4.2.3) are constrained using additional multiplicative Gaussian pull terms in the likelihood function of Equation 31. Using the known exposure and detection efficiency, the number of detected geoneutrinos and reactor antineutrinos can then be expressed in the units of TNU.

Systematic sources of uncertainty
The systematic uncertainty on the measured geoneutrino signal is very small compared to the statistical uncertainty. The total uncertainties on the geoneutrino and reactor antineutrino signals were evaluated to be +5.2 −4.0 % and +5.1 −5.5 %, respectively. The different systematic sources are discussed below and summarised in Table 9.
• Atmospheric neutrinos: Atmospheric neutrinos as a source of background for geoneutrinos were discussed in Section 4.2.2. The uncertainty on the expected number of atmospheric neutrino events is estimated to be around 50%. Therefore, the atmospheric neutrino PDF was included in the standard spectral fit, in order to evaluate the impact of this systematic source.
• Shape of reactor antineutrino spectra: The standard spectral fit was performed using the MC PDF of reactor antineutrinos without the "5 MeV excess" which was discussed in Section 4.2.2. The effect of this assumption was studied by using PDFs with and without the "5 MeV excess" in the spectral fit.
• Inner vessel shape reconstruction: The systematic uncertainty on the FV definition due to the 5 cm error on the IV shape reconstruction is negligible. However, there is a systematic uncertainty associated with the selection of the IBD candidates using the DFV cut, which was evaluated by smearing the  Table 9: Summary of the different sources of systematic uncertainty in the geoneutrino and reactor antineutrino measurement. Different contributions are summed up as uncorrelated. From [7].
distance-to-IV of each IBD candidate with a Gaussian function (σ = 5 cm). Then, the spectral fit was performed on the newly selected candidates with the DFV cut and repeated 50 times.
• MC efficiency: The major source of uncertainty for the MC efficiency arises from the event losses close to the IV edges, especially near the south pole, because of the combined effect of a large number of broken PMTs and the IV deformation.
• Position reconstruction: Since the events are selected inside the DFV based on the reconstructed position, the uncertainty in the position reconstruction of events affects the error on the fiducial volume, and thus, on the resulting exposure.

Results and geological interpretations
In the period between December 9, 2007 and April 28, 2019, corresponding to 3262.74 days of data acquisition, 154 IBD candidates were found [7]. The exposure of (1.29 ± 0.05) × 10 32 protons × year represents an increase by a factor of two with respect to the previous analysis in 2015 [98]. The time, spatial, and energy distributions of the IBD candidates were compatible with the expectations. The following subsections discuss the geoneutrino signal measured by Borexino at LNGS (Section 4.3.1), the extraction of the mantle geoneutrino signal (Section 4.3.2), the calculation of the radiogenic heat and Urey ratio from Borexino measurement and its comparison to different BSE model predictions (Section 4.3.3), and the upper limits on a hypothetical georeactor (Section 4.3.4).

Geoneutrino signal at LNGS
The spectral fit described in Section 4.2.5 was performed on the prompt charges of the 154 IBD candidates. The three major non-antineutrino backgrounds, namely, the cosmogenic 9 Li background, the (α, n) background from the scintillator, and accidental coincidences were included in the fit and were constrained according to values in Table 8 with Gaussian pull terms. The geoneutrino and reactor antineutrino contributions were left free in the fit and the geoneutrino MC PDF was constructed assuming a Th/U chondritic ratio of 3.9. This resulted in 52.6 +9.4 For comparison, the star shows the best fit performed assuming the 238 U and 232 Th contributions as free and independent fit components. From [7].
spectral fit and the 1, 3, 5, and 8σ contours for the number of geoneutrinos versus reactor antineutrinos are shown in Figure 15. Compatible results were obtained when the spectral fit was performed by leaving the 238 U and 232 Th contributions free, without fixing them to the chondritic ratio. Unfortunately, Borexino does not have the sensitivity to measure the Th/U ratio. The geoneutrino signal measured by Borexino over the years is compared with the latest measurement in Figure 16(a). The signal along with the 68% C.I. is compared to different BSE models in Figure 16(b). The bulk lithosphere contribution for all these models is the same as discussed in Section 4.2.1, while the mantle signal varies according to the different model predictions.

Extraction of mantle signal
The mantle signal was extracted from the spectral fit after constraining the contribution from the bulk lithosphere to be 28.8 +5.5 −4.6 events corresponding to the expected signal of 25.9 +4.9 −4.1 TNU mentioned in Section 4.2.1. The corresponding MC PDF was constructed from the PDFs of 232 Th and 238 U geoneutrinos. They were scaled with the lithospheric Th/U signal ratio equal to 0.29, that is based on geological observations. The MC PDF used for the mantle was also constructed from the 232 Th and 238 U PDFs, but the applied Th/U signal ratio was 0.26. This procedure maintains the chondritic Th/U mass ratio of 3.9 for the bulk Earth and infers this ratio in the mantle to be 3.7 (Section 4.2.1). The mantle signal, as well as the reactor antineutrino contribution were kept free in the fit. The data points and the different PDFs of the fit are shown in Figure 17. After considering the systematic uncertainties, the final mantle signal obtained was 21.2 +9.5 −9.0 (stat) +1.1 −0.9 (sys) TNU. The statistical significance of the mantle signal was studied using MC pseudo-experiments with and without a generated mantle signal. This way, the null-hypothesis of the mantle signal has been rejected with 99.0% C.L. (corresponding to 2.3σ significance).

Radiogenic heat and convective Urey ratio
The mantle geoneutrino signal S mantle (U+Th) is linearly proportional to the mantle radiogenic power H mantle rad (U+Th), that in turn scales with the 238 U and 232 Th masses in the mantle M mantle (U) and M mantle (Th), respectively. This relationship is expressed by the following equation: where h(U) = 98.5 µW/kg and h(Th) = 26.3 µW/kg [110] are the 238 U and 232 Th specific heats. The value 3.7 is the assumed Th/U mass ratio in the mantle as discussed in the previous subsection. The scaling factor β depends only on 238 U and 232 Th distributions in the mantle. For Borexino, a central value of β centr = 0.86 TNU/TW represents the average between the assumptions of a homogeneous mantle (β high = 0.98 TNU/TW) and of the HPE-rich layer just above the CMB (β low = 0.75 TNU/TW). This relation is shown in Figure 18 for the range of mantle radiogenic power predicted by various groups of BSE models discussed in Section 4.   Figure 17: Spectral fit of the charge spectrum of prompt IBD candidates to extract the mantle signal after constraining the contribution of the bulk lithosphere. The grey shaded area shows the summed PDFs of all the signal and background components. From [7].
In order to compare the Borexino estimate of the radiogenic power with the Earth's total surface heat flux of H tot = (47 ± 2) TW [75], the contribution from 40 K must be considered. Firstly, assuming the contribution from 40 K to be 18% of the total mantle radiogenic heat (Section 4.1.1), the total radiogenic mantle signal can be expressed as H mantle rad (U + Th + K) = 30.0 +13.5 −12.7 TW. Through the further addition of the lithospheric contribution H LSp rad (U+Th+K) = 8.1 +1.9 −1.4 TW, the 68% C.I. for the Earth's radiogenice heat can be obtained as H rad (U+Th+K) = 38.2 +13.6 −12.7 TW. Figure 19(a) compares the decomposition of the Earth's total surface heat flux into radiogenic heat and the heat coming from secular cooling H SC , when performed for various BSE models and the Borexino measurement. A preference is found for models with relatively high radiogenic power, which indicates a cool initial environment at early formation stages of Earth and assumes the fraction of heat coming from secular cooling H SC to be small. However, no model can be excluded at 3σ level.
The total radiogenic heat estimated by Borexino can be used to extract the convective Urey ratio, according to Equation (21) in Section 4.1.1. The resulting value of UR CV = 0.78 +0. 41 −0.28 is compared to the UR CV predicted by different BSE models in Figure 19(b). The Borexino geoneutrino measurement constrains the mantle radiogenic heat power at a 90(95)% C.L. to be H mantle rad (U+Th) > 10(7) TW and H mantle rad (U+Th+K) > 12.2(8.6) TW and the convective Urey ratio to be UR CV > 0.13(0.04).

Limits on a hypothetical georeactor
Limits on a hypothetical georeactor, mentioned in Section 4.2.2, can be obtained by constraining the number of expected reactor antineutrino events in the spectral fit, since both these components have similar spectral shapes. The georeactor PDFs for the fit were generated for three different positions. However, their shapes were practically identical and Borexino does not have any sensitivity to distinguish them. The geoneutrino (Th/U fixed to chondritic mass ratio of 3.9) and georeactor contributions were left free in the fit. This resulted in an upper limit of 18.7 TNU, at 95% C.L., on a georeactor signal. Considering the predicted georeactor signal expressed in TNU, for a 1 TW georeactor in different locations [7], Borexino has excluded the existence of a georeactor with a power greater than 0.5/2.4/5.7 TW at 95% C.L., assuming its distance from the detector to be 2900 km (CMB below LNGS) /6371 km (Earth's centre) /9842 km (CMB on opposite hemisphere), respectively.  Figure 18: Mantle geoneutrino signal expected in Borexino as a function of 238 U and 232 Th mantle radiogenic heat: the area between the red and blue lines denotes the full range allowed between a homogeneous mantle (high scenario) and a unique rich layer just above the CMB (low scenario). The slope of the central inclined black line (β centr = 0.86 TNU/TW) is the average of the slopes of the blue and red lines. The blue, green, red, and yellow ellipses are calculated using the 238 U and 232 Th mantle radiogenic power (with 1σ error) according to different BSE models: Cosmochemical (CC) model, Geochemical (GC) model, Geodynamical (GD) model, and Fully radiogenic (FR) model. For each model darker to lighter shades of respective colours represent 1, 2, and 3σ contours. The black horizontal lines represent the mantle signal measured by Borexino: the median mantle signal (solid line) and the 68% coverage interval (dashed lines). From [7].  Figure 19: (a) Decomposition of the Earth's total surface heat flux H tot = (47 ± 2) TW (horizontal black lines) into its three major contributions: lithospheric radiogenic heat H LSp rad (brown), mantle radiogenic heat H mantle rad (orange), and secular cooling H SC (blue). The labels on the x axis identify different BSE models, while the last bar labeled BX represents the Borexino measurement. The lithospheric contribution is the same for all bars. (b) Comparison of Borexino constraints (horizontal band) with predictions of the BSE models (points with ±3σ error bars) for the convective Urey ratio UR CV (Equation 21), assuming the total heat flux H tot = (47 ± 2) TW. The blue, green, and red colours represent different BSE models (CC, GC, and GD, respectively). From [7].

Conclusions and outlook
In this paper the latest Borexino measurements on neutrinos from the Sun and Earth were discussed, from highlighting the key elements of the analyses up to the discussion and interpretation of the results. The success of Borexino in the measurement of low-energy neutrinos is primarily based on the extreme radio-purity of the liquid scintillator.
Borexino measures solar neutrinos via neutrino-electron elastic scattering. The direct result of the solar neutrino analysis are the interaction rates (given for zero-threshold) for neutrinos produced by different nuclear reactions. Borexino has performed a complete spectroscopy of the pp chain solar neutrinos [1]. By performing a multivariate spectral fit in the Low Energy Region, the rates of neutrinos from initial proton-proton fusion pp( +8.7 −10.6 %), the three-body proton-electron-proton fusion pep( +16.0 −17.4 % for HZ-SSM and +14.7 −16.3 % for LZ-SSM), and the electron-capture decay of 7 Be( +2.4 −2.7 %) were extracted. The numbers in parentheses indicate the total experimental error summing quadratically the statistical and systematic uncertainties. The constraint on the CNO rate used in the fit, based on the HZ-SSM and LZ-SSM predictions, influences only the pep neutrinos. In both cases however, the absence of the pep reaction in the Sun was rejected with > 5σ significance. The rate of scattered electrons above 3 MeV due to 8 B( +7.1 −7.7 %) solar neutrinos interactions was extracted via a radial fit without any assumption on the form of the energy spectrum and thus, oscillation parameters. Borexino is the only experiment that measured all (with the exception of hep neutrinos) neutrinos from the pp chain. While the Low Energy Region results are the most precise measurements existing in the world, the 8 B result is less precise than the measurements of large volume water Cherenkov detectors, but is compatible with them. Borexino measurements provide a direct determination of the relative intensity of the two primary terminations of the pp chain (pp-I and pp-II) and, assuming standard three-neutrino oscillations, an indication that the temperature profile in the Sun is more compatible with SSM models that assume high surface metallicity. Assuming the SSM prediction for the solar fluxes to hold, the survival probability of solar electron neutrinos can be determined from the measured rates: by comparing its values at different energies, Borexino probes simultaneously the neutrino flavour-conversion paradigm, both in vacuum and in matter-dominated regimes. The vacuumonly hypothesis is disfavoured at 98.2% C.L. Borexino also confirmed the solar origin of the measured signal assigned to 7 Be solar neutrinos, by observing the expected seasonal variation of the respective rate [2], induced by the Earth elliptical orbit around the Sun. Measurement of solar neutrinos also help in constraining the non-standard neutrino interactions [5] and placing a stringent limit on the effective neutrino magnetic moment [6].
Borexino provided the first experimental evidence at 5σ significance of neutrinos produced in the CNO fusion cycle in the Sun [3]. This was achieved by a multivariate spectral fit performed with a constraint on the rates of pep solar neutrinos and 210 Bi internal background. The pep rate was constrained with 1.4% precision based, on theoretical expectations and a global fit to all solar data, excluding the Borexino dataset used in this analysis. An upper limit constraint was placed on 210 Bi rate, obtained via a fit of the 210 Po distribution in the Low Polonium Field. This procedure was made possible, thanks to the thermal stabilisation of the Borexino detector during the Phase III, that minimised the convective currents bringing 210 Po from the nylon vessel holding the scintillator to the fiducial volume of the analysis. Thus, in the region of the scintillator free from convection, 210 Po rate approaches the rate of its parent 210 Bi nuclei. In the CNO cycle, the fusion of Hydrogen is catalysed by Carbon, Nitrogen, and Oxygen, and so its rate as well as the flux of emitted CNO neutrinos depend directly on the abundance of these elements in the solar core. Borexino result is compatible with both HZ-SSM and LZ-SSM predictions, but paves the way towards a direct measurement of the solar metallicity using CNO neutrinos. In addition, Borexino measurements quantify the relative contribution of CNO fusion in the Sun to be of the order of 1%. In massive stars, however, CNO is the dominant process of energy production and thus, the primary mechanism for the stellar conversion of Hydrogen into Helium in the Universe.
Geoneutrinos, antineutrinos from the decays of long-lived radioactive elements inside the Earth, can be exploited as a new and unique tool to study our planet. Only two experiments in the world, KamLAND and Borexino, have detected geoneutrinos so far. Geoneutrinos are detected via the Inverse Beta Decay (IBD) on protons. The 1.8 MeV kinematic threshold of this interaction allows to measure the high-energy part of geoneutrinos emitted along 238 U and 232 Th chains, while it leaves 40 K geoneutrinos completely unreachable for the present-day technology. Inverse Beta Decay has about two orders of magnitude higher cross section than the elastic scattering, a huge advantage for the detection of geoneutrinos with the flux about four orders of magnitude lower than the pp solar neutrino flux. In addition, it provides a unique background-suppressing topology of the fast space-time coincidences between a prompt and a delayed signal. Borexino has presented its latest geoneutrino measurement in 2019 [7]. The geoneutrino signal was extracted through a spectral fit of the prompt events, related to the energy of incident antineutrinos. Geoneutrino and reactor antineutrino contributions were kept as free fit parameters, while non-antineutrino backgrounds were constrained to the values estimated independently. Thanks to both more acquired data and improved analysis techniques in an enlarged fiducial volume, the updated measurement reached +18.3 −17.2 % total precision, assuming the same Th/U mass ratio as found in chondritic CI meteorites. Antineutrino background from reactors was fit unconstrained and found compatible with the expectations. The nullhypothesis of observing a geoneutrino signal from the Earth's mantle was excluded at a 99.0% C.L., when exploiting detailed knowledge of the local crust near the experimental site. The measured mantle signal was then converted to mantle radiogenic heat from decays of Uranium and Thorium, assuming a range of geological models describing the distribution of these elements in the mantle: from a homogeneous distribution up to an assumption of an enriched layer at the core-mantle boundary. The measured mantle signal is compatible with different geological predictions, however there is a ∼2.4σ tension with those Earth models, which predict the lowest concentration of heat-producing elements in the mantle. In addition, by constraining the number of expected reactor antineutrino events, the existence of a hypothetical georeactor at the centre of the Earth having power greater than 2.4 TW was excluded at 95% C.L.
To summarize, this paper reviewed the latest Borexino measurements of solar and geo neutrinos-results that certainly mark the history of neutrino physics. The Borexino collaboration is analysing the latest data taken with the detector, featuring exceptional radio-purity and thermal stability, promising conditions for an improved CNO solar neutrino measurement. In spite of that, Borexino is expected to stop data taking within 2021.

Funding
This work is funded by the recruitment initiative of the Helmholtz Association of German Research Centres (HGF).

Data availability
Data used to obtain presented results are provided at https://bxopen.lngs.infn.it/, shared and owned by the Borexino collaboration.