COMPREHENSIVE GEONEUTRINO ANALYSIS WITH BOREXINO

This paper presents a comprehensive geoneutrino measurement using the Borexino detector, located at Laboratori Nazionali del Gran Sasso (LNGS) in Italy. The analysis is the result of 3262.74 days of data between December 2007 and April 2019. The paper describes improved analysis techniques and optimized data selection, 1 ar X iv :1 90 9. 02 25 7v 1 [ he pex ] 5 S ep 2 01 9 which includes enlarged fiducial volume and sophisticated cosmogenic veto. The reported exposure of (1.29 ± 0.05)×1032 protons × year represents an increase by a factor of two over a previous Borexino analysis reported in 2015. By observing 52.6+9.4 −8.6 (stat) +2.7 −2.1 (sys) geoneutrinos (68% interval) from 238U and 232Th, a geoneutrino signal of 47.0+8.4 −7.7 (stat) +2.4 −1.9 (sys) TNU with +18.3 −17.2% total precision was obtained. This result assumes the same Th/U mass ratio found in chondritic CI meteorites but compatible results were found when contributions from 238U and 232Th were both fit as free parameters. Antineutrino background from reactors is fit unconstrained and found compatible with the expectations. The null-hypothesis of observing a geoneutrino signal from the mantle is excluded at a 99.0% C.L. when exploiting detailed knowledge of the local crust near the experimental site. Measured mantle signal of 21.2+9.6 −9.0 (stat) +1.1 −0.9 (sys) TNU corresponds to the production of a radiogenic heat of 24.6 +11.1 −10.4 TW (68% interval) from 238U and 232Th in the mantle. Assuming 18% contribution of 40K in the mantle and 8.1+1.9 −1.4 TW of total radiogenic heat of the lithosphere, the Borexino estimate of the total radiogenic heat of the Earth is 38.2+13.6 −12.7 TW, which corresponds to the convective Urey ratio of 0.78+0.41 −0.28. These values are 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 fitting with a constraint on the number of expected reactor antineutrino events, the existence of a hypothetical georeactor at the center of the Earth having power greater than 2.4 TW at 95% C.L. is excluded. Particular attention is given to the description of all analysis details which should be of interest for the next generation geoneutrino measurements using liquid scintillator detectors. †Deceased in August 2019 1Present address: Dipartimento di Fisica, Sapienza Università di Roma e INFN, 00185 Roma, Italy 2Present address: Dipartimento di Fisica, Università degli Studi Federico II e INFN, 80126 Napoli, Italy 3Present address: Universidad Autónoma de Madrid, Ciudad Universitaria de Cantoblanco, 28049 Madrid, Spain 4Present address: Fermilab National Accelerato Laboratory (FNAL), Batavia, IL 60510, USA 5Present address: INFN Laboratori Nazionali del Gran Sasso, 67010 Assergi (AQ), Italy 6Also at: MTA-Wigner Research Centre for Physics, Department of Space Physics and Space Technology, Konkoly-Thege Miklós út 29-33, 1121 Budapest, Hungary Preprint submitted to Elsevier September 6, 2019

which includes enlarged fiducial volume and sophisticated cosmogenic veto. The reported exposure of (1.29 ± 0.05)×10 32 protons × year represents an increase by a factor of two over a previous Borexino analysis reported in 2015. By observing 52.6 +9. 4 −8.6 (stat) +2.7 −2.1 (sys) geoneutrinos (68% interval) from 238 U and 232 Th, a geoneutrino signal of 47.0 +8.4 −7.7 (stat) +2.4 −1.9 (sys) TNU with +18.3 −17.2 % total precision was obtained. This result assumes the same Th/U mass ratio found in chondritic CI meteorites but compatible results were found when contributions from 238 U and 232 Th were both fit as free parameters. Antineutrino background from reactors is fit unconstrained and found compatible with the expectations. The null-hypothesis of observing a geoneutrino signal from the mantle is excluded at a 99.0% C.L. when exploiting detailed knowledge of the local crust near the experimental site. Measured mantle signal of 21.2 +9. 6 −9.0 (stat) +1.1 −0.9 (sys) TNU corresponds to the production of a radiogenic heat of 24.6 +11.1 −10.4 TW (68% interval) from 238 U and 232 Th in the mantle. Assuming 18% contribution of 40 K in the mantle and 8.1 +1.9 −1.4 TW of total radiogenic heat of the lithosphere, the Borexino estimate of the total radiogenic heat of the Earth is 38.2 +13. 6 −12.7 TW, which corresponds to the convective Urey ratio of 0.78 +0. 41 −0. 28 . These values are 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 fitting with a constraint on the number of expected reactor antineutrino events, the existence of a hypothetical georeactor at the center of the Earth having power greater than 2.4 TW at 95% C.L. is excluded. Particular attention is given to the description of all analysis details which should be of interest for the next generation geoneutrino measurements using liquid scintillator detectors.

INTRODUCTION
Neutrinos, the most abundant massive particles in the universe, are produced by a multitude of different processes. They interact only by the weak and gravitational interactions, and so are able to penetrate enormous distances through matter without absorption or deflection. Thus, they represent a unique tool to probe otherwise inaccessible objects, such as distant stars, the Sun, as well as the interior of the Earth.
The present availability of large neutrino detectors has opened a new window to study the deep Earth's interior, complementary to more conventional direct methods used in seismology and geochemistry. For example, atmospheric neutrinos can be used as probe of the Earth's structure [1]. This absorption tomography is based on the fact that the Earth begins to become opaque to neutrinos with energies above ∼10 TeV. Thus, the attenuation of the neutrino flux, as measured by the signals in large Cherenkov detectors, provides information about the nucleon matter density of the Earth. Recently, IceCube determined the mass of the Earth and its core, its moment of inertia and verified that the core is denser than the mantle using data obtained from atmpospheric neutrinos [2]. A complementary information about the electron density could, in principle, be inferred by exploiting the flavor oscillations of atmospheric neutrinos in the energy range from MeV to GeV [3].
An independent method to study the matter composition deep within the Earth, can be provided by geoneutrinos, i.e. (anti)neutrinos emitted by the Earth's radioactive elements. Their detection allows to assess the Earth's heat budget, specifically the heat emitted in the radioactive decays.
It took several decades to prove these ideas feasible. Currently, large-volume liquid-scintillator neutrino experiments KamLAND [12,13,14,15] and Borexino [16,17,18] have demonstrated the capability to efficiently detect a geoneutrino signal. These detectors are thus offering a unique insight into 200 years long discussion about the origin of the Earth's internal heat sources.
The Borexino detector, located in hall-C of Laboratori Nazionali del Gran Sasso in Italy (LNGS), was originally designed to measure 7 Be solar neutrinos. However thanks to the unprecedented levels of radiopurity, Borexino has surpassed its original goal and has now measured all 8 the pp-chain neutrinos [20,21,22]. We report here a comprehensive geoneutrino measurement based on the Borexino data acquired during 3262.74 days (December 2007 to April 2019). Thanks to an improved analysis with optimized data selection cuts, an enlarged fiducial volume, and a sophisticated cosmogenic veto, the exposure of (1.29 ± 0.05)×10 32 protons × year represents a factor 2 increase with respect to the previous Borexino analysis [18].
A detailed description of all the steps in the analysis is reported, and should be important to new experiments measuring geoneutrinos, e.g. SNO+ [23], JUNO [24], and Jinping [25]. Hanohano [26] is an interesting, additional proposal to use a movable 5 kton detector resting on the ocean floor. As the oceanic crust is particularly thin and relatively depleted in HPEs, this experiment could provide the most direct information about the mantle. Finally, it is anticipated that using antineutrinos to study the Earth's interior will increase in the future based on the availability of new detectors and the continued development of analysis techniques. This paper is structured as follows: Section 2 introduces the fundamental insights on what the geoneutrino studies can bring to the comprehension of the Earth's inner structure and thermal budget. Section 3 details a description of the Borexino detector and the structure of its data. In Section 4, theν e detection reactionthe Inverse Beta Decay of the free proton, that will be abbreviated as IBD through the text -is illustrated. It is shown that only geoneutrinos above 1.8 MeV kinematic threshold can be detected, leaving 40 K and 235 U geoneutrinos completely unreachable with present-day detection techniques. Section 5 deals with the estimation of the expected antineutrino signal from geoneutrinos, through background from reactor and atmospheric neutrinos, up to a hypothetical natural georeactor in the deep Earth. Section 6 describes the non-antineutrino backgrounds, e.g. cosmogenic or natural radioactive nuclei whose decays could mimic IBD. The criteria to selectively identify the best candidates in the data, are discussed in Sec. 7, which involves the optimization of the signal-to-background ratio. Section 8 shows how the signal and background spectral shapes, expressed in the experimental energy estimator (normalized charge), were constructed and how the detection efficiency is calculated. Both procedures are based on Borexino Monte Carlo (MC) [27], that was tuned on independent calibration data. Section 9 introduces the analyzed data set and discusses the number of expected signal and background events passing the optimized cuts, based on Secs. 5 and 6. In Section 10, the Borexino sensitivity to extract geoneutrino signals is illustrated. Finally, Sec. 11 discusses our results. The golden IBD candidate sample is presented (Sec. 11.1) together with the spectral analysis (Sec. 11.2) and sources of systematic uncertainty (Sec. 11.3). The measured geoneutrino signal at LNGS is compared to the expectations of different geological models in Sec. 11.4. The extraction of the mantle signal using knowledge of the signal from the bulk lithosphere is discussed in Sec. 11.5. The consequences of the new geoneutrino measurement with respect to the Earth's radiogenic heat are discussed in Sec. 11.6. Placing limits on the power of a hypothetical natural georeactor, located at different positions inside the Earth, is discussed in Sec. 11.7. Final summary and conclusions are reported in Sec. 12.

WHY TO STUDY GEONEUTRINOS
Our Earth is unique among the terrestrial planets 9 of the solar system. It has the strongest magnetic field, the highest surface heat flow, the most intense tectonic activity, and it is the only one to have continents composed of a silicate crust [28]. Understanding the thermal, geodynamical, and geological evolution of our planet is one of the most fundamental questions in Earth Sciences [29].
The Earth was created in the process of accretion from undifferentiated material of solar nebula [30,31]. The bodies with a sufficient mass undergo the process 9 Mercury, Venus, Earth, and Mars The Earth has a concentrically layered structure with an equatorial radius of 6378 km. The metallic core includes an inner solid portion (1220 km radius) and an outer liquid portion which extends to a depth of 2895 km, where the core is isolated from the silicate mantle by the core-mantle boundary (CMB). Seismic tomography suggests a convection through the whole depth of the viscose mantle, that is driving the movement of the lithospheric tectonic plates. The lithosphere, subjected to brittle deformations, is composed of the crust and continental lithospheric mantle. The mantle transition zone, extending from a depth of 400 to 700 km, is affected by partial melting along the mid-oceanic ridges where the oceanic crust is formed. The continetal crust is more complex and thicker than the oceanic crust.
of differentiation, i.e. a transformation from a homogeneous object into a body with a layered structure. The geophysical picture of a concentrically layered internal structure of the present Earth ( Fig. 1), with mass M E = 5.97 ×10 24 kg, is relatively well established from its density profile, which is obtained by precise measurements of seismic waves on its surface. During the first differentiation, metallic segregation occurred and the core (∼32.3% of M E ) separated from the silicate Primitive Mantle or Bulk Silicate Earth (BSE). The latter further differentiated into the present mantle (∼67.2% of M E ) and crust (∼0.5% of M E ). 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 [32]. The outer core has an approximate 10% admixture of lighter elements and plays a key role in the geodynamo process generating the Earth's magnetic field. The core-mantle boundary (CMB) seismic discontinuity divides the core from the mantle. The mantle reaches temperature of about 3700 K at its bottom, while being solid but viscose on long time scales, so the mantle convection can occur. The latter drives the movement of tectonic plates at the speed of few cm per year.
A whole mantle convection is supported by high resolution seismic tomography [33], which proves existence of material exchange across the mantle, as in the zones of deeply subducted lithosperic slabs and mantle plumes rooted close to the CMB. At a depth of [400 -700] km, the mantle is characterized by a transition zone, where a weak seismic-velocity heterogeneity is measured. The upper portion of the mantle contains the viscose 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 underlying the CC included between the Moho discontinuity 10 and a seismic and electromagnetic transition at a typical depth of ∼175 km [34]. The CC with a thickness of (34 ± 4) km [35] has the most complex history being the most differentiated and heterogeneous layer. It consists of igneous, metamorphic, and sedimentary rocks. The OC with (8 ± 3) km thickness [35] is created along the mid-oceanic ridges, where the basaltic magma differentiates from the partially melting mantle up-welling towards the ocean floor.
Traditionally, direct methods to obtain information about the deep Earth's layers, from where there are few or no direct rock samples, are limited to seismology. Seismology provides relatively precise information about the density profile of the deep Earth [36], but it lacks direct information about the chemical composition and radiogenic heat production. Geoneutrinos come into play here: their small interaction cross-section (∼ 10 −42 cm 2 at MeV energy for the IBD, Sec. 4), on one hand, limits our ability to detect them, on the other hand, it makes them a unique probe of inaccessible innermost parts of the Earth. In the radioactive decays of HPEs, the amount of released geoneutrinos and radiogenic heat are in a well-known ratio (Eqs. 1 to 5). Thus, a direct measurement of the geoneutrino flux provides useful information about the TABLE I. Integrated terrestrial surface heat fluxes H tot estimated by different authors. The lower limit estimation [37] is due to the approach based only on direct heat flow measurements in contrast with the thermal model of half space cooling adopted by the remaining references.

Reference
Earth's heat flux [TW] Williams & Von Herzen (1980) [38] 43 Davies (1980) [39] 41 Sclater et al. (1980) [40] 42 Pollack et al. (1993) [41] 44 ± 1 Hofmeister et al (2005) [37] 31 ± 1 Jaupart et al. (2007) [42] 46 ± 3 Davies & Davies (2010) [43] 47 ± 2 composition of the Earth's interior [5]. Consequently, it also provides an insight into the radiogenic heat contribution to the measured Earth's surface heat flux. The heat flow from the Earth surface to space results from a large temperature gradient across the Earth [44]. Table I shows estimations of this heat flow, H tot , integrated over the whole Earth's surface. Different studies of this flux are based on several thousands of inhomogeneously distributed measurements of the thermal conductivity of rocks and the temperature gradients within deep bore holes. The existence of perturbations produced by volcanic activity and hydrothermal circulations, especially along the mid-ocean ridges (where the data are sparse), requires the application of energy-loss models [45]. Except for Ref. [37], the papers account for the hydrothermal circulation in the young oceanic crust by utilizing the half-space cooling model, which describes ocean depths and heat flow as a function of the oceanic lithosphere age. The latter is unequivocally correlated with the distance to mid-ocean ridges, where the oceanic crust is created [46] and from where the older crust is pushed away by a newly created crust. This approach leads to an H tot estimation between 41 and 47 TW, with the oceans releasing ∼70% of the total escaping heat. The most recent models [42,43] are in excellent agreement and provide a value of (46-47) TW with (2-3) TW error. However, Ref. [37], based only on direct measurements and not applying the half-space cooling model, provides a much lower H tot of (31 ± 1) TW. We assume a H tot = (47 ± 2) TW as the best current estimation.
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 of the Earth, i.e. cooling from the time of the Earth's formation when gravita-6 tional binding energy was released due to matter accretion H SC , and (ii) radiogenic heat H rad from HPEs' radioactive decays in the Earth. The relative contribution of radiogenic heat to the H tot is crucial in understanding the thermal conditions occurring during the formation of the Earth and the energy 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, as the following ratio [44]: where H CC rad is the radiogenic heat produced in the continental crust. The secular cooling of the core is expected in the range of [5 -11] TW [45], while no radiogenic heat is expected to be produced in the 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 [45]. Additionally, HPEs' abundances, and thus H rad of Eq. 6, are globally representative of BSE models, defining the original chemical composition of the Primitive Mantle. 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) isotopical abundances in the solar photosphere and in the CI chondrites ( Fig. 2 in [32]). Such correlations can be then assumed also for the material from which the Earth was created. The BSE models agree in the prediction of major elemental abundances (e.g. O, Si, Mg, Fe) within 10% [49]. Uranium and Thorium are refractory (condensate at high temperatures) and lithophile (where silicates over metal are preferred) elements. The relative abundances of the refractory lithophile elements are expected to be stable to volatile loss or core formation during the early stage of the Earth [57]. The content of refractory lithophile elements (e.g. U and Th), which are excluded from the core 11 , are assumed based on relative abundances in chondrites, and dramatically differ between different models. In Table II global masses of HPEs and their corresponding radiogenic heat are reported, covering a wide spectrum of BSE 11 Recent speculations [58] about possible partitioning of some lithophile elements (including U and Th) into the metallic core are still debated [59,60]. This would explain the anomalous Sm/Nd ratio observed in the silicate Earth and would represent an additional radiogenic heat source for the geodynamo process. compositional models.
Three classes of BSE models are adopted in this work: the Cosmochemical, Geochemical, and Geodynamical models, as defined in [32,54]. The Cosmochemical (CC) model [34] is characterized by a relatively low amount of U and Th producing a total H rad = (11 ± 2) TW. This model bases the Earth's composition on enstatite chondrites.
The Geochemical (GC) model class predicts intermediate HPEs'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 [49,56].
The Geodynamical (GD) model shows relatively high U and Th abundances. It is based on the energetics of mantle convection and the observed surface heat loss [53]. Additionally, an extreme model can be obtained following the approach described in [55], where the terrestrial heat H tot of 47 TW is assumed to be fully accounted for by radiogenic production H rad . When keeping the HPEs' abundance ratios fixed to chondritic values and rescaling the mass of each HPE component accordingly, one obtains estimates for Fully Radiogenic (FR) model (Table II).
A global assessment of the Th/U mass ratio of the Primitive Mantle could hinge on the early evolution of the Earth and its differentiation. The most precise estimate of the planetary Th/U mass ratio reference, having a direct application in geoneutrino analysis, has been refined to a value of M Th /M U = (3.876 ± 0.016) [61]. Recent studies [62], based on measured molar 232 Th/ 238 U values and their time integrated Pb isotopic values, are in agreement estimating M Th /M U = 3.90 +0. 13 −0.08 . Significant deviations from this average value can be found locally, especially in the heterogeneous continental crust. This fact is attributable to many different lithotypes, which can be found surrounding the individual geoneutrino detectors [63,64]. In the local reference model for the area surrounding the Borexino detector (see also Sec. 5.2), the reservoirs of the sedimentary cover, which account for 30% of the geoneutrino signal from the regional crust, are characterized by a Th/U mass ratio ranging from ∼0.8 (carbonatic rocks) to ∼3.7 (terrigenous sediments) [65].
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, often based on indirect arguments. In this puzzle, direct U and Th geoneutrino  [53]. The Cosmochemical (CC), Geochemical (GC), and Geodynamical (GD) BSE models correspond to the estimates reported in [54]; the Fully Radiogenic (FR) model is defined adopting the approach of [55], assuming that the total heat H tot = (47 ± 2) TW is due to only the radiogenic heat production H rad (U+Th+K). The CC and GD models correspond to the estimates based on predictions made by Javoy et al., 2010 [34] and Turcotte & Schubert, 2002 [53], respectively. The GC is based on estimates reported by McDonough & Sun, 1995 [49] with K abundances corrected following [56]. The radiogenic heat H rad released in the radioactive decays of HPEs is calculated adopting the element specific heat generation h (H rad = h × M with h(U) = 98.5 µW/kg, h(Th) = 26.3 µW/kg, and h(K) = 3.33 × 10 −3 µW/kg taken from [55]). It is assumed that the uncertainties on U, Th, and K abundances are fully correlated. measurements are candidates to play a starring role. 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 litospheric contribution, as we show in Sec. 11.5. The lithospheric contribution would be particularly small and easily determined on a thin, HPEs depleted oceanic crust. This would make the ocean floor an ideal environment for geoneutrino detection. Geoneutrino measurements can also contribute to the discussion about possible additional heat sources, which have been proposed by some authors. For example, stringent limits (Sec. 11.7) can be set on the power of a hypothetical Uranium natural georeactor suggested in [66,67,68,69] and discussed in Sec. 5.5. In future, by combining measurements from several experiments placed in distant locations and in distinct geological environments, one could test whether the mantle is laterally homogeneous or not [54], as suggested, for example, by the Large Shear Velocity Provinces observed at the mantle base below Africa and Pacific ocean [70].
In future, detection of 40 K geoneutrinos might be possible [71,72]. This would be extremely important, since Potassium is the only semi-volatile HPE. Our planet seems to show ∼1/3 [49] to ∼1/8 [34] Potassium when compared to chondrites, making its expected bulk mass span of a factor ∼2 across different Earth's models. Two theories on the fate of the mysterious "missing K" include loss to space during accretion [49] or segregation into the core [73], but no experimental evidence has been able to confirm or rule out any of the hypotheses, yet. As a consequence, the different BSE class models predict a K/U ratio in the mantle in a relatively wide range from 9700 to 16000 [54]. According to these ratios, the Potassium radiogenic heat of the mantle varies in the range [2.6 -4.3] TW, which translates to an average contribution of 18% to the mantle radiogenic power. We will use this value in the evaluation of the total Earth radiogenic heat from the Borexino geoneutrino measurement (Sec. 11.6).

THE BOREXINO DETECTOR
Borexino is an ultra-pure liquid scintillator detector [74] operating in real-time mode. It is located in the hall-C of the Gran Sasso National Laboratory in central Italy at a depth of some 3800 m.w.e. The rock above the detector provides shielding against cosmogenic backgrounds such that the muon flux is decreased to ( 0.003) · 10 −4 m −2 s −1 [75]. The general scheme of the Borexino detector is shown in Fig. 2. The detector has a concentric multi-layer structure. The outer layer (Outer Detector (OD)) serves as a passive shield against external radiation as well as an active Cherenkov veto of cosmogenic muons. It consists of a steel Water Tank (WT) of 9 m base radius and 16.9 m height filled with approximately 1 kt of ultra-pure water. Cherenkov light in the water is registered in 208 8" photo-multiplier tubes (PMTs) placed on the floor and outer surface of a Stainless Steel Sphere (SSS, 6.85 m radius), which is contained within the WT. The Inner Detector (ID) within the SSS comprises three layers and it is equipped with 2212 8" PMTs mounted on the inner surface of the SSS. Over time, the number of working PMTs in the ID has decreased, from 1931 in December 2007 to 1183 by the end of April 2019. The three ID layers are formed by the insertion of the two 125 µm thick nylon "balloons", the Inner Vessel (IV) and the Outer Vessel (OV) with the radii 4.25 m and 5.50 m, respectively. The two layers between the SSS and the IV, separated by the OV, form the Outer Buffer (OB) and the Inner Buffer (IB).
The antineutrino target is an organic liquid scintillator (LS) confined by the IV. The scintillator is composed of pseudocumene (PC, 1,2,4-trimethylbenzene, C 6 H 3 (CH 3 ) 3 ) solvent doped with a fluorescent dye PPO (2,5-diphenyloxazole, C 15 H 11 NO) in concentration of 1.5 g/l. The scintillator density is (0.878 ± 0.004) g cm −3 , where the error considers the changes due to the temperature instabilities over the whole data acquisition period. The nominal total mass of the target is 278 ton and the proton density is (6.007 ± 0.001) × 10 28 per 1 ton. A careful selection of detector materials, accurate assembling, and a complex radio-purification of the liquid scintillator guaranteed extremely low contamination levels of 238 U and 232 Th. After the additional LS purification in 2011, they achieved < 9.4 × 10 −20 g/g (95% C.L.) and < 5.7 × 10 −19 g/g (95% C.L.), respectively.
The buffer liquid, consisting of a solution of the dimethylphthalate (DMP, C 6 H 4 (COOCH 3 ) 2 ) light quencher in PC, shields the core of the detector against external γs and neutron radiation. The OV and IV themselves block the inward transfer of Radon emanated from the internal PMTs and SSS. The quencher concentration has been varied twice, changing it from the initial 5 g/l to 3 g/l and then to 2 g/l. These operations have reduced the density difference between the buffer and the scintillator in order to minimize the scintillator leak (appeared in April 2008) from the central volume through the small hole in the IV to the IB as much as possible. This campaign was mostly successful, but the IV shape has become non-spherical and changing in time. We are able to reconstruct the IV shape from the data itself, as it will be described in Sec. 3.3.
Borexino has a main data acquisition system (the main DAQ) and a semi-independent Fast Wave Form Digitizing (FWFD) or Flash Analog-to-Digital Converter (FADC) sub-system designed for energies above 1 MeV. Both systems process signals from both the ID and the OD PMTs, but in different ways. Every ID PMT is AC-coupled to an electronic chain made by an analogue front end (so-called FE boards, FEBs) followed by a digital circuit (so-called Laben boards 12 ). While the main DAQ treats every PMT individually, the FADC sub-system receives as input the sums of up to 24 analogue FEB outputs. More details about the Borexino data structure are given in Sec. 3.1.
The effective light yield in Borexino is approximately 500 detected photoelectrons (p.e.) per 1 MeV of deposited energy.
This results in the 5%/ √ E (MeV) energy resolution. Borexino is a position sensitive detector. For point-like events, the vertex is reconstructed based on the time-of-flight technique [20] with ∼10 cm at 1 MeV resolution at the center of the detector. For other positions with larger radii, the resolution decreases on average by a few centimeters.
A comprehensive calibration campaign [76] was performed in 2009.
It served as a base for understanding of the detector's performance and for tuning a custom, Geant4 (release 4.10.5)-based, Monte Carlo (MC) code called G4Bx2. This MC simulates all processes after the interaction of a particle in the detector, including all known characteristics of the apparatus [27]. Since the Borexino MC chain results in data files with the same format as real data, the same software can be applied to both of them. During the calibration, radioactive sources were employed in approximately 250 points through the IV scintillator volume. Using seven CCD cameras mounted inside the detector, the positions of the sources could be determined with a precision better than 2 cm. Several gamma sources with energies between 0.12 to 1.46 MeV were used for studying the energy scale. 222 Rn source, emitting alpha particles characterized by point-like interactions, was applied to study the homogeneity of the detector's response as well as position reconstruction.
For geoneutrino studies, employment of the 241 Am-9 Be source is of particular interest, since the emitted neutrons closely represent the delayed signal of an Inverse Beta Decay (Sec. 4), which is the interaction used to detect geoneutrinos. In addition, a 228 Th source emitting 2.615 MeV gammas was placed in 9 detector inlets, constructed in a way that the sources were practically positioned at the SSS. This calibration was fundamental in the optimization of the biasing technique used in the MC simulation of the external background.
Along with the special calibration campaign at the beginning of data acquisition, there are constant offline checks of the detector's stability and regular online PMTs' calibration. The time equalization among PMTs is performed once a week with a special laser (λ = 394 nm, 50 ps wide peak) calibration run. This procedure is of utmost importance for position and muon track reconstructions, as well as for the α/β discrimination (Sec. 3.4), based on their different fluorescence time profiles. The charge calibration of the single photoelectron response of each PMT is preformed typically 4 times a day. It is based on plentiful β − decays of 14 C (Q = 156 keV), inevitably present in each organic liquid scintillator and dominating the triggering rate, which varied between 20-30 s −1 (above roughly 50 keV threshold) during the analyzed period. A new Borexino Trigger Board (BTB) was installed in May 2016, in place of the old module which began to have failures. The thorough tests have proven unbiased performance and further improved stability of the detector.

Borexino data structure
Borexino electronics must handle about 10 6 events per day, which are dominated by 14 C decays. The residual flux of cosmogenic muons results in the detection of approximately 4300 per day internal muons which cross IV scintillator and/or the buffer, and approximately the same number of external muons which cross only the OD but not the ID [77]. The details of the read-out system, electronics, and trigger can be found in [74].
The key features relevant for the geoneutrino analysis are presented in this Section.

The main DAQ
The main DAQ reads individually all channels from both the ID and the OD when the BTB issues a global trigger. A global trigger condition occurs when at least one of the two sub-detectors has a trigger.
The ID trigger threshold, set to 25-20 PMTs triggered in a selected time window (typically 90 ns wide), corresponds to a deposited energy of approximately 50 keV. The trigger threshold in the Muon Trigger Board (MTB) is 6 hits in a 150 ns time window. The information of the OD trigger is stored in the 2 2 bit of the trigger word, and is referenced as the BTB4 condition, or Muon Trigger Flag (MTF). This is described in Sec. 3.2. In addition, the BTB processes special calibration triggers such as random, electronic pulse, and timing-laser triggers. These triggers are used to monitor the detector status and are regularly generated, typically every few seconds. Service interrupts, used in the generation of calibration triggers and synchronized with 20 MHz base clock, can raise other bit fields of the trigger word. Each event is associated with an absolute time read from a GPS receiver. Listed below are different Trigger Types (T T ), their names and main purpose. Each trigger (or event) has a 16 µs long DAQ gate, with the exception of a 1.6 ms long TT128 associated with the detection of cosmogenic neutrons.
• TT1 & BTB0 -the main trigger type for point-like events inside the ID, when OD did not see a signal. No bit of the trigger word is raised. After each event, there is approximately a 2-3 µs dead time window, during which no trigger can be issued.
• TT1 & BTB4 -the main category of internal muons which did trigger both the OD as well as the ID.
• TT128 -special 1.6 ms long trigger issued shortly (order 100 ns) after every TT1 & BTB4 event to guarantee a high detection efficiency of cosmogenic neutrons, sometimes created with very high multiplicity. The duration of this trigger type corresponds to about six times the neutron capture time [77].
• TT2 & BTB4 -the main category of external muons detected by the OD only.
• TT8 -calibration trigger of the timing laser used to monitor the quality of the laser pulse.
• TT32 -calibration trigger with the electronics pulse, used to monitor the number of working electronic channels.
• TT64 -the random trigger, used to monitor the dark noise and the 14 C-dominated energy spectrum below the BTB threshold.
During each trigger, the raw hits are the actual hits recorded during the event, while the decoded hits are all valid raw hits. A cluster is defined as an aggregation of decoded hits in the DAQ gate, well above the random dark-noise coincidence (typically, about 1 dark noise hit per 1 µs in the whole detector is observed). Each cluster represents a physical event and its visible energy is parametrized by the energy estimators, each normalized to 2000 working channels: • N P -number of triggered PMTs; • N h -number of hits within the cluster, including possible multiple hits from the same PMT; • N pe -number of photoelectrons, calculated as a sum of charges of all individual hits contributing to N h .
The event structure, i.e. the depiction of the time distribution of the decoded hits in the DAQ gate, is shown in Fig. 3 for different trigger types and for both the ID and the OD, as an integral of many events acquired during a typical 6 hours run. Instead, Fig. 4 shows in an analogous way the typical structure of individual events. The structure of trigger TT128 is shown in Fig. 5.

The FADC DAQ sub-system
There is an auxiliary data acquisition system based on 34 Fast Wave Form Digitizing (FWFD or FADC), 400 MHz, 8 bit VME boards with 3 input channels on each board. This additional read-out system was created to extend the Borexino energy range to ∼50 MeV, which is important to detect supernova neutrinos. The FADC DAQ energy threshold is ∼1 MeV. The system has been essentially continuously operational since it was started in December 2009, with the exception in 2014 when the system had technical problems and was not operating properly for about 5 months.
Each FADC channel receives a summed signal from up to 24 PMTs. The system works independently from the main DAQ if the PMTs and FEBs are operating. The FADC DAQ has a separate trigger, implemented through the programmable FPGA unit. The FADC event time window (DAQ gate) is 1.28 µs long. The trigger module receives additional trigger flags such as Muon Trigger Flag (MTF, will be explained in Sec. 3.2), the output of the OD analog sum discriminator, the calibration pulser, and laser flags. The trigger unit produces the permissions and prohibitions of signals allowing recording of physical events and moderating the rate of the calibration signals. Figure 6 shows examples of typical waveforms, in particular for a point-like event, a muon, and a calibration pulser signal. As it will be discussed in Sec. 3.2, the FADC system allows for an accurate pulse shape discrimination which is of paramount importance to achieve high muon detection efficiency.
The main and FADC DAQ systems are synchronized and merged offline, based on the GPS time of each trigger, using a special software utility. Typically, FADC waveforms are extended up to 16 µs (or 1.6 ms for T T 128 events) with a simple unperturbed baseline. When multiple FADC events correspond to different clusters of the same main DAQ event, they are merged together and eventual gaps are substituted with simple baselines.

Muon detection
The overall muon detection efficiency with the main DAQ system has been evaluated on 2008 -2009 data to be at least 99.992% [77]. When using the additional FADC system for muon detection, this efficiency increases to 99.9969%. The high performance of muon tagging over the 10 years period was demonstrated in a recent study of the seasonal modulation of the muon signal [75]. In this Section we review the muon tagging Internal muon (TT1&BTB4) methods in Borexino, and also provide updated analysis of the overall muon tagging efficiency using the main DAQ and evaluate its stability over the analyzed period. The principal muon tagging is performed by a specifically designed Water-Cherenkov OD. Additionally, ID pulse-shape and several special muon flags have been designed particularly for the geoneutrino analysis, where undetected single muons could become an important background among the approximate 15 antineutrino candidates per year. Below is a summary of different categories of muon detection in Borexino.

Muon Trigger Flag (MTF)
MTF muons trigger the OD and set the bit BTB4 of the trigger board, as described above in Sec. 3.1.

Muon Cluster Flag (MCF)
The MCF flag is set by a software reconstruction algorithm using the hits acquired from the OD. It considers separately two subsets of OD PMTs: those mounted on the SSS and those on the floor of the Water Tank. The MCF condition is met if 4 PMTs of either subset are fired within 150 ns.

Inner Detector Flag (IDF)
The IDF identifies muons based on different time profiles of hits originating from muon tracks with respect to those from point-like events. These time profiles are characterized by the peak time and the mean time of the cluster of hits from the ID. The mean time is the mean value of the times of decoded hits that belong to the same cluster, while the peak time is the time at which most of the decoded hits are deposited. The IDF is op- timized in three different energy ranges. For muons, the N h energy estimator is proportional to the track length across the ID rather than to the muon energy. Events with N h > 2100 and with mean time greater than 100 ns are considered as muons. In the lower energy interval 100 < N h < 2100, events with peak time >30(40) ns are tagged as muons above (below) N h = 900, respectively. In addition, the Gatti α/β discrimination parameter (see Sec. 7.4) is used in IDF to reject electronics noise, re-triggering of muons, as well as scintillation pulses coming from the buffer. In the very low energy interval N h < 80 dominated by 14 C background, IDF is not applied due to the limited performance of pulseshape identification. A detailed explanation of the IDF flag can be found in [77].

Strict Internal Muon Flag
The three above mentioned muon flags (MTF, MCF, IDF) are optimized in order to maximize the muon tagging efficiency while keeping the muon sample as clean as possible. Events that deposited energy in the ID, i.e. those that have a cluster of hits in the ID data, and are tagged by any of these three tags, are identified as Strict Internal Muons. These events are largely dominated by T T 1 events. However, we include in this category also the rare T T 2 events, which did not trigger the ID, but have a cluster with more than 80 hits in the ID-data. The mutual efficiencies of the three strict muon flags for muons with N h > 80 were studied over the years, as shown in Fig. 7. These are crucial in order to estimate the number of untagged muons which affect the geoneutrino candidate sample (Sec. 9.3).
The MTF efficiency has been mostly stable over the years. The MCF efficiency was slowly increasing since 2008 and has reached its optimal stable performance in 2011. The lowered MCF efficiency during the first years was due to occasional instability of the OD in the main DAQ data stream, that however did not influence neither the trigger system nor the MTF flag. In this case the OD data was not acquired. The IDF efficiency has been slowly decreasing since 2014 due to the decreasing number of active PMTs in the ID. This efficiency decrease is limited only to low energy ranges.
The average mutual efficiencies over the whole analyzed period are shown in Table III. The highest mutual inefficiency is 1.19% (the inefficiency of the IDF flag with respect to the MCF flag) and the highest mutual efficiency is 99.89% (the efficiency of the MTF flag with respect to the MCF flag). This can be used to calculate the overall inefficiency of the three muon flags as 0.0119

Special Muon Flags
In addition to the Strict Internal Muon Flags, there are also six kinds of Special Muon Flags. These are designed to tag the small, remaining fraction of muons at the cost of decreased purity of the muon sample. In fact, special tags mark as muons also noise events. In the antineutrino analysis, all special muons are conservatively treated as internal muons (see Sec. 6.1), regardless of the fact whether they have a cluster of hits in the IDdata.
• Special Flag 1 This flag was introduced to tag muons when the electronics was too saturated to correctly read out all photons, i.e. there was sufficiently high amount of raw hits (NR) but very few decoded hits (ND). Therefore, this flag tags events with NR > 200 when only less than 5% of them are decoded (ND/NR < 0.05). These events can have trigger type T T 1 or T T 2, without any condition regarding the BTB trigger word.
• Special Flag 2 This flag tags events of trigger type T T 1 or T T 2 with ND > 100 that have the 2 2 bit of the trigger word raised, but the OD triggered in coincidence with the service interrupt which generated service triggers. In this case, additional bits of the trigger word can be raised. For these events, only part of the muon data can be present in the DAQ gate or additional calibration pulses are possibly present in the event. • Special Flag 5 There is an extremely small number of T T 1 events, that by definition triggered the ID, but have zero clusters in the ID-data. Conservatively, if they are tagged by the MTF or MCF, they are considered as special muons.
• Special Flag 6 The events of trigger type T T 8, T T 32, or T T 64 (examples of these events are shown in the three lower rows of Fig. 4) that have unexpectedly high number of ND hits can be due to a muon occurring inside a service event. This is graphically shown in Fig. 8.

External Muon Flag
External muons are those that did not deposit energy in the ID and passed only through the OD. For all of them we require that they do not have a cluster in the ID-data.
Most of these muons are T T 2 events, very few also T T 1 triggers, that satisfy MTF or MCF (slightly modified) conditions.

FADC muon identification
The FADC DAQ improves the efficiency of muon detection in the ID. The advantage of the FADC system over the main electronics is the availability of detailed pulse shape information of the event. The selection of muons is obtained using a special algorithm which includes 8 different independent tests to classify a muon event. In addition to checking information about the triggers of and the data from the OD, four different classifiers are deployed. Three of the classifiers are based on machine learning and one of them considers formal and simple quantitative characteristics of the pulse shape, as the rise of the leading edge and the pulse amplitude. The following approaches are used for machine learning classifiers: a neural network based on the Multi-Layer Perceptron (MLP) [78], a Support Vector Machine (SVM) [79], and a Boosted Decision Tree (BDT) [80]. These were implemented to make a decision using the Toolkit for Multivariate Data Analysis (TMVA) with ROOT [78]. The classifiers use the event pulse-shape and additional parameters that determine the distribution of the digitized waveform. Different tests have different tagging efficiencies. There are three levels of reliability of tagging, defined according to the number of classifiers which tag an event as a muon. The lowest level of reliability is suitable for the antineutrino analysis. In this case, the maximum tagging efficiency is achieved with the price of the greatest over-efficiency. The combined muon tagging efficiency of the main and the FADC DAQ systems is 99.9969% [81].

Inner Vessel shape reconstruction
The shape of the thin nylon inner vessel (IV) holding the Borexino scintillator changes with time, deviating from a spherical shape. This deformation is a consequence of a small leak in the IV, with a location esti- mated as 26 • < θ < 37 • and 225 • < φ < 270 • [20]. The leak developed approximately in April 2008 and was detected based on a large amount of events reconstructed outside the IV. In order to minimize the buoyant force between the buffer and the scintillator liquids, their density difference was reduced by partial removal of DMP from the buffer by distillation, with negligible consequences on the buffer's optical behaviour. The evolution of the IV shape needs to be monitored and is also crucial for the antineutrino analysis.
In the geoneutrino analysis, the so-called Dynamical Fiducial Volume (DFV) (Sec.7.6) is defined along the time-dependent reconstructed IV shape. The reconstruction method, introduced in [20], is based on events in the 800 -900 keV energy range (N pe = 290 -350 p.e.) reconstructed on the IV surface (Fig. 9). These events originate in the radioactive contamination of the nylon and are dominated by 210 Bi, 40 K, and 208 Tl. The reconstructed position of selected events is fit assuming uniform azimuthal symmetry (xy plane) so that the θ-dependence of the vessel radius R can be determined. Three weeks of data provide sufficient statistics for this analysis.
A combination of a high-order polynomial, a Fourier series, and a Gaussian distribution was used as the function to fit the 2-dimensional distribution (R, θ). Figure 10(a) shows examples of the reconstructed IV shapes, in particular the reconstructed radius R as a function of θ. Fixed parameters of the fit are the two end points (θ = 0, θ = π) of the distribution where the vessel radius is imposed to be 4.25 m, since the IV is fixed to the end caps that are held in place by rigid support.
This procedure was cross-checked and calibrated over several ID pictures with an internal CCD camera system, which were taken throughout the data collection. The precision of this method is found to be ∼1% (±5 cm). An estimation of the active volume of scintillator is then calculated by revolving the (R, θ) function around the z-axis. Figure 10(b) shows the results of the rotational integration as a function of time. This evolution is very important to monitor the status and the stability of the Borexino inner vessel. The errors in the plot represent the goodness of the 2D fit only. The effect of the IV reconstruction precision is included in the systematic uncertainty of the geoneutrino measurement and will be discussed in Sec. 11.3.

α / β discrimination
The time distribution of the photons emitted by the scintillator depends on the details of the energy loss, and consequently on the particle type that produced the scintillation. For example, α particles have high specific energy loss due to their higher charge and mass. The energy deposition of a particle provides a way to characterize the pulse shape which can be used for particle identification [20]. Thus Borexino has different pulse-shape discrimination parameters which aid in the distinction of α-like and β-like interactions, and even more generally, to discriminate highly ionizing particles (α, proton) from particles with lower specific ionization (β − , β + , γ). These parameters were tuned using the Radoncorrelated 214 Bi(β − ) -214 Po(α) coincidence sample that was present in the detector during the Water-Extraction (WE) cycles. This occurred between June 2010 and August 2011 as a part of the scintillator purification process. The α/β discrimination parameters are important in the geoneutrino analysis since they help in distinguishing the nature of the delayed signals, as it will be shown in Sec. 7.4.

Gatti Optimal Filter
The Gatti optimal filter (G) is a linear discrimination technique, which allows to separate two classes of events with different time distributions [82,20]. First, using the typical β and α time profiles after the time-of-flight subtraction, the so-called weights w(t n ) are defined for time bins t n : where P α (t n ) and P β (t n ) are the probabilities that a photoelectron is detected at the time bin t n for α and β events, respectively. The Gatti parameter G for an event with the hit time profile f (t n ), after the time-of-flight subtraction, is then defined as: Figure 11(a) shows the distributions of the Gatti parameter for β particles from 214 Bi (G < 0) and α particles from 214 Po (G > 0).

Multi-Layer Perceptron
The Multi-Layer Perceptron (MLP) is a non-linear technique developed using deep learning for supervising binary classifiers, i.e. functions that can decide whether an input (represented by a vector of numbers) belongs to one class or another. In Borexino this technique was applied for α/β discrimination [83], and uses several pulseshape variables, parametrizing the event hit-time profile, as input. Among these variables are, for example, tailto-total ratio for different time bins t n , mean time of the hits in the cluster, their variance, skewness, kurtosis, and so on. The α-like events tend to have an MLP value of 0, while the β-like events tend to have an MLP value of 1, as it can be seen in Fig. 11(b).

ANTINEUTRINO DETECTION
Antineutrinos are detected in liquid scintillator detectors through the Inverse Beta Decay (IBD) reaction ( Fig. 12):ν e + p → e + + n, (9) in which the free protons in hydrogen nuclei, that are copiously present in hydrocarbon (C n H 2n ) molecules of organic liquid scintillators, act as target. IBD is a charge-current interaction which proceeds only for electron flavoured antineutrinos. Since the produced neutron is heavier than the target proton, the IBD interaction has a kinematic threshold of 1.806 MeV. The cross section of the IBD interaction can be calculated precisely with an uncertainty of 0.4% [84].
In this process, a positron and a neutron are emitted as reaction products. The positron promptly comes to rest and annihilates emitting two 511 keV γ-rays, yielding a prompt signal, with a visible energy E p , which is directly correlated with the incident antineutrino energy Eν e : The offset results mostly from the difference between the 1.806 MeV, absorbed from Eν e in order to make the IBD kinematically possible, and the 1.022 MeV energy released during the positron annihilation. The emitted neutron initially retains the information about theν e direction. However, the neutron is detected only indirectly, after it is thermalized and captured, mostly on a proton. Such a capture leads to an emission of a 2.22 MeV γ ray, which interacts typically through several Compton scatterings. These scattered Compton electrons then produce scintillation light that is detected in a single coincident delayed signal. In Borexino, the neutron capture time was measured with the 241 Am-9 Be calibration source to be (254.5 ± 1.8) µs [77]. During this time, the directional memory is lost in many scattering collisions. Figure 13 shows the N pe spectrum of delayed signals due to the gammas from captures of neutrons emitted from the 241 Am-9 Be calibration source placed in the center of the detector. In addition to the main 2.22 MeV peak due to the neutron captures on protons, higher energy peaks are clearly visible. The 4.95 MeV γs originate from the neutron captures on 12 C present in LS which occurs with about 1.1% probability. The higher energy peaks are from neutron captures on stainless steel nuclei (Fe, Ni, Cr) used in the source construction.
The pairs of time and spatial coincidences between the prompt and the delayed signals offer a clean signa- ture ofν e interactions, which strongly suppresses backgrounds. In the following sections, we refer to these signals as prompt and delayed, respectively.

EXPECTED ANTINEUTRINO SIGNAL
This section describes the expected antineutrino signals at the LNGS location ( r = 42.4540 • N, 13.5755 • E). We express them in Terrestrial Neutrino Units (TNU). This unit eases the conversion of antineutrino fluxes to the number of expected events: 1 TNU corresponds to 1 antineutrino event detected via IBD (Sec. 4) over 1 year by a detector with 100% detection efficiency containing 10 32 free target protons (roughly corresponds to 1 kton of LS).
Since we detect only electron flavour of the total antineutrino flux (Sec. 4), neutrino oscillations affect the expected signal expressed in TNU. Thus, the neutrino oscillations and the adopted parameters are discussed in Sec. 5.1. The evaluation of the expected geoneutrino signal from the Earth's crust and mantle is described in Sec. 5.2. Section 5.3 details the estimation of the signal from antineutrinos from the world reactors, the most important background for geoneutrino measurements. Atmospheric neutrinos, discussed in Sec. 5.4, also represent a potential background source for geoneutrinos. The existence of a georeactor, a naturally occurring Uranium fission in the deep Earth, was suggested by some authors. We present this idea as well as the expected signal from such a hypothetical source in Sec. 5.5. The final number of the expected events from each of these sources, expected in our data set (Sec. 3.1) and passing all optimized selection cuts (Sec. 7.8) will be presented in Sec. 9.2.

Neutrino oscillations
The presently accepted Standard Model of elementary particles describes neutrinos existing in three flavours (electron, muon, and tau) with masses smaller than 1/2 of the Z 0 boson mass.
The experiments with solar, atmospheric, as well as reactor antineutrinos observed that the neutrino flavour can change during the travel between the source and the detector. The process of neutrino oscillations has been established and confirmed that neutrinos have a non-zero rest mass. At present, most experimental results on neutrino flavor oscillation agree with a three neutrino scenario, where the weak neutrino eigenstates, i.e. flavor eigenstates (ν e , ν µ , ν τ ) mix with the mass eigenstates (ν 1 , ν 2 , ν 3 ) via the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix, parametrized with the three mixing angles (θ 12 , θ 13 , θ 23 ) and possible CP-violating and Majorana phases.
Therefore, to establish the expected electron antineutrino flux at a given site, it is necessary to consider the survival probability P ee of the electron flavoured neutrinos, which depends on the PMNS mixing matrix, as well as on the differences between the squared masses of the mass eigenstates, neutrino energy Eν, and the travelled baseline L. In the calculation of P ee for MeV antineutrinos, Eq. (37) from [85] is adopted, using the neutrino oscillation parameters as in Table IV, obtained by NU-FIT 3.2 (2018) [86] from a global fit to data provided by different experiments: where and L and Eν e are expressed in natural units ( = c = 1). We assume Normal Hierarchy for neutrino mass eigenstates (m 1 < m 2 < m 3 ) and neutrino oscillations in vacuum. In addition we assume ∆m 2 = ∆m 2 31 since |∆m 2 31 | ≈ |∆m 2 32 | and in Normal Hierarchy both difference squared masses are positive.
The size of matter effects on the P ee , when the neutrinos cross the Earth, is discussed individually for each antineutrino source in the following subsections. It depends on the baseline lenght in matter, the matter electron density N e , as well as on the antineutrino energy. Equation 62 of [85] is adopted in this calculation: P matter ee (L, Eν e , N e ) =c 4 13 (1 − 4s 2 12c 2 12 sin 2δ ) + s 4 13 , (13) where "tilde" denotes the mixing parameters (θ 12 ,δ) in matter, related with the vacuum oscillation parameters (θ 12 , δ) through relations: sin 2θ 12 = sin 2θ 12 (1 − µ 12 cos 2θ 12 ) δ =δ(1 + µ 12 cos 2θ 12 ) with where G F if the Fermi coupling constant.

Geoneutrinos
The Earth is a planet shining essentially in a flux of antineutrinos with a luminosity L ∼ 10 25 s −1 . For a detector placed on the continental crust, the expected U and 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 (7.41 · 10 7 kg −1 s −1 for 238 U and 1.62 · 10 7 kg −1 s −1 for 232 Th). Eν is geoneutrino energy. The geoneutrino energy spectra dn(i;Eν) dEνand , discussed in Sec. 5.2.1, are normalized 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 (Eq. 11). The matter effect (Sec. 5.1) is estimated to be of the order of 1% [87], i.e. much less than other uncertainties involved in the geoneutrino signal prediction. The average survival probability P ee = 0.55. The ρ( r ) is the density of the voxel emitting geoneutrinos and it is taken from geophysical models of lithosphere [35] and mantle [88]. 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 HPEs reservoirs, as discussed in Sec. 5.2.2.
To convert the differential geoneutrino flux dΦ(i;Eν, r) dEν to geoneutrino signal S (i) expressed in TNU (given in Sec. 5.2.2), it is necessary to account for the detection process via IBD on free protons and to perform integration over the geoneutrino energy spectra: where N p = 10 32 target protons, t is 1 year measuring time, and σ(Eν) is the IBD cross section [84]. Note that for a reference oscillated flux of 10 6 cm −2 s −1 , the geoneutrino signals from U and Th are S (U) = 12.8 TNU and S (Th) = 4.04 TNU, respectively. Considering the specific antineutrino production rates ε ν (U, Th), one can calculate the signal ratio R S for a homogeneous reservoir characterized by a fixed a(Th)/a(U) ratio: 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, we get R s = 0.27. Geophysical and geochemical observations of the lihtosphere constrain the a(Th)/a(U) = 4.3 (Table V), implying a signal ratio of 0.29 (Table VI). As a consequence, maintaining the global chondritic ratio of 3.9 for the bulk Earth, the inferred mantle ratio a(Th)/a(U) = 3.7, which correspond to a R S = 0.26.

Geoneutrino energy spectra
The expected geoneutrino signal depends on the shape and rates of the individual decays. The isotopes of HPEs, i.e. 238 U, 235 U, 232 Th, and 40 K release geoneutrinos with different energy spectra reported in Fig. 15(a)) for one decay of the head element of the chain. The number of emitted antineutrinos per decay is 6 for 238 U, 4 for 235 U and 232 Th, and 0.89 for 40 K. Note that the maximal energy of both 40 K and 235 U antineutrinos is below the IBD threshold (Sec. 4), while 0.38 and 0.15 antineutrinos per one decay are above this threshold for 238 U and 232 Th, respectively. The effective transitions producing detectable antineutrinos are given by 234m Pa and 214 Bi in 238 U decay chain and 228 Ac and 212 Bi in the 232 Th decay chain. Neglecting 210 Tl, having branching probability <0.1%, the 238 U and 232 Th antineutrino maximal energies are 3.27 MeV and 2.25 MeV, produced from 214 Bi and 212 Bi, respectively. In the energy distribution of U and Th antineutrinos reported in Fig. 15(a) and Fig. 15(b), the spectral structures of β and (β + γ) decays are clearly visible. Note that only 214 Bi decay spectral shape has been studied on the basis of experimental measurements [90]. The other energy distributions of antineutrinos are given with unknown uncertainties, since they are generally calculated assuming a well-known universal shape distribution. Figure 15(b)) shows the oscillated geoneutrino spectra expected at LNGS considering the geophysical and geochemical inputs as discussed in the following Sec. 5.2.2.

Geological inputs
The geoneutrino signal together with its uncertainty can be calculated considering the observational data concerning U and Th abundances in the lithosphere [35], the density profile of the Earth [88], and the BSE constraints on the global amounts of HPEs (Table II). In Table V the masses of the main reservoirs of the lithosphere are reported together with the HPEs' masses and the released radiogenic heat. Note that, although the mass of the bulk crust is less than 1% of the BSEs mass, it contains ∼35% of U and Th masses predicted by the GC model (Table II). The HPEs' radiogenic heat of the whole lithosphere is 8.1 +1.9 −1.4 TW. Following the scheme reported in Fig. 14, the expected geoneutrino signal in Borexino S (U+Th) can be expressed as the sum of three components: • S LOC (U+Th), the local crust (LOC) signal produced from the 6 • × 4 • 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 signal expected from the bulk lithosphere, as the sum of LOC and FFL contributions, is given in  [35] for the far-field lithosphere and from [65] for the local crust. The flux from the mantle is calculated assuming a two-layer distribution and adopting HPEs' abundances in BSE according to the GC model. The vertical dashed lines in both plots represent the kinematic threshold of the IBD interaction.  (Table II), in Table VII.

Local crust contribution
The S LOC (U+Th) is estimated adopting the local refined model based on specific geophysical and geochemical data described in [65]. The 492 km × 444 km region of continental crust surrounding the LNGS is divided in a Central Tile (CT) and the Rest of the Region (RR) (Fig. 14c). For the CT, which includes the crustal portion within ∼100 km from the Borexino detector, a 3D model with a typical resolution of (2.0 km × 2.0 km × 0.5 km) is built. The crustal structure of the CT is based on a simplified tectonic model that includes the main crustal thrusts and near vertical reflection seismic profiles of the CROP project [91]. The ∼35 km thick crust has a layered structure typical of Central Apennines, characterized by thick sedimentary cover (∼13 km) which is not reported in any global crustal model. It is constituted by three Permo-Mesozoic carbonatic successions and a unit of the Cenozoic terrigenous sediments. Since the local seismic sections do not highlight any evidence of middle crust, the crystalline basement is subdivided into upper crust (∼13 km) and lower crust (∼9 km). The U and Th mass abundances are obtained by ICP-MS and gamma spectroscopy measurements of the rock samples collected within 200 km from the LNGS and from representative outcrops of upper and lower crust of the south Alpine basement. It's relevant to note that ∼75% of the sedimentary cover volume of CT is constituted by Mesozoic carbonates particularly poor of U and Th. It implies that the overall U and Th abundances of sediments are a(U) = (0.8 ± 0.2) µg/g and a(Th) = (2.0 ± 0.5) µg/g to compare with a(U) = (1.73 ± 0.09) µg/g and a(Th) = (8.10 ± 0.59) µg/g [92] used for the global crustal estimations. A geophysical model with a lower spatial resolution (0.25 • × 0.25 • ) is built for the RR, which treats the sedimentary cover as a single and homogeneous layer with the same U and Th abundances of CT sediments. The geoneutrino signal of the LOC is S LOC (U+Th) = (9.2 ± 1.2) TNU 13 (Table VI) where 77% of the signal originates from U and Th distributed in the CT. The maximal and minimal excursions of various input values and uncertainties reported in [65] are taken as the ±3σ error range. The U and Th signal errors are conservatively considered fully positively correlated. Note that the reduction of ∼6 TNU with respect to the estimations of the global reference model [35] is due to presence of thick TABLE V. Total masses, HPEs' masses M, and total radiogenic heat H rad of the lithosphere and its components, the CLM and the crust (bulk crust = CC + CC). The error propagation assumes no correlation among different components and is performed using a Monte Carlo simulation to propagate the asymmetrical uncertainties of the non-Gaussian distributions. The median values and the 1σ uncertainties are shown. Due to the positive asymmetry of the distributions, in some cases the median of the output matrices is not coincident with the sum of the medians of the individual components.    sedimentary deposits composed primarily of U-and Th-poor carbonate rocks. The signal ratio R s (Eq. 18) for the local crustal contribution is 0.24 (Table VI).

Far field lithosphere contribution
The FFL includes the CLM and the remaining crust after subtracting the LOC (Fig. 14). The geoneutrino signal S FFL is calculated adopting the 1 • × 1 • geophysically based, 3D global reference model [35], which provides the abundances and distributions of HPEs in the lithosphere, together with their uncertainties.
The crust is subdivided in 64,800 cells labeled with their thickness, density, and velocity of compressional and shear waves for eight layers (ice, water, three sediment layers, upper, middle, and lower crust). The total crustal thickness and the associated uncertainty correspond, respectively, to the mean and the half range of three crustal models: • CRUST 2.0 [93,94], a global crustal model with (2 • ×2 • ) resolution based on refraction and reflec-tion seismic experiments and on extrapolations from geological and tectonic settings for regions lacking field measurements.
• CUB 2.0 [95], a (2 • × 2 • ) resolution model obtained with shear velocities reporting the uncertainties of crustal thickness. The model is obtained by inverting surface wave dispersion data performing a Monte Carlo multi-step process, using a priori constraints.
• GEMMA [96], a high-resolution (0.5 • × 0.5 • ) map of Moho depth obtained by inverting satellite gravity field data collected by GOCE. Additional external information (e.g. topography, bathymetry, and ice sheet models) and prior hypotheses on crustal density relative to the main geological provinces are taken also into account.
The relative thickness of the crustal layers are incorporated from CRUST 2.0, while the information about the sedimentary cover is adopted from [97]. The HPEs' abundances in the sediments, OC, and upper crust layers are taken from published reference values reporting the uncertainties [35], while U and Th abundances in the deep crust are inferred using seismic velocity arguments. The distinctive ultrasonic velocities reported in geophysical databases can be related to acidity (SiO 2 content) of igneous rocks, which is generally correlated with the U and Th abundances.
The CLM is geophysically and geochemically distinct from the rest of the mantle (the sublithospheric mantle), and it is characterized by abundances of a(U) = 0.03 +0.05 −0.02 µg/g and a(Th) = 0.15 +0.28 −0.10 µg/g taken from a database of ∼500 xenolithic peridotite samples, representing the typical rock types of the CLM. The CLM geoneutrino signal is The geoneutrino signal of the FFL is S FFL (U+Th) = 16.3 +4.8 −3.7 TNU (Table VI) and it constitutes the 63% of the signal of the bulk lithosphere.

Mantle contribution
Earth scientists have debated the picture of the mantle convection over the last decades. Some geochemical arguments support a two-layer convection, while geophysical reasoning affirms a whole-mantle convection. The main arguments for a layered mantle are based on (i) chemical and isotopic differences between Mid-Ocean-Ridge Basalts (rocks differentiated from mantle transition zone depleted in incompatible elements) and Ocean Island Basalts (rocks melted out from the deeper enriched mantle), (ii) isotope variations between continental and oceanic crust, and (iii) the missing radiogenic heat source paradox [56].
Beyond this controversy, all models agree that U and Th abundances are basically spherically distributed and non-decreasing with depth. This is an important point which permits us to keep the masses of HPEs in the unexplored mantle as free parameters and to provide constraints on the mantle contribution to the geoneutrino signal. There are three different predictions for the mantle signal given by: • Low scenario (LS) (Fig. 16a): after subtracting high lithospheric masses of U and Th (median + 1σ, Table V) from the BSE masses (Table II), the remaining masses are placed in a layer just above CMB, • Intermediate scenario (IS) (Fig. 16b) assumes HPEs' distribution in two layers , a lower enriched mantle (EM) and an upper depleted mantle (DM) separated at 2180 km of depth (for more details, see section 2.4 of [35]).
• High scenario (LS) (Fig. 16c): after subtracting low lithospheric masses of U and Th (median -1σ, Table V) from the BSE masses (Table II), the remainder is distributed homogeneously in the mantle.
Assuming these low and high scenarios, the predicted mantle geoneutrino signals S mantle , for each of the BSE models (Table II), are reported in Table VII, together with the corresponding HPEs' masses (M mantle ), the radiogenic heat power (H mantle rad ), and the convective Urey ratio (UR CV , Eq. 6). Since the BSE model which is based on enstatitic chondrites composition [55] is poor in HPEs, the Th mass calculated after lithosphere subtraction shows slightly negative values which have been set to zero.

Total geoneutrino signal at Borexino
The total geoneutrino signal from the bulk lithosphere expected at Borexino is a crucial piece of information for extracting the mantle signal from the Borexino measurement (Sec. 11.5). Since LOC and FFL are modelled independently, for each element the signal contributions are summed as linearly independent. The obtained bulk lithosphere signal is S LS (Th+U) = 25.9 +4.9 −4.1 TNU (Table VI).
Considering the intermediate scenario ( Fig. 16b) for each mantle model (Table VII), the total expected geoneutrino signal can cover a wide range from S CC tot = 28.5 +5.5 −4.8 TNU to S GD tot = 45.6 +5.6 −4.9 TNU passing through S GC tot = 34.6 +5.5 −4.8 TNU (Table VIII). The highest signal S FR tot = 55.3 +5.7 −5.0 TNU is given by a Fully genic Earth. The estimated 1σ error of the mantle signals in Table VIII corresponds to [S HS mantle − S LS mantle ]/6 and is conservatively summed to lithospheric uncertainty as fully positive correlated.
Plotting the cumulative geoneutrino signal, as a function of the distance from Borexino (Fig. 17), we observe that 40% of the total signal comes from U and Th in the regional crust that lies within 550 km of the detector. Up to a distance of ∼150 km from Borexino, 100% of the geoneutrino signal is generated from the LOC.
The geoneutrino spectrum expected at LNGS is presented in Fig. 15(a)). Figure 18 shows instead the geoneutrino spectrum as expected to be detected via the IBD interaction, showing explicitly the contributions from 238 U and 232 Th, as well as those from the bulk lithosphere and the mantle.  (Table II): the low and high mantle values are obtained by subtracting the 1σ high and low values of the lithosphere (Table V), respectively. The range of the expected geoneutrino signal is obtained according to the low scenario (minimal value) and the high scenario (maximal value). The convective Urey ratio UR CV values are obtained assuming a total heat flux H tot (U+Th+K) = 47 TW and taking into account the radiogenic heat produced by the continental crust, H CC rad = 6.8 +1.4 −1.1 TW (Table V).

Reactor antineutrinos
The main source of background in geoneutrino detection is the production of electron antineutrinos by nuclear power plants, the strongest man-made antineutrino source. Many nuclei, produced in the fission process of nuclear fuel decay through β-processes with the consequent emission of electron antineutrinos, the socalled reactor antineutrinos. Their energy spectrum extends up to 10 MeV, well beyond the end point of the geoneutrino spectrum (3.27 MeV). As a consequence, in the geoneutrino energy window (1.8 -3.27 MeV), there is an overlap between geoneutrino and reactor antineutrino signals.
At present, there are approximately 440 nuclear power reactors in the world, providing, nominally, a total amount of about 1200 GW thermal power, corresponding to approximately 400 GW of electrical power. With ∼200 MeV average energy released per fission and 6ν e produced along the β-decay chains of the neutron-rich unstable fission products, a reactor with a typical thermal power of 3 GW emits 10 20ν e s −1 . 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. The spectrum of the reactor antineutrino events expected to interact in a detector with efficiency ε, the number of target protons N p , and during the time t:  [35] for FFL and from [65] for the LOC. The signal from the mantle is calculated assuming a two-layer distribution and adopting HPEs' abundances in BSE according to the GC model.
where the index r cycles over N rea reactors considered: P r is its nominal thermal power, L r is the reactor-to-detector distance (the core positions are taken from [98] and we assume a spherical Earth with radius R = 6371 km), < LF > r indicates the weighted average of monthly thermal load factors (LF mr ). The index i stands for the different components of nuclear fuel ( 235 U, 238 U, 239 Pu, and 241 Pu), p i is the power fraction of the component i, Q i is the energy released per fission of the component i taken from [99] with a 0.2% quoted uncertainty, φ i (Eν e ) is the antineutrino spectrum originating from the fission of the i th component, σ(Eν e ) is the IBD cross section [84], and P ee is the survival probabili use in Eq. 11. The < LF > r is calculated for each reactor r: where E m are the monthly Borexino exposures in kton × year during the 137 months period from December 2007 to April 2019 (Sec. 9.1). In our calculation the nominal thermal power P r and the monthly load factors L mr originate from the Power Reactor Information System (PRIS), developed and maintained by the International Atomic Energy Agency (IAEA) [100]. Each year in summer, the PRIS produces documents containing information about the nuclear power reactor performance relative to the previous year 14 . The L mr reported are defined as the ratio between the net electrical energy produced during a reference period (after subtracting the electrical energy taken by auxiliary units) and the net electrical energy that would have been supplied to the grid if the unit were operated continuously at the nominal power during the whole reference period. In our calculation we assume that such electrical load factors are equal to the thermal ones, which are not available at present. We also consider an uncertainty on thermal power P r of the order of 2%.
Concerning power fractions p ri , one has to take into account that throughout the years, several technologies in building nuclear power plants have been developed. Different core types are characterized by different fuel compositions which give rise to different isotope contributions to total thermal power. In addition, during the power cycle of a nuclear reactor, the composition of the fuel changes since Pu isotopes are bred and U is consumed. Thus, the power fractions p ri are in principle quantities which depend on the reactor type and time. At the moment, we do not know the exact time-dependent fuel composition in each core operating in the world, so, as in [87], we assume some representative values for the power fractions.  [87] is considered as the uncertainty on p ri . The antineutrino energy spectra φ i (Eν e ) deserve particular attention. Experimental results from the reactor antineutrino experiments Daya Bay [101], Double CHooz [102], RENO [103], NEOS [104] coherently show that the measured IBD positron energy spectrum deviates significantly from the spectral predictions of Mueller et al. 2011 [105] in the energy range between 4 -6 MeV: the so-called 5 MeV excess.
In addition, an overall deficit is observed with respect to the prediction of [105]. The origin of this effect is being still debated [106,107,108]. In order to take into account this effect, we proceed in the following way: first, we calculate the neutrino spectra corresponding to all four isotopes according to the parametrization from [105]. Then we multiply the total spectrum by an energy-dependent correction factor based on the Daya Bay high precision measurement (extracted from lower panel of Fig. 3 in [101]). The shapes of the antineutrino spectra expected in Borexino in the period from December 2007 to April 2019, without and with the 5 MeV excess, are compared in Fig. 19. The effect of this shape difference on the precision of the geoneutrino measurement is discussed in Sec. 11.3. Finally, the expected signal from reactor antineutrinos S rea , expressed in TNU, can be obtained by integrating Eq. 19 and assuming 100% detection efficiency for a detector containing N p = 10 32 target protons and operating continuously for t = 1 year. The The quoted errors are at 1σ level, where the uncertainties related to reactor antineutrino production, propagation, and detection processes are estimated using a Monte Carlo-based approach discussed in [87]. As expected, by introducing the "5 MeV excess", the predicted signal varies by about 6% consistently with the normalization factor R = 0.946 found in [101]. In [87] the authors investigated the matter effects concerning the antineutrino propagation from the reactor to several experimental sites including LNGS. Since the size of the effect is proportional to the baseline travelled in matter, a maximum effect is found for the location in Hawaii, far away from all the reactors, and amounts to 0.7%. For Borexino, it can be considered negligible with respect to the overall uncertainties on the reactor antineutrino signal.

Atmospheric neutrinos
We have studied atmospheric neutrinos as a potential background source for the geoneutrino measurement. Atmospheric neutrinos originate in sequential decays of π(K) ± mesons and µ ± muons produced in cosmic rays' interactions with atmospheric nuclei. The flux of atmospheric neutrinos, the energy spectrum of which is shown in Fig. 20, contains both neutrinos and antineutrinos, and the muon flavour is roughly twice abundant than 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.
For energies above 100 MeV, the atmospheric neutrino fluxes are taken from the HKKM2014 model [109], while below 100 MeV the fluxes from the FLUKA code [110] are adopted. The resulted energy spectrum is the one shown in Fig. 20. We consider the neutrino fluxes averaged from all directions. We calculated the flavour oscillations during neutrino propagation through the Earth, including the matter effects, with the modified Prob3++ software [111] that comprises 1 km wide constant-density layers according to the PREM Earth's model [88]. The neutrino interactions with 12 C, 13 C, and 1 H nuclei were simulated with the GENIE Neutrino Monte Carlo Generator (version 3.0.0) [112]. GENIE output final state particles are used as primary particles for the G4Bx2 Borexino MC code [27]. The number of expected IBD-like interactions due to atmospheric neutrinos passing all the optimized selection cuts (Sec. 7.8) will be given in Sec. 9.2, Table XIV.

Georeactor
A possible existence of georeactor, i.e. natural nuclear fission reactor in the Earth interior, was first suggested by Herndon in 1993 [66]. Since then, several authors have discussed its possible existence and the characteristics. Different models suggest the existence of natural nuclear reactors at different depths: at the center of the core [67], at the inner core boundary [68], and the core-mantle boundary [69]. These models predict that the georeactor output power sufficient to explain terrestrial heat flow measurements (Table I) and helium isotope ratios in oceanic basalts [113].
In principle, characteristics of the antineutrino spectrum observed at the Earth's surface could specify the location and the power of a georeactor, discriminating among these models [114]. Some authors [115] even suggested an alternative interpretation of the KamLAND data through the existence of a ∼30 TW georeactor at the boundary of the liquid and solid phases of the Earth's core. Borexino 2013 results [17] set a 4.5 TW upper bound at 95% C.L. for a georeactor in the Earth's center. KamLAND has also studied this hypothesis and obtained an upper limit of 3.7 TW at 95% C.L. [15].
In this paper, we set the new upper bounds on the power of a potential georeactor in Sec. 11.7. In order to be able to set such limits for different hypothetical locations of the georeactor, we have calculated the expected antineutrino spectra of a 1 TW point-like georeactor, operating continuously during the data taking period (December 2007 -April 2019) with the power fractions of fuel components as suggested in [116] ( 235 U : 238 U ≈ 0.76 : 0.23). The energy released and the antineutrino spectra per fission are as in Sec. 5.3. Using only the flux parametrisation of [105], we neglect the question of 5 MeV excess. As shown in Fig. 21(a), we consider three different depths as extreme georeactor locations: 1) GR2: the Earth center (d = R Earth ), 2) GR1: the CMB placed 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 (Sec. 5.1), the matter effect (Eq. 13) is taken into account by assuming an Earth constant density ρ = 5.5 g/cm 3 . With respect to oscillations in vacuum, we observe 1.4% increase of the signal. The expected oscillated spectra for 1 TW georeactor at the three locations GR1, GR2, and GR3 are shown in Fig. 21(b). Since the oscillation length of MeV neutrinos is much smaller with respect to the studied baselines, we observe very fast oscillations in all three spectra. As it will be discussed in Sec. 8, Fig. 34, Borexino energy resolution does not allow to distinguish the spectral differences due to oscillations. The total expected signal S geo is reported in Table VIII. Concerning error estimation, the same procedure as in Sec. 5.3 in the determination of the contribution due to oscillation parameters, IBD cross section, and energy released per fission were adopted. Uncertainties due to power fractions were not taken into account, since the georeactor is assumed not to have a significant change of isotopic composition on a scale of few years. The systematic uncertainties due to our simplified treatment of neutrino oscillations in matter were estimated: the density variation in the range between 2.2 g cm −3 (crust) and 13 g cm −3 (core) causes about 3% (1%) variation of the expected signal in the total (geoneutrino) energy window. The errors reported in Table VIII have been calculated by adding in quadrature the errors discussed above.  VIII. Summary of the antineutrino signals expected at LNGS in different energy windows from geoneutrinos (for the whole Earth and from the mantle separately), reactor antineutrinos, and from 1 TW georeactor located in three different positions. The geoneutrino signals are calculated summing the lithospheric contribution S LS (Th + U) = 25.9 +4.9 −4.1 TNU and the mantle signal S mantle . The latter is predicted according to Cosmochemical (CC), Geochemical (GC), Geodynamical (GD), and Fully genic (FR) models (Table VII). The central value represents the intermediate scenario estimates (Fig. 16) 3.24 ± 0.10

NON-ANTINEUTRINO BACKGROUNDS
In this section we describe the origin of non-antineutrino backgrounds for geoneutrino measurements.
In particular, the cosmogenic background in Sec. 6.1, background due to accidental coincidences in Sec. 6.2, due to the (α, n) and (γ, n) reactions in Sec. 6.3 and Sec. 6.4, respectively, while the radon-correlated background in Sec. 6.5 and 212 Bi-212 Po coincidences in Sec. 6.6. Typically, the rates of these backgrounds depend on the selection cuts for IBD events, whose optimization is described in Sec. 7.
Therefore, the final evaluation of the non-antineutrino background levels is given in the following Sec. 9.

Cosmogenic background
One of the most important non-antineutrino backgrounds for geoneutrino detection is the cosmogenic background due to the residual muon flux. It is necessary to eliminate muons along with their spallation daughters as they can imitate the IBD signals. The various cosmogenic backgrounds relevant in geoneutrino measurement are explained in this section. The muon detection methods and efficiencies were already explained in Sec. 3.2.
The actual evaluation of the cosmogenic background passing the IBD selection cuts is explained in detail in Sec. 9.3.

Hadronic background
Table IX summarizes the lifetimes (τ), Q-values, branching ratios of the relevant decays, and the production rates of the three hadronic spallation products relevant for geoneutrino measurement. 9 Li and 8 He isotopes are the most important cosmogenic background, since they can produce (β − + n) pairs during their decays: 9 Li → e − +ν e + n + 2α 8 He → e − +ν e + n + 7 Li, indistinguishable from IBD's products. Generally speaking, 1 to 2 s veto after all internal muons can sufficiently suppress this kind of background. Due to the very small residual muon flux at LNGS, this kind of approach was applied in the previous Borexino geoneutrino analyses and is further improved in this work, as will be shown in Sec. 7.1. Borexino has studied the production of hadronic cosmogenic backgrounds [117]. For the 8 He production rate, only an upper limit was set (Table IX). When taking the 3σ CL limit as the production rate itself, the ratio of products [production rate × branching ratio] for 9 Li: 8 He is 6:1 in Borexino. Considering also the fact that in the far detector of Double Chooz, where the 9 Li production rate is 575 times higher than in Borexino, the 8 He production rate is observed to be compatible with zero [119], we conclude that the (β − + n) decaying hadronic background is largely dominated by 9 Li. In addition, two β − s from two decays of cosmogenic 12 B can imitate IBD signals: 12 B → e − +ν e + 12 C.
The second 12 B decay should happen within the IBD coincidence time window (1.28 ms, Sec. 7) in order to be a background for geoneutrino analysis. Since this time window is much smaller when compared to 12  In addition, the 13.4 MeV Q value of 12 B decay largely extends over the end-point of geoneutrino as well as reactor antineutrino spectra, which further suppresses the importance of this background for geoneutrinos.

Untagged muons
The small amount of muons that go undetected can cause muon-related background in the geoneutrino sample. This kind of background can consist of: • µ + µ: two untagged muons close in time. In particular muons with short tracks inside the buffer, where the scintillation light yield is strongly suppressed, could fall within the IBD selection cuts.
• µ + n: an untagged muon, again especially a buffer muon, can mimic a prompt, while muon-generated spallation neutron is indistinguishable from an IBD delayed signal. Also, the capture time of cosmogenic and IBD-produced neutrons is the same.
• Muon daughters: obviously, after undetected or unrecognized muons, no vetoes are applied. Thus, their spallation products can mimic IBDs. They can be, for example, neutron -neutron pairs or other hadronic backgrounds or fast neutrons described below.

Fast cosmogenic neutrons
Cosmic muons can produce neutrons that have a very hard energy spectrum extending to several GeV. These In principle, some level of pulseshape discrimination between the scattered proton and the positron signal, prompt signal of a real IBD, is possible.
In Borexino, 2 ms veto is sufficient to suppress this background after muons that cross the detector and are detected either as internal or external muons. Possible background of this kind can arise from undetected muons crossing the water tank and from the muons producing spallation products in interactions outside the detector, e.g. in the surrounding rock. The details of the estimation of both categories of this kind of background is explained in Sec. 9.3.

Accidental coincidences
In Borexino, the rate of prompt-and delayed-like events is 0.03 s −1 and 0.01 s −1 , respectively, in the whole scintillator volume. The accidental coincidences of these events in time and space can imitate IBD signals. These events are dominated by the external background, thus their reconstructed positions are mostly at large radii and close to the IV. This means that increasing the Fiducial Volume (FV) for the geoneutrino analysis necessarily increases also the rate of accidental coincidences. These coincidences are searched in an extended off-time window of 2 ms -20 s and scaled back to the geoneutrino search time window of 1.28 ms for a certain set of selection cuts. The evaluation of this background is explained in detail in Sec. 9.4.

(α, n) background
The decays along the chains of 238 U, 235 U, and 232 Th generate α particles that can initiate (α, n) reactions, a possible background for geoneutrinos. Thanks to the ultra-radiopurity levels achieved in Borexino, the only relevant source of α-particles in this terms is 210 Po (τ = 199.6 days), a product of 222 Rn, found fully out of equilibrium with the rest of the 238 U chain. 210 Po decays into a stable 206 Pb emitting α with the energy of 5.3 MeV: Thanks to the α/β discrimination techniques developed in Borexino (Sec. 3.4), it is possible to count 210 Po decays via an event-by-event basis, as it will be shown in Sec. 9.5, along with the estimation of this background passing the IBD selection cuts. The (α, n) reactions can occur on different nuclides, but it takes place mostly on 13 C in the Borexino scintillator: 13 According to the recent revision [120] of the previous data [121,122], the respective cross section equals to 200 mb. The produced neutron, with energies up to 7.3 MeV, is indistinguishable from the neutron that originates in the IBD interaction (Sec. 4). The relevance of this interaction for the antineutrino measurement, first studied by KamLAND [13], arises from the fact that there are also three possibilities for the generation of the prompt, as it is schematically shown in Fig. 22: • prompt I: γ-emission with energy of 6.13 MeV or 6.05 MeV, as a result of 16 O * deexcitation, if 16 O is produced in such an excited state; • prompt II: recoil proton appearing after scattering of the fast neutron on proton; • prompt III: 4.4 MeV γ-ray that is a product of the two-stage process: First, 12 C is excited into 12 C * in an inelastic scattering off fast neutron. Then, 12 C * transits to the ground state, accompanied by the γ emission: n + 12 C −→ 12 C * + n,

(γ, n) interactions and fission in PMTs
Energetic γ rays are produced by neutron capture reactions in the detector materials and in the surrounding rocks or in natural radioactive decays. These γ's may in turn give (γ, n) reactions with the nuclei of the LS or the buffer in Borexino. In particular, if the γ-ray makes a Compton scattering before being absorbed, its interaction and the following neutron capture gives a coincidence that almost perfectly mimics the IBD signal and in such a case, the pulse shape discrimination is not effective. Table X shows the (γ, n) reactions' thresholds for the most abundant isotopes in the Borexino scintillator.
The (γ, n) interaction having the lowest energy threshold (2.22 MeV) is the one on 2 H. Taking into account that the IBD's prompt energy starts at 1 MeV (Eq. 10) and considering the energy resolution, we conclude that only gammas with energies higher than 3 MeV could be a source of background in our analysis. Although there are some uncertainties in the branching ratios listed in the isotopes' tables, the natural U and Th chains do not emit sizeable γ lines at energies above 3 MeV. Thus, the intrinsic scintillator and vessel contaminations are not an important source of background, while muon and neutron induced gammas have to be carefully studied. Natural radioactivity of detector materials can be a source of an additional type of background. Spontaneous fission can generate neutrons up to several MeV. These can mimic antineutrino signals in the target mass. 238 U has by far the largest fission probability among the nuclei we can consider and equals to 5.45 · 10 −7 . The PMTs can be considered as the most important source of this background.
The final evaluation of the background due to (γ, n) interactions as well as to spontaneous fission in PMTs will be given in Sec. 9.6.

Radon background
In 2010-11, the Borexino scintillator was subjected to several purifications with water extraction (WE) procedures followed by nitrogen stripping. These operations were mainly aimed to remove the 85 Kr contamination and to lower as much as possible the 210 Bi one. During the purification campaigns, some radon entered the detector due to the emanation in the purification system. The 222 Rn has a lifetime of τ = 5.52 days and, after the operations, the correlated backgrounds typically disappear in a time window of a couple of weeks, leaving the corresponding amount of 210 Pb in the detector. An easy approach to precisely evaluate the Rn contamination level during the operations is to observe the 214 Bi(β − ,Q = 3.272 MeV)-214 Po(α, Q = 7.686 MeV) coincidence, among the radon daughters. By following this approach, we have estimated a contamination factor 100-1000 larger than the typical amount in the Borexino scintillator, during the WE period. For this reason, these transition periods are not used in solar neutrino studies, but, with some care can be used forν e analysis.
It is precisely the 214 Bi-214 Po coincidence that can become a potential source of background for the IBD selection. This coincidence has a time constant of τ = (236.0 ± 0.5) µs ( 214 Po lifetime [123]), very close to the neutron capture time in PC. The α-particles emitted by 214 Po usually have a visible energy well below the neutron capture energy window. However, in 1.04 × 10 −4 and in 6 × 10 −7 of cases, the 214 Po decays to excited states of 210 Pb and the α is accompanied by the emission of prompt gammas of 799.7 keV and of 1097.7 keV, respectively (see Table XI). In liquid scintillators, the γ of the same energy produces more light with respect to an α-particle. Therefore, for these (α + γ) decay branches, the observed light yield is higher with respect to pure αdecays and is very close to the neutron capture energy window. As already mentioned in Sec. 3.4, the Borexino liquid scintillator offers a possibility to discriminate highly ionizing particles (α, proton) from particles with lower specific ionization (β − , β + , γ) by means of pulse shape analysis. In the case of these rare branches, the scintillation pulses from the α and γ decays are so close in time that they practically overlap and result to be partially α-like. In order to suppress this background to negligible levels, during the purification periods we have increased the lower limit on the charge of delayed signal and applied a slight pulse shape cut as described below in Sec. 7. The final evaluation of the radon correlated background passing the optimized IBD selection cuts will be given in Sec. 9.7.
characterized with τ = (425.1 ± 1.5) ns of the 212 Po decay [123]. The 212 Bi is a β − emitter with Q = 2.252 MeV, while the α of 212 Po decay has 8.955 MeV energy. Given the short time coincidence, they could be a potential source of background for the IBD candidates, searched among the double cluster events, when both the prompt and delayed fall within one 16 µs DAQ gate (Sec. 3.1). This kind of events was included in the geoneutrino analysis for the first time (Sec. 7.2). Fortunately, the 212 Po is not giving any (α + γ) decay branch, so its effective energy distribution is below the neutron capture peak. Moreover, being a pure α decay, it can be easily recognised and rejected with a proper pulse shape analysis. The final evaluation of this background passing the optimized selection cuts will be given in Sec. 9.8.

DATA SELECTION CUTS
In this section we describe the cuts for the selection of antineutrino candidates and the process of their optimization. The vetoes applied after different muon categories are described in Sec. 7.1. Sections 7.2 and 7.3 deal with the definitions of the time and spatial correlation windows between the prompt and delayed IBD candidates. The application of the α/β discrimination techniques (Sec. 3.4) on the delayed, in order to suppress the radon background (Sec. 6.5), is shown in Sec. 7.4. Optimization of the energy cuts for the prompt and delayed signals is shown in Sec. 7.5, while the selection of the Dynamic Fiducial Volume in Sec. 7.6. The so-called multiplicity cut to suppress the background due to undetected muons with multiple neutrons and some noise is explained in Sec. 7.7. Finally, Sec 7.8 summarizes all the optimized values for all cuts.

Muon vetoes
Muons and related spallation products represent an important background for geoneutrino measurement, as described in Sec. 6.1. In the data selection, we first remove all categories of detected muons, as defined in Sec. 3.2. After different types of detected muons, different kinds of vetoes are applied. They are described in this section and are schematically shown in Fig. 23.

Veto after external muons
Among different kinds of spallation products of external muons, only cosmogenic neutrons can  penetrate inside the scintillator and thus present a background for geoneutrino analysis. The neutron capture time in Borexino is (254.5 ± 1.8) µs [77]. Therefore, a 2 ms veto, i.e. about 8 times the neutron capture time, is applied after all external muons detected by the OD.

Veto after internal muons
For internal muons, in addition to fast neutrons, 9 Li, 8 He, and 12 B isotopes (Table IX) are potential background sources for IBD selection as well. In the past analyses, a 2 s veto of the whole detector has been applied after all categories of internal strict and special muons. Since 2 s is nearly 8 times the lifetime of the longest-lived isotope ( 9 Li), this background was effectively eliminated for the price of about 10-11% loss of exposure. In this analysis, we reduce the exposure loss to 2.2% by introducing different categories of vetoes.

Veto after internal muons • 2 s veto of the whole detector
We apply a conservative 2 s veto of the whole detector after internal strict muons, special muons, and FADC muons that did not trigger the OD (not tagged by MTF flag). These represent 6.3% of all internal large muon category, which includes some noise events as well (Sec. 3.2). In addition, we apply this kind of veto after the so-called (µ + n) muons, i.e. those internal muons that triggered the OD (TT1 & BTB4, Sec. 3.1) and were followed by at least one neutron observed in the following event of TT128. These muons represent only 1.8% of all internal muons and have a higher probability to produce 9 Li events with a detectable decay neutron. The (µ + n) muon sample was used to characterize the veto parameters after the independent sample of (µ -n) muons, i.e. MTF/BTB4 internal muons that are not followed by any neutron, as described below.
• 2 ms veto of the whole detector It was observed that (µ + n) muons that also produce IBD-like hadronic background, always have more than 8000 decoded hits (Fig. 24). In the following text, we apply the notation (µ + n) ≥8000 for this and similar muon types. Muons producing less than 8000 decoded hits are typically not crossing the scintillator and are passing only through the buffer, where neutrons cannot be effectively detected due to low light yield. Thus, (µ -n) <8000 have little chance to produce 9 Li events with a detectable decay neutron. We apply a conservative 2 ms dead time after them, suppressing potential fast neutrons with negligible exposure loss. These muons represent 57.8% of all internal muons. Muons passing through the scintillator, i.e. with ≥8000 decoded hits, have high probability that neutron from a potential 9 Li decay would be detected. Thus, for these muons, a 2 ms veto is not sufficient. For the (µ -n) ≥8000 muon category that represents 34.1% of all internal muons, we apply a veto reduced from 2 s to 1.6 s, after which only 0.2% of 9 Li candidates survive. In addition, only 10% of the observed 9 Li background is produced due to (µ -n) ≥8000 muons, as it will be explained in detail in Sec. 9.3.
Muons in Borexino can be tracked based on their reconstructed entry and exit points in the OD and ID [77]. We consider the muon track to be reliable only when all the four track points are reconstructed. When this is not the case, the 1.6 s veto is applied to the whole detector, for 17.1% of all internal muons.

• 1.6 s cylindrical veto
The application of a cylindrical veto around the muon track, as schematically shown in Fig. 25(a), instead of vetoing the whole detector, can further increase the exposure for geoneutrino analysis. This kind of veto is applied to (µ -n) ≥8000 muons for which all the four muon track points are reconstructed.
The radius of the cylindrical veto is set by studying the lateral distance between the muon and IBD-like prompt and delayed observed within 2 s after the passage of a (µ + n) muon with four track points ( Fig. 25(b)). Within the observed statistics, this distance is very similar for the prompt and the delayed. 97.7% of the prompts (on which the DFV cut is applied) lie within a 3 m radius from the muon track in the IBD selection (Sec. 7.6). Since the lateral distribution of the muon daughters is expected to be the same for (µ -n) muons as well, a cylindrical veto of 1.6 s duration with 3 m radius is applied for (µ -n) ≥8000 muons, which constitute 27% of all internal muons.
The resulting exposure loss of 2.2% after all muon vetos was calculated using a Monte Carlo simulation. In total, 74 million point-like events were generated homogeneously in the dynamical fiducial volume for this study (Sec. 7.6), following the changing shape of the IV (Sec. 3.3). After considering the GPS times of all the muons and the track geometry reconstructed for (µ -n) ≥8000 muons, the relative exposure loss was calculated as the fraction of the events removed by all vetoes.

Time coincidence
The coincidence time window between the prompt and the delayed (dt) is an important background-suppressing cut. It is implemented based on the neutron capture time that was measured during the calibration campaign with the 241 Am- 9 Be neutron source to be (254.5 ± 1.8) µs [77]. Considering the 16 µs DAQ window, followed by an electronics dead time of 2-3 µs (Sec. 3.1), one has to consider separately the case when the prompt and delayed are either two separate events/triggers with single cluster each (similar to the event in top left of Fig. 4), or are represented by two clusters in a single event, as shown in Fig. 26.

Single cluster events
For this category of IBD candidates, the coincidence time window is between dt min = 20 µs and dt max = 1280 µs. The lower threshold guarantees that the delayed signal can trigger after the dead-time of the prompt trigger. The dt max corresponds to about five times the measured neutron capture time. This time window covers 91.8% of all IBD interactions.

Double cluster events
This category of IBD candidates is included in the present analysis for the first time. In this case, the coincidence time window is between dt min = 2.5 µs and dt max = 12.5 µs. The inclusion of double cluster events led to a 3.8% increase in the IBD tagging efficiency.
The lower threshold was optimized by studying the cluster duration of prompts from the MC of reactor antineutrinos (Sec. 8), which spectrum extends up to about 10 MeV. It guarantees that even for the highest energy prompt, after the dt min , there is no light that could alter the hit pattern of the delayed.
The dt max was set considering the variable position of the trigger-generating cluster (prompt) inside the DAQ window, that occurred during the analyzed period due to the changes in the trigger system. At the same time, it guarantees that the delayed can always have cluster duration of up to 2.5 µs before the end of the DAQ window.

Space correlation
Similar to dt, the spatial distance between the prompt and the delayed (dR) is also an important background-suppressing cut. In the previous analyses, dR = 1 m was used. The reconstructed distance between the prompt and the delayed is larger with respect to the distance between their respective points of production. This is caused predominantly by these effects: • Interaction of gammas: the gammas, from the positron annihilation and the neutron capture, interact in the LS mostly through several Compton scatterings. Thus, their interaction is intrinsically not point-like and the barycenter of the cloud of these Compton electrons is not identical to the point of the generation of the gammas.
• Position reconstruction: the position reconstruction (10 cm at 1 MeV) further smears the reconstructed positions.
In the optimization of the dR cut, two main aspects have to be considered.
On one hand, in the prompt-delayed reconstructed distance, shown in Fig. 27 for the MC sample of geoneutrinos (Sec. 8), the efficiency is quickly dropping below 1 m. On the other hand, with the increasing dR, the accidental background is also strongly increasing (Fig. 28). In order to find the optimized value, we have generated thousands of MC pseudo-experiments as described in Sec. 10. The cuts were set to optimized values, while the dR cut was varied. The dR cut was then set to 1.3 m, within the interval, where the variation in the expected precision of the geoneutrino measurement is small, as shown in Fig. 28.

Pulse shape discrimination
In the previous geoneutrino analyses, a Gatti cut G < 0.015 was applied on the delayed for efficient rejection of α-like events from radon-correlated background (Sec. 6.5). The cut value was set using gammas from neutron-captures from the 241 Am-9 Be calibration source and cross checked with the MC. In this analysis, an MLP cut > 0.8 was applied to the delayed using the better α/β discrimination power of the MLP (Sec. 3.4) when compared to the Gatti. The cut threshold was chosen based on the 241 Am-9 Be calibration data as shown in Fig. 29. It was found that only 3.7·10 −3 of the neutron-capture γs remains below the threshold of 0.8.

Energy cuts
Energy cuts were applied to the the prompt and delayed to aid the identification of IBD signals. The analysis was performed using the charge energy estimator N pe (Sec. 3.1. The intervals in the respective charges Q p and Q d are optimized as explained in this Section.

Charge of prompt signal
The energy spectrum of the prompt E p starts at ∼1 MeV, which corresponds only to the two 511 keV annihilation gammas (Eq. 10). Therefore, the threshold on the prompt charge was set to Q min p = 408 p.e. This threshold value corresponds to approximately 0.8 MeV and remains unchanged from previous analyses. No upper limit is set on the charge of the prompt candidate.

Charge of delayed
The delayed signal can be either due to a 2.22 MeV (n-capture on 1 H), or a 4.95 MeV gamma (n-capture on 12 C) with about 1.1% probability, as described in Sec. 4. The corresponding values in the N pe variable are 1090 p.e. and 2400 p.e., respectively, as measured in the detector center and shown in Fig. 13. However, at large radii, gammas can partially deposit their energy in the buffer, which decreases the visible energy. Consequently, the γ-peak develops a low-energy tail and even the peak position can shift to lower values (Fig. 30).
In the last geoneutrino analysis [18], Q min d was set to 860 p.e., a conservatively large value because of the radon-correlated background, particularly due to 214 Po(α + γ) decays. This was discussed in Sec. 6.5. In this analysis, Q min d was decreased to 700 p.e. based on the improved performance of α/β separation with MLP (Sec. 7.4), which improved the ability to suppress 214 Po(α + γ) decays. This cut was applied to all data with the exception of the water-extraction period, which had an increased radon contamination. In this case the Q min d 860 p.e. was retained. The Q min d cut was not decreased below 700 p.e. for the following reasons: 1. The N pe spectrum of accidental background increases at lower energies, as shown in Fig. 37. 2. The end point of the α peak from the main 214 Po decay in the radon-correlated background is ∼600 p.e., as described in Sec. 9.7.
In this analysis, we include the neutron captures on 12 C (4.95 MeV) as well. Consequently, Q max d was increased from 1300 p.e. (2.6 MeV) to 3000 p.e. (≈5.5 MeV).

Dynamical fiducial volume cut
The shape of the Borexino IV, that is changing due to the presence of a small leak, can be periodically reconstructed by using the data (Sec. 3.3). A Dynamical Fiducial Volume (DFV) cut, i.e. a cut on the distance of the prompt from the IV, is applied along the reconstructed IV shape. In the previous analysis [18] a conservative cut of 30 cm has been applied to account for the uncertainty of the IV shape reconstruction and the potential background coincidences near the IV. In this analysis, the DFV cut was decreased to 10 cm, which leads to a 15.8% relative increase in exposure. This choice is justified below.
The geoneutrino sensitivity studies (Sec. 10) were performed for different combinations of the DFV cut and the Q d min , as shown in Fig. 31  that no excess of the IBD candidates is observed at large radii (close to the IV), as it will be shown in Sec. 11.1. For a given DFV cut, the choice of Q min d , shown on the x-axis, does not strongly influence the geoneutrino expected precision. Also note, that for 10 cm (final choice) and no DFV cut, the sensitivity to geoneutrinos is nearly the same.

Multiplicity cut
The multiplicity cut requires that no additional "high-energy" (N pe > 400 p.e.) event is observed within ±2 ms around either the prompt or the delayed candidate.
This cut is designed to suppress the background from undetected muons, for example, neutron-neutron or buffer muon-neutron pairs. This justifies the selected time window which is nearly 8 times the neutron capture time. The charge cut was lowered to account for those neutrons that are depositing their energy partially in the buffer. Thanks to a high radiopurity of the LS, the corresponding exposure loss due to accidental coincidences of IBD candidates with N pe > 400 p.e. events within 2 ms is of the order of 0.01%, which is negligible.

Summary of the selection cuts
The summary of all the optimized selection cuts is listed in Table XII. The spectral fit of the prompt charge Q p (Sec. 10.1), relies extensively on the use of MC-constructed probability distribution functions (PDFs) which represent the shapes of geoneutrino signal and, with the exception of accidental coincidences (Sec. 9.4), all backgrounds.
The construction of these PFDs is described in Sec. 8.1. The Geant-4 based MC of the Borexino detector was tuned on independent data acquired during an extensive calibration campaign with radioactive sources [76] and is described in detail in [27]. For the antineutrino analysis, the calibration with 241 Am-9 Be neutron source is of particular importance, since the delayed IBD (Sec. 4) signal is represented by a neutron. The comparison of the neutron spectra from the 241 Am- 9 Be calibration source at the detector center and at (0, 0, -4) m position inside the detector is shown in Fig. 30.

Monte Carlo spectral shapes
Once the full G4Bx2 based MC code is working reliably and the origin of the signal and backgounds is known, it is, in principle, easy to simulate the PDFs that incorporate the detector response and that can be used directly in the spectral fit (Sec. 10.1).
The simulated signal and backgrounds follow the same experimental conditions as observed in real data, including the number of working channels, the shape of the IV, and the dark noise, as described in [27]. Each run of the complete data set from December 2007 to April 2019 is simulated individually. After the simulation, the optimized geoneutrino selection cuts (Sec. 7.8) are applied as in the real data.
For antineutrinos, pairs of positrons and neutrons were simulated. The neutron energy spectrum is taken from [124], while for the positrons the energy spectra as discussed in Sec. 5 are used. The antineutrino energy spectra are transformed to positron energy spectra following Eq. 10. In particular, for geoneutrinos, the energy spectra as in Fig. 18 are used. Individual spectra from 232 Th and 238 U chains were also simulated, so that they can be weighted according to the expected R s ratio (Eq. 18) for different geological contributions. For reactor antineutrinos, calculated energy spectra "with and without 5 MeV excess" as in Fig. 19 are used as MC input. The resulting PDFs are shown in Fig. 32.
The MC-based PDFs for non-antineutrino backgrounds are shown in Fig. 33. The inputs are discussed for (β + n)-decays of cosmogenic 9 Li in Sec. 6.1, for (α, n) background in Sec. 6.3, while for atmospheric neutrinos in Sec. 5.4.
The PDFs of prompts due to antineutrinos from a hypothetical georeactor (Sec. 5.5) are shown in Fig. 34. We compare the shapes (in all cases normalized to one) for the non-oscillated spectrum and for the oscillated cases with the georeactor placed at three different positions GR1, GR2, and GR3 as shown in Fig. 21(a). As it can be seen, the shapes of the prompt energy spectra are almost identical.

Detection efficiency
The detection efficiencies for geoneutrinos (ε geo , ε Th , ε U ), reactor antineutrinos (ε rea ), and antineutrinos from a hypothetical georeactor (ε georea ) are summarized in Table XIII. They represent a fraction of MC events passing all the optimized data selection cuts (Sec. 7.8) from those generated in the FV of this analysis (10 cm DFV cut). The errors due to the FV definition and the position reconstruction resolution are included in the calculation of the systematic uncertainty, as it will be discussed in Sec. 11.3.

EVALUATION OF THE EXPECTED SIGNAL AND BACKGROUNDS WITH OPTIMIZED CUTS
In this Section the evaluation of the number of expected antineutrino signal and background events after the optimized selection cuts (Table XII) is described. The data set and the total exposure are presented in Sec. 9.1.
The number of expected antineutrino events from different sources, based on the estimated antineutrino signals as in Table VIII, is presented in Sec. 9.2. The next sections treat the non-antineutrino background, following the structure of Sec. 6, where the physics of each of these background categories is described. In particular, cosmogenic background is discussed in Sec. 9.3, accidental background in Sec. 9.4, (alpha, n) interactions in Sec. 9.5, (γ, n) and fission in PMTs in Sec. 9.6, radon correlated background in Sec. 9.7, and finally 212 Bi-212 Po background in Sec. 9.8. Section 9.9 summarizes the expected total number of non-antineutrino background events.

Data set and exposure
In this analysis the data taken between December 9, 2007 and April 28, 2019, corresponding to t DAQ = 3262.74 days of data taking, are considered. The average life-time weighted IV volume and the FV used in this analysis (10 cm DFV cut) are V IV = (301.3 ± 10.9) m 3 and V FV = (280.1 ± 10.1) m 3 , respectively, after taking the changing shape of the IV into account (Sec. 3.3). These correspond to the average IV and FV mass of m IV = (264.5 ± 9.6) ton and m FV = (245.8 ± 8.7) ton, respectively, considering the scintillator density of ρ LS = (0.878 ± 0.004) g cm −3 . The total exposure after the cosmogenic veto (Sec. 7.1) and in the FV is E = (2145.8 ± 82.1) ton × yr, considering the systematic uncertainty on position reconstruction (Sec. 11.3). This can be expressed as E p = (1.29 ± 0.05) ×10 32 protons × yr, using the proton density in Borexino LS of N p = (6.007 ± 0.001) × 10 28 protons ton −1 . Applying the geoneutrino detection efficiency described in Sec. 8.2, the effective exposure for the geoneutrino detection reduces to E = (1866.4 ± 78.4) ton × yr and E p = (1.12 ± 0.05) ×10 32 protons × yr.

Antineutrino events
This Section summarizes the number of expected antineutrino events detected with the optimized selection cuts, assuming the expected antineutrino signals as given in Table VIII. An overview is given in Table XIV.

Cosmogenic background
The various cosmogenic backgrounds in Borexino which affect the geoneutrino analysis were explained in Sec. 6.1, while the analysis and technical evaluation of these backgrounds is explained in this Section.

Hadronic background
The hadronic background, expected to be dominated by 9 Li (Sec. 7.1), which remains after the detected muons out of the vetoed space and time is evaluated here. We first study the time and spatial distributions of the detected 9 Li candidates with respect to the parent muon. After that, this background is evaluated for the different kinds of internal muons according to the respective vetoes applied after them, as previously described in Sec. 7.1 and summarized in Fig. 23.

Time distribution dt Li−µ
First, a search for IBD-like signals, passing the optimized selection cuts (Table XII), is performed after the category of muons for which 1.6 or 2.0 s veto is applied. We perform this search starting from 2 ms after each muon, in order to remove cosmogenic neutrons. In total, we found 305 such IBD-like candidates, dominated by 282 candidates afer (µ + n) muons. The Q Li p charge energy spectrum of the prompts is compatible with the expected MC spectrum of 9 Li, as it is shown in Fig. 35(a). The distribution of the time differences between the prompt and the preceding muon, dt Li−µ , is shown in Fig. 35(b). As it can be seen, no events are observed after dt Li−µ > 1.6 s. The decay time τ is extracted by performing an exponential fit to the dt Li−µ distribution and is found to be (0.260 ± 0.021) s. This is compatible with τ9 Li = 0.257 s decay time of 9 Li.

Spatial distribution dR Li−µ
The distance of the 9 Li prompt from the muon track, dR Li−µ , is studied for (µ + n) muons with reliably reconstructed tracks. It is shown for 85 candidates in Fig. 35(c). This distribution is fit with the convolution of an exponential (with a characteristic length λ) with a Gaussian with parameters µ and σ, and a normalisation factor n: The fit results in σ = (0.35 ± 0.11) m, which represents well the combined position reconstruction of the muon track and the prompt, and in λ = (0.68 ± 0.24) m.
Considering these values and the fit function in Eq. 30, 2.75% of 9 Li prompts would be reconstructed out of the cylinder with 3.0 m radius around the muon track. Below, the hadronic background after different muon-veto categories (Fig. 23) is evaluated, considering the dt Li−µ and dR Li−µ parametrisations described above. For 27.0% of the muons with a reliably reconstructed track, only a cylinder with 3 m radius is vetoed for 1.6 s. In the [2 ms, 1.6 s] time window after these muons, 7 IBD-like events were observed in the whole detector. Thus, after 1.6 s, we expect 0.015 +0.013 −0.011 events distributed in the whole detector. Of these, 2.75% would be reconstructed out of the cylindrical veto, as mentioned before. This means, our expected background of this category can be conservatively set to (0.20 ± 0.08) events. By restricting the veto from the whole detector to only the cylindrical volume, the exposure increases by 1.47%, corresponding to (2.2 ± 0.2) expected IBD events. We observe one additional candidate.
• 2 ms veto of the whole detector: For 57.8% of all muons which have a lower probability to produce detectable hadronic background ((µ -n) <8000 ), we restrict the time veto of the whole detector to 2 ms, to veto only the cosmogenic neutrons. Consequently, the exposure increases by 6.6%, corresponding to (9.7 ± 0.8) expected IBD candidates. Seven candidates were observed, well within the expectation. However, this does not guarantee that we did not introduce additional 9 Li background, that is estimated as follows.
Muons with less than 8000 hits produce less light because they pass mostly through the buffer region, where the neutrons are typically not detected or are below the threshold of this analysis. It is reasonable to assume that the production ratios for 9 Li and cosmogenic neutrons are the same for µ <8000 and µ >8000 muon categories, since the muons have typically the same energies and the traversed media (LS and the buffer) have nearly the same density. It is also reasonable to assume, that for the µ <8000 , the detection efficiency of the corresponding cosmogenic neutrons and neutrons from 9 Li decays are the same. Thus, the equality of the ratios of the number of observed 9 Li candidates (with decay neutron, N Li ) and cosmogenic neutrons (N n ) should hold: (N Li = 305) >8000 (N n = 8.6 × 10 5 ) >8000 = (N Li ) <8000 (N n = 9181) <8000 .
(31) From this equality, the expected number (N Li ) <8000 is 3.2 events. This is the number of expected 9 Li events produced by muons with less than 8000 hits, independent of whether this muon was followed by a neutron or not. This is a conservative number for the background estimation due to the (µ − n) <8000 muons, after which we apply the reduced 2 ms cut. Even if they represent about 99% of all µ <8000 muons (Fig. 24), (µ + n) <8000 muons can be expected to have a higher probability to produce observable 9 Li than (µ − n) <8000 muons. However, we do not observe any 9 Li candidate for (µ + n) <8000 muons. Thus to summarize, our expected 9 Li background for (µ − n) <8000 muons is 3.2 ± 1.0 events, which includes larger systematic error.
In total, after summing all the contributions, the ex-pected 9 Li background within our golden IBD candidates is (3.6 ± 1.0) events.

Untagged muons
The (0.0013 ± 0.0005)% mutual inefficiency of the strict muon flags shown in Sec. 3.2 corresponds to (195 ± 75) undetected muons in the entire data set. Following the discussion in Sec. 6.1, these muons could eventually cause background of three types: • µ + µ: Considering the small amount of undetected muons in the entire data set, the probability that two undetected muons would fall in a 1260 µs time window of the delayed coincidence is completely negligible.
• µ + n: The most dangerous are pairs of buffer muons (possibly fulfilling the Q p cut) followed by a single neutron (multiple neutrons are removed by the multiplicty cut). The probability that a (µ + n) pair falls within the IBD selection cuts, evaluated on the subset of MTB muons followed by the TT128 trigger dedicated to neutron detection, is found to be (9.7 ± 0.003) × 10 −5 . Hence, there will be (0.019 ± 0.007) events of this kind in the IBD sample due to the untagged muons.
• Muon daughters: After 1.497 × 10 7 internal muons we have observed 305 IBD-like background events in a [2 ms, 1.6 s] time window, that covers 99.02% of IBD-like candidates of the same type. Therefore, after (195 ± 75) undetected muons, we can estimate to have (0.0040 ± 0.0015) IBD-like events created any time after these muons and falling within the selection cuts.
Summing all the three components, the estimated background originating from untagged muons is (0.023 ± 0.007) events in the IBD sample. We note that this is a very small number.

Fast neutrons
As described in Sec. 6.1, undetected muons that pass the WT or the surrounding rocks, can produce fast neutrons that can give IBD-like signals. Fast neutrons from cosmic muons were simulated according to the energy spectrum from [125]. We have found that the eventual signal from a scattered proton follows in nanosecond time scale after the neutron production, that is simultaneous with the muon signal. Considering the data structure detailed in Sec. 3.1, this time range dictates the data selection cuts, as described below, in order to search for fast-neutron related IBD-like signals after the detected external muons pass the WT. Knowing the fraction of these muons creating IBD-like background, we can estimate the fast neutron background from the undetected muons that pass the WT. MC simulations are used to obtain an estimation of fast neutron background due to the muons passing through the surrounding rock and not the detector. Both estimations are given below.
Water Tank muons: In this search, the signal in the ID should correspond to a scattered proton, which is not tagged by the muon Inner Detector Flag (IDF). Without the proton signal in the ID, the external muon would be a TT2 & BTB4 event. The presence of the ID signal can, with lower than 100% efficiency, turn the muon to be a TT1 & BTB4 event. Therefore, we search for two kinds of coincidences: • The prompt signal is an internal muon TT1 & BTB4 that is not tagged by the IDF. The delayed signal is a neutron cluster found in the TT128 gate which is opened immediately after the muon.
• The prompt signal is an external muon TT2 & BTB4 that has a cluster inside the ID, and is not tagged by the IDF. The delayed signal follows within 2 ms as a point-like TT1 & BTB0 event.
This search was done with relaxed energy, dt, and dR cuts, without any DFV or multiplicity cuts. This yielded 25 coincidences of the first kind and 12 coincidences of the second kind. However, only one coincidence satisfied all the geoneutrino selection cuts. The amount of these coincidences in the IBD data sample can be due to muons that go undetected by the OD. The average inefficiency of MTF with respect to the MCF and IDF flag is 0.27%, as shown in Table III. This gives an upper limit of 0.013 IBD-like coincidences at 95% C.L. due to the fast neutrons from undetected muons crossing the WT.

Surrounding rocks:
In order to study the fast neutron background due to muons passing through the rocks surrounding the detector, we used the Borexino MC with the initial flux and energy spectrum of neutrons and their angular distributions taken from [125] for the specific case of LNGS. The total statistics of MC-generated neutrons corresponds to 3.3 times the exposure of this analysis. Fast neutrons with energies in diapason 1 MeV -3.5,GeV were simulated on the surface of the Borexino outer water tank. Full simulation with tracking of each scintillation photon was done for the fast neutrons and other particles penetrating inside the ID. Finally, we obtain only one IBD-like event for neutrons from the rock passing the optimized IBD selection cuts, which corresponds to an upper limit of the corresponding background in our geoneutrino analysis of <1.43 events at 95 C.L.

Accidental coincidences
In order to evaluate the amount of accidental coincidences in the antineutrino sample, coincidence events were searched for in the off-time interval dt = [2 s, 20 s] and were then scaled to the 1270 µs duration of the geoneutrino selection time window (dt = [2.5 µs, 12.5 µs] + [20 µs, 1280 µs], Sec. 7.8). In this scaling, a suppression factor due to the muon veto must be considered, as explained below.
To evaluate the accidental background rate which is not biased by the cosmogenics, it is required not only for the prompt, but also for the delayed not to be preceded by a muon within 2 s. This means, that once the prompt is accepted, there is no preceding muon within 2 s before the prompt. After the prompt, as dt between the prompt and a potential delayed increases, so does the probability that the delayed will be discarded due to the muon falling in between the prompt and the delayed. For time intervals longer than 2 s, this probability becomes constant, because the muon veto is of 2 s. This behavior is illustrated in Fig. 36(a), which shows the time distribution between the prompt and delayed accidental signals in a time window dt = [2 ms, 4 s]. One can see that until 2 s, there is a decrease, while after 2 s the distribution is flat. Note that this plot was constructed with relaxed selection cuts and serves only to demonstrate the suppression factor that depends only on the muon rate r µ and the muon veto time. The fit function for this distribution in the interval dt = [2 ms, 2 s] is as follows: where r acc 0 would be the rate of accidental background with relaxed cuts and without the muon veto supression factor. After 2 s, the exponential suppression factor becomes constant and consequently the fit function for dt > 2 s acquires a constant form: When fitting the spectrum with relaxed cuts, r µ = (0.0501 ± 0.0016) s −1 was obtained. This is compatible with the measured rate of internal muons of (0.05311 ± 0.00001) s −1 . The validity of this behaviour has been also verified by a MC study. The suppression factor exp(−r µ · dt) for the dt of the real IBD selection is larger than 0.99993 and thus can be neglected. However, for the times dt > 2 s, the suppression factor is 0.896 ± 0.003, conservatively considering also the difference between the r µ resulting from the fit in Fig. 36(a) and just by measuring the rate of internal muons.
In order to determine the rate of accidental coincidences r acc 0 for the geoneutrino measurement, the dt = [2 s, 20 s] distribution of 49004 events selected with optimized IBD selection cuts was constructed, as shown in Fig. 36(b). This distribution is, as expected, flat and is fit with a function: The exponential suppression factor is set to 0.896 ± 0.003, the value discussed above, since the muon veto conditions are the same as in the accidental search with relaxed cuts. The resulting r acc 0 is (3029.0 ± 12.7) s −1 . This means that the number of accidental coincidences among our IBD candidates can be estimated as r acc 0 × 1270 µs, that is (3.846 ± 0.017) events.
The N pe spectra of the prompt and delayed signals of the accidental coincidences, selected with optimized geoneutrino cuts in dt = [2 s, 20 s] time window, are shown in Fig. 37.

(α, n) background
The (α, n) background evaluation is done in three stages. First, the amount of α particles that could initiate this interaction is estimated. In Borexino, the only relevant isotope is 210 Po that is found out of equilibrium with the rest of 238 U chain [20]. In the energy region of 210 Po (N pe = 150 -300 p.e.), α-like particles (MLP < 0.3) reconstructed in the DFV of the geoneutrino analysis are selected. The evolution of the weekly rates of such events for the whole analyzed period is shown in Fig. 38.
The mean rate of R DFV ( 210 Po) = (12.75 ± 0.08) counts/(day·ton) is used to evaluate the (α, n) background from the 210 Po contamination of the LS.
In the second stage, the neutron yield, i.e. the probability that 210 Po α would trigger an (α, n) reaction in the LS, was calculated with the NeuCBOT program [126,127,128,129], which is based on the TALYS software for simulation of nuclear reactions [130,131,122]. Only PC was considered as a target material.
The contribution from PPO is negligible, as its relative mass fraction is small. According to a recent article [120], the analytical calculation of the (α, n) cross section with TALYS provides a result similar to the experimental data for energies of the 210 Po α particles. This fact permits to apply as a relative uncertainty of our calculation the 15% uncertainty of the experimental data [120]. Taking this into account, the neutron yield Y n of the (α, n) reaction in the PC is found to be (1.45 ± 0.22) × 10 −7 neutrons per a single 210 Po decay. Note that the corresponding value used in the previous Borexino geoneutrino analysis, based on [132], was approximately three times smaller.
The final calculation of the number of IBD-like coincidences N (α,n) triggered by the neutrons from the (α, n) reaction over the whole analysis period can be computed using the following formula: where E = (2145.8 ± 82.1) ton × yr (Sec. 9.1) is the exposure and ε (α,n) = 56% is the probability of the (α, n) interaction to produce an IBD-like signal passing all selection cuts, obtained with a full G4Bx2 MC study. Based on this evaluation, the expected (α, n) background due to the 210 Po contamination of the LS is (0.81 ± 0.13) events.
Another potential source of background are (α, n) reactions due to 210 Po decays in the buffer. Based on a G4Bx2 MC study, we have found that these interactions occurring in the outer buffer have negligible probability to create IBD-like background. However, for the interactions occurring in the inner buffer, this probability was estimated to be 0.23% and the energy spectrum of prompts is very similar to the (α, n) from the LS (middle panel of Fig. 33). It is extremely difficult to determine the 210 Po contamination of the buffer, since the α peak is completely quenched below the detection threshold. In 2009 we have estimated this contamination as <0.67 mBq/kg [16] by employing the samples of buffer liquids in the center of the Counting Test Facility of Borexino [133], that not any more operational. This limit is several orders of magnitude above the contamination of the LS. DMP quencher, that is only present in the buffer, is considered to be the main source of the 210 Po contamination in the buffer. In January 2010, the DMP concentration in the buffer was reduced to 2 g/l (the original contamination was 5 g/l), as discussed in Sec. 3. Since then, no further operations have been performed with the buffer and the 210 Po contamination is expected only to decay (τ = 199.6 day) and to be sup- pressed in April 2019 by a factor 3.9 × 10 −8 . In the present analysis, the estimated upper limit for this contamination is 0.14 mBq/kg, which corresponds to an upper limit of 2.6 background events (from which only 0.3 events in the period from January 2010). We note however, that the original estimate of the 210 Po rate in the buffer is very conservative, because of high risk of contamination of the samples during their handling. As it will be discussed in Sec. 11.1, the golden IBD candidates are evenly distributed in time and no excess close to the IV is observed.

(γ, n) interactions and fission in PMTs
In order to obtain an upper limit to the possible background from (γ, n) reactions in the Borexino scintillator or in surrounding materials, we counted all the registered events with energies higher that 3 MeV and we made the conservative assumption that they are only due to γ-ray interactions. Since the energy response of the detector is not uniform in space and time, an energy release of 3 MeV does not correspond to a unique value of the registered charge N pe . To consider this effect, 3 MeV γs have been generated with the G4Bx2 MC code following the detector status during the whole analyzed period. According to Fig. 39, a conservative charge threshold of 1200 p.e. was chosen and a correction of 5.4% for the inefficiency of the cut was then applied: 589,917 events were selected above 1200 p.e., resulting in 623,571 hypothetical γ-rays after the correction. Each of them can only interact with the deuterons that meets along its path before being absorbed: an estimation for this background is then obtained by multiplying the numbers of gammas for the deuteron density, the interaction cross section, and the gamma's absorption length. According to the γ-ray attenuation coefficients calculated for the Borexino scintillator, the absorption length λ for a 3 MeV gamma is 29 cm and the capture cross section on 2 H is σ D = 1.6 mb. Since the deuteron density is ρ D = 7.8 × 10 18 atoms/cm −3 , the upper limit on the number of events N (γ,n) due to this background, taking into account the estimated detection efficiency ε (γ,n) = 50%, is: An attenuation length of 3λ was chosen to obtain the 95% C.L. A possible contribution of neutron capture on 13 C and 12 C nuclei was also considered, but it was found to be more that a factor 10 smaller and therefore, neglected. In PMTs we have determined a 238 U contamination of (31 ± 2) ppb in the glass and (60 ± 4) ppb in the dynodes. PMTs are located at about 6.85 m from the center of the detector. To estimate the background induced by spontaneous fission, we consider that for each PMT the glass accounts for about 0.3 kg and the dynodes for 0.05 kg. In addition, we take into account the subtended solid angle by the IV and the neutron attenuation while propagating from the PMTs to the IV, which is of the order of 1 m. The estimated number of neutrons reaching the scintillator is (0.057 ± 0.004) for the current exposure. The corresponding neutron-induced background will be negligibly small and we set for it a conservative upper limit of 0.057.

Radon background
In Section 6.5 we have discussed how the radon contamination of the LS can induce IBD-like background. Figure 40 demonstrates the increased Radon contamination during the water extraction period. A proper choice of the IBD selection cuts is extremely useful to reduce this kind of background and to safely include the WE period in the geoneutrino analysis.
The energy scale of the (α + γ) decays of 214 Po (Table XI) was evaluated. The energy scale in G4Bx2, including the overall light yield and the Birk's constant kB important in the description of quenching (see Sec. VII of [20]), is tuned based on the calibration with γ sources. Since the kB is particle dependent, the energy scale for αs must be further adjusted. For this purpose, the dominant pure α decay of 214 Po was simulated and compared to the Radon events from the data (selected via the 214 Bi 214 Po delayed coincidence tag), as demonstrated in Fig. 41. With both γ and α energy scales fixed, the spectra for (α + γ) 214 Po decays were simulated, as it is shown in Fig. 41. Since the overall radon statistics amounts to 1.1 × 10 5 decays, the events due to the 10 −7 branch can be neglected, while ∼11 events are expected from the 10 −4 branch, when 214 Po decays to 210 Pb in the first excited state. In this case, the de-excitation gamma is emitted along with an α-particle. In order to suppress these events, we keep the Q min d threshold fixed to 860 p.e. during the WE period, which effectively reduces this background by a factor of 10 3 . Application of the pulse shape cut  (Table XI) obtained using MC with the α energy scale tuned on pure α-decays (Fig. 41). The α decay (red line) and the two (α + γ) branches are shown: 10 −4 (green line) and 10 −7 (blue line) branch. The three spectra are normalised to have the same area.
(MLP > 0.8) on the delayed (Table XII) further reduces the background by a factor of 5-6. During the analysis of non-WE period, the 10 −4 branch also becomes negligible, hence we can lower Q min d to 700 p.e., safely above the 214 Po α-peak (Fig. 41). The total number of background events correlated with radon contamination is expected to be (0.003 ± 0.0010), which is completely negligible.

212 Bi-212 Po background
A MC study proves that the cut on Q min d = 860 p.e., adopted to reject the radon contamination during the WE periods, is effective in removing the 212 Bi-212 Po fast coincidences, as shown in Fig. 42. It can be seen that the end point of the 212 Po α peak is around 700 p.e. Therefore, a Q min d of 700 p.e., combined with the MLP pulse  shape cut on the delayed, makes the overall 212 Bi-212 Po background fully negligible in geoneutrino analysis.
9.9 Summary of the estimated non-antineutrino background events Table XV summarizes the expected number of events from all non-antineutrino backgrounds passing the optimised selection cuts listed in Table XII.

SENSITIVITY TO GEONEUTRINOS
This Section describes the Borexino sensitivity to geoneutrinos and the MC based procedure with which it was evaluated. In Sec. 10.1 the description of the basic ingredients of the analysis focuses on the spectral fit of the Q p spectrum. Section 10.2 describes the sensitivity tool that performs such fits on 10,000 Q MC p spectra, each 51 corresponding to a MC generated pseudo-experiment. The expected precision for the Borexino geoneutrino measurement as well as its sensitivity to the mantle signal is discussed in Sec. 10.3. This approach was also used in the optimization of the selection cuts, as mentioned in Sec. 7. The systematic uncertainties given in Sec. 11.3 are not included in the sensitivity studies and are only considered in the final results of Sec. 11. As it will be shown, the error on the geoneutrino measurement is largely dominated by the statistical error.

Geoneutrino analysis in a nutshell
The geoneutrino signal is extracted from the spectral fit of the charges of the prompts of all selected IBD candidates. Since the number N IBD of selected candidates is relatively small (in this analysis, N IBD = 154 candidates, see Sec. 11.1), an unbinned likelihood fit is used: where Q p is the vector of individual prompt charges Q i p , and index i runs from 1 to N IBD . The symbol θ indicates the set of the variables with respect to which the function is maximized, namely the number of events corresponding to individual spectral components with known shapes. In particular, we fit the number of geoneutrino and reactor antineutrino events as well as the number of events from several background components. The shapes of all spectral components are taken from the MC-constructed PDFs (see Figs. 32,33), with the exception of the accidental background, which can be measured with sufficient precision as shown in Fig. 37 (prompt spectrum). Some of the spectral components are kept free (typically geoneutrinos and reactor antineutrinos), while others (typically other than reactor antineutrino backgrounds) are constrained using additional multiplicative Gaussian pull-terms in the likelihood function of Eq. 37.
Naturally, the number of geoneutrinos is always kept free. One way of doing it is by having one free fit parameter for geoneutrinos, when we use 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 R s signal ratio of 0.27 (Sec. 5.2). Alternatively, 232 Th and 238 U contributions can be fit as two independent contributions. Additional combinations are of course possible. For example, in the extraction of the geoneutrino signal from the mantle (Sec. 11.5), we constrain the expected lithospheric contribution, while keeping the mantle contribution free.
The number of reactor antineutrino events is typically kept free. It is an important cross-check of our ability to measure electron antineutrinos, when we compare the unconstrained fit results (Sec. 11.2.1) with the relatively-well known prediction of reactor antineutrino signal (Sec. 5.3). In addition, an eventual constraint on reactor antineutrino contribution does not significantly improve the precision of geoneutrinos, as verified and discussed below. The constrained reactor antineutrino signal is however used when extracting the limit on the hypothetical georeactor (Sec. 5.5), as it will be discussed in Sec. 11.7.
Typically, we include the following non-antineutrino backgrounds in the fit: cosmogenic 9 Li, accidental coincidences, and (α, n) interactions. Atmospheric neutrinos are included in the calculation of systematic uncertainties, as it will be described in Sec. 11.3.
These background components are constrained in the fit, since independent analyses can yield the well constrained estimates of their rates, as they are summarized in Table XV for non-antineutrino backgrounds and in Table XIV given for atmospheric neutrinos.

Sensitivity study
A Monte Carlo approach was used in order to estimate the Borexino sensitivity to geoneutrinos, as well as to optimize the IBD selection cuts (Sec. 7). This socalled sensitivity study can be divided in the following four steps: • The arrays of charges of prompts for signal and backgrounds are generated from the PDFs including the detector response that were either created by the full G4Bx2 MC code (Sec. 8.1, Figs. 32 and 33) or measured, as for accidental background (Fig. 37, prompt spectrum). For each component, the number of generated charges is given by the expectations, as shown for antineutrino signals in Table XIV and for non-antineutrino backgrounds in Table XV.
• The generated spectra are fit in the same way as the data (Sec. 10.1), using in the fit the same PDFs that were used for the generation of these pseudoexperiments. This means, uncertainty due to the shape of the spectral components is not considered. This is justified by the fact, that Borexino's sensitivity to geoneutrinos is by far dominated by the statistical uncertainty.   (Table XIV) and represent the expected statistical uncertainty of the measurement.
• The procedure is repeated 10,000 times for each configuration. In each pseudo-experiment, the number of generated events for signal and all backgrounds is varied according to the statistical uncertainty.
• The distributions of ratios of the resulting fit value over the MC-truth (expected) value in each individual fit are constructed for the parameters of interest. For example, such a distribution for the ratio of the number of geoneutrinos resulting from the fit over the number of generated geoneutrino prompts should be centered at one (when there is no systematic bias), while the width of this distribution corresponds to the expected statistical uncertainty of the measurement.

Expected sensitivity
Using the sensitivity tool as explained in Sec. 10.2, the expected statistical uncertainty of the Borexino geoneutrino measurement in the presented analysis varies from (13.76 ± 0.10)% to (23.09 ± 0.17)%, depending on the expected signal for different geological models (Table XIV), as demonstrated in Fig. 43. This study assumes the Th/U chondritic ratio to hold. In the previous 2015 Borexino geoneutrino analysis [18], the statistical error was ∼26.2%.
The sensitivity of Borexino to measure the 232 Th/ 238 U ratio was also studied. As it is shown in Fig. 44, Borexino does not have any sensitivity to determine this ratio. Despite the input ratio assuming the chondritic value (considering the statistical fluctuations), the 232 Th/ 238 U ratio resulting from the fit has nearly a flat distribution for the 10,000 pseudo-experiments. This will be also confirmed by large 232 Th versus 238 U contours shown in Fig. 48(d) for the fit of the data with free 238 U and 232 Th components.
The sensitivity of Borexino to measure the mantle signal was studied using the log-likelihood ratio method [134] for the expectations according to four different geological models (CC, GC, GD, and FR, Table XIV). For each geological model, we have generated a set of 10,000 pseudo-experiments with the mantle geoneutrino component included. In addition, we have generated 1.2 million pseudo-experiments without the mantle contribution. In each data set, we have included the relatively-well known lithospheric contribution (Table VI), as well as the reactor antineutrino "without 5 MeV excess" (Table XIV)  non-antineutrino backgrounds (Table XV). Each pseudo-experiment from all five data sets (one without the mantle and four with mantle signal according to four geological models), are fit twice: with and without the mantle contribition. The best fit with the mantle contribution fixed to zero corresponds to the likelihood L{0}.
The fit with the mantle component left free results in the likelihood L{µ}. Obviously, for the data set without the mantle being generated, the two likelihoods tend to be the same. For the data sets with the mantle included, the L{0} tends to be worse than L{µ}: the bigger this difference, the better the sensitivity to observe the mantle signal.
We define the test statistics q (q ≥ 0): that we call q 0 for the data set without the mantle generated. The q 0 and the four q distributions for different geological models are shown in Fig. 45. The q 0 corresponds to the theoretical f (q|0) distribution: The four q distributions we fit with the f (q|µ) where Φ stands for a cumulative Gaussian distribution with mean µ and standard deviation σ. For high statistical significance, µ/σ is very large and Φ( µ σ ) → 1. In Figure 45 we show also q med = (µ/σ) 2 , the median value of f (q|µ), for the four different geological models. We express the Borexino sensitivity to measure the mantle geoneutrino signal, according to these four geological models, in terms of the p-value, which is given by p = ∞ q med f (q|0). The q obs from the data fit should be used to obtain the final statistical significance of the mantle signal, which will be described in Sec. 11.5.

RESULTS
This Section describes the results of our analysis. In Sec. 11.1 the final IBD candidates selected with the optimized selection cuts are presented. In Sec. 11.2 the analysis, and in particular the spectral fit with the 238 U/ 232 Th ratio fixed to the chondritic value (Sec. 11.2.1) or left free (Sec. 11.2.2), is described. The systematic uncertainties are discussed in Sec. 11.3. A summary of the geoneutrino signal as measured at the LNGS is given in Sec. 11.4. Considering the expected signal from the bulk lithosphere (Table VI), we estimate the geoneutrino signal from the mantle in Sec. 11.5. The consequences with regard to the Earth radiogenic heat are then presented in Sec. 11.6. Finally, in Sec. 11.7 the constraints on the power of a hypothetical georeactor (Sec. 5.5) are set.

Golden candidates
In the period between December 9, 2007 and April 28, 2019, corresponding to 3262.74 days of data acquisition, N IBD = 154 golden IBD candidates were observed to pass the data selection cuts described in Sec. 7.
The events are evenly distributed in time ( Fig. 46 (a)) and radially in the FV (Fig. 46 (b)). The charge distributions of the prompt and delayed signals are also compatible with the expectations, as shown in Figs. 46(c) and 46(d).
The distance to the inner vessel of the prompt signal was also studied. This test would be particularly sensitive to a potential background originated from the vessel itself or from the buffer: in the radial distribution of Fig. 46 (b), due to the changing IV shape, a small excess of this origin could have been smeared. As it is shown in Fig. 47, this test was done for all candidates, as well as separately for the geoneutrino energy window (below 1500 p.e.) and above. No excess was observed.

Analysis
An unbinned likelihood fit, as described in Sec. 10.1, was performed with the prompt charge of the 154 golden candidates shown in Sec. 11.1. 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 using the PDFs shown in Fig. 33 and Fig. 36(b), respectively. These components were constrained according to values in Table XV with Gaussian pull terms. Reactor antineutrinos were unconstrained in the fit, using the PDF as in the bottom left of Fig. 32. The differences in the shape of the reactor antineutrino spectra "with 5 MeV excess" and "without 5 MeV excess" (bottom right in Fig. 32) are included in the systematic uncertainty calculation (Sec. 11.3). Obviously, geoneutrinos were also kept unconstrained. The fit was performed in two different ways with respect to the relative ratio of the 232 Th and 238 U contributions, as detailed in the next two sub-sections.
The presented fit results are obtained following the recommendations given under the Statistics Chapter of [135] for cases, when there are physical boundaries on the possible parameter values. In our case, all the fit parameters must have non-negative values. For the main parameters resulting from the fit, the profiles of the likelihood L (Eq. 37) are provided, and in addition to the best fit values, the mean, median, as well as the 68% and 99.7% coverage intervals for non-negative parameter values, are provided in the summary Table XVII.

Th/U fixed to chondritic ratio
The fit was performed assuming the Th/U chondritic ratio and using the corresponding PDF shown in the top-left of Fig. 32. The resulting spectral fit is shown in Fig. 48(a). The likelihood profile for the number of geoneutrinos N geo (Fig. 49(a)), yields the best fit value N best geo = 51.9, the median value N med geo = 52.6, and the 68% coverage interval I 68stat Ngeo = 44.0 -62.0 events. For reactor antineutrinos, N best rea = 92.5 was obatined with the median value N med rea = 93.4, and the 68% coverage interval I 68stat Nrea = 82.6 -104.7 events. This is compatible with the reactor antineutrino expectation of 97.6 ± 1.7 events (without "5 MeV excess") as well as 91.9 ± 1.6 events (with "5 MeV excess"), given in Table XIV. The contour plot for N geo versus N rea is shown in Fig. 48(c). The fit was also performed by constraining the expected number of reactor antineutrino events to 97.6 ± 1.7 events (Table XIV) 61.1 events) is nearly unchanged with respect to that obtained when leaving the reactor antineutrino contribution free. The best fit value is shifted by about 1.5% and the error is only marginally reduced. This fit stability is due to the fact that above the geoneutrino energy window there is almost no non-antineutrino background, and thus the data above the geoneutrino end point well constrain the reactor antineutrino contribution also in the geoneutrino window. The fact that without any constraint on N rea the fit returns a value compatible with expectation is an important confirmation of the Borexino ability to measure electron antineutrinos.

Th and U as free fit parameters
The second type of fit was performed by treating 238 U and 232 Th contributions as free and independent fit components. The corresponding MC PDFs from topright of Fig. 32 were used. The spectral fit is shown in Fig. 48  tineutrinos, N best rea = 95.8 and I 68stat Nrea = 85.2 -109.0 events were obtained, which is also compatible with the expectation. The contour plot for N U versus N Th is shown in Fig. 48(d). The results obtained after constraining the expected N rea were again fully compatible with the results obtained when leaving the reactor antineutrino component free and without any significant reduction on error.

Systematic uncertainties
This section discusses the different sources of systematic uncertainty in the geoneutrino and reactor antineutrino measurement. They are detailed below and summarized in Table XVI.

Atmospheric neutrinos
Atmospheric neutrinos as the source of background were discussed in Sec. 5.4, while the expected number of IBD-like events passing the geoneutrino selection cuts in different energy regions was given in Table XIV. The uncertainty of this prediction is large, estimated to be 50%. However, there is an indication of some over-estimation of this background, since above the end-point of the reactor spectrum, where we would expect (3.3 ± 1.6) atmospheric events, no IBD candidates are observed. In the estimation of the systematic uncertainty due to atmospheric neutrinos, two fits were preformed which were similar to that shown in Fig. 48(a), but with additional contribution due to atmospheric neutrinos. These are represented by the PDF shown in the top part of Fig. 33. The fit was performed in two energy ranges and the number of  Table XVII. events from atmospheric neutrino background was constrained according to the values in Table XIV. Firstly, the fit was performed for data to 4000 p.e. and thus up the end point of the reactor antineutrino spectrum, the interval containing 63% of atmospheric neutrino background.
The resulting number of atmospheric neutrino events is 4.6 ± 3.2 and is compatible with the expectation of 6.7 ± 3.4. The geoneutrino signal is almost unchanged, while N rea decreased to (89.0 ± 11.3) events. Secondly, we have performed the fit up to the end-point of the atmospheric-neutrino background passing our IBD selection criteria (7500 p.e.). In this case, due to the fact that no IBD candidates are observed above the reactor antineutrino energy window, the resulting number of atmospheric neutrino background events is very low and with a large error, (1.2 ± 4.1) events. Fortunately, the resulting N geo and N rea are nearly unchanged. To summarize, we estimate the respective systematic uncertainty on geoneutrinos as +0.00 −0.38 % and on reactor antineutrinos as +0.0 −3.90 %.

Shape of the reactor spectrum
The likelihood fit, as described in Sec. 11.2 was performed using the MC PDF of reactor antineutrinos without any "5 MeV excess", based on the flux prediction of [105], as discussed in Sec. 5.3. In order to study the changes that might arise due to the observed "5 MeV excess", the fit was also performed using the corresponding MC PDF as shown in Fig. 32, based on the measured Daya Bay spectrum [101]. Since there is no constraint on N rea and the two spectral shapes are relatively similar, the change in N rea is very small: we observe an increase

Inner vessel shape reconstruction
We consider a conservative 5 cm error on the IV position (Sec. 3.3). This means that the function defining our DFV 10 cm inward from the IV is inside the scintillator with high probability. This implies that the systematic uncertainty on the FV defintion due to the IV shape reconstruction is negligible. However, there is a systematic uncertainty due to the selection of the IBD candidates using the DFV cut, which was evaluated by smearing the distance-to-IV of each IBD candidate with a Gaussian function with σ = 5 cm. Consequently, the DFV cut was applied on the smeared distances and the spectral fit was performed on newly selected candidates. This procedure was repeated 50 times. The distributions of the differences between the resulting N geo and N rea values with respect to the default fit have positive offsets, which were then conservatively taken as the systematic uncertainty due to the IV shape reconstruction. We estimate the respective systematic uncertainty on geoneutrinos as +3.46 −0.00 % and on reactor antineutrinos as +3.25 −0.00 %.

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. The trigger efficiency for the 2.2 MeV gamma from 241 Am-9 Be calibration source compared to MC simulations for different source positions was studied. The uncertainty in the efficiency was then set to a conservative limit of 1.5%.

Position reconstruction
The position of events in Borexino is calculated using the photon arrival times. 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. This uncertainty is obtained using the calibration campaign performed in 2009 [76]. Data from the 222 Rn and 241 Am- 9 Be sources placed at 182 and 29 positions in the scintillator, respectively, was used for this. The reconstructed position of the source was compared to the nominal source position measured by the CCD camera inside the detector. The uncertainty in position reconstruction for the geoneutrino analysis was calculated using the shift in the positions for the 241 Am-9 Be source. The maximal resulting uncertainty in the position was observed to be 5 cm. Considering the nominal spherical radius of our FV of 4.15 m, this gives an uncertainty of 3.6% in the fiducial volume and consequently, in the corresponding exposure.

Geoneutrino signal at LNGS
This Section details the conversion of the number of geoneutrino events N geo , resulting from the spectral fits described in Sec. 11.2, to the geoneutrino signal S geo expressed in TNU, the unit introduced in Sec. 5.2: where the detection efficiency ε geo = 0.8698 ± 0.0150 (Table XIII) Fig. 50. Table XVII summarizes the signals, expressed in TNU, for geoneutrinos and reactor antineutrinos obtained with the two fits, assuming Th/U chondritic ratio and keeping U and Th contributions as free fit parameters, as described in  Table VI), while the mantle signal is obtained considering an intermediate scenario ( Fig. 16(b)). The error bars represent the 1σ uncertainties of the total signal S (U+Th). The horizontal solid back line represents the geoneutrino signal S med geo , while the grey band the I 68full S geo interval as measured by Borexino.
Sec. 11.2. It was shown in Sec. 10.3 that Borexino does not have any sensitivity to measure the Th/U ratio with the current exposure. Therefore, the ratio obtained from the fit when U and Th are free parameters is not discussed.

Extraction of mantle signal
The mantle signal was extracted from the spectral fit by constraining the contribution from the bulk lithosphere according to the expectation discussed in Sec. 5.2 and given in Table XIV as (28.8 ± 5.6) events. The corresponding MC PDF was constructed from the PDFs of 232 Th and 238 U geoneutrinos shown in Fig. 32. They were scaled with the lithospheric Th/U signal ratio equal to 0.29 (Table VI). 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, the value discussed in Sec. 5.2. The mantle signal, as well as the reactor antineutrino contribution were free in the fit. The best fit is shown in Fig. 51(a). It resulted in a mantle signal of N best mantle = 23.1 events, with the median value N med mantle = 23.7 events, and the 68% coverage interval I 68stat Nmantle = (13.7 -34.4) events. The likelihood profile of the mantle signal is shown in Fig. 51(b).
After considering the systematic uncertainties, the final mantle signal can be given as S best mantle = 20.6 TNU, with the median value N med mantle = 21.2 TNU, and the 68% coverage interval I 68full S mantle = 12.2 -30.8 TNU, as shown also in Table XVII. The statistical significance of the mantle signal was studied using MC pseudo-experiments with and without a generated mantle signal as described in Sec 10.3. The q obs obtained from the spectral fit is 5.4479, and it is compared with the theoretical function f (q|0), described in Sec. 10.3, Eq. 39, as shown in Fig. 52. The corresponding p-value is 9.796 × 10 −3 . Therefore, in conclusion the null-hypothesis of the mantle signal can be rejected with 99.0% C.L. (corresponding to 2.3σ significance). The Borexino mantle signal can be compared with calculations according to a wide spectrum of BSE models (Table VII) suming for the mantle a Th/U mass ratio of 3.7.

Estimated radiogenic heat
The global HPEs' masses in the Earth are estimated by matching geophysical, geochemical, and cosmochemical arguments. Direct samplings of the accessible lithosphere constrain the radiogenic heat of H LS rad (U+Th+K) = 8.1 +1.9 −1.4 TW (Table V), corresponding to ∼17% of the total terrestrial heat power H tot = (47 ± 2) TW. The radiogenic heat from the unexplored mantle could embrace a wide range of H mantle rad (U+Th+K) = (1.2 -39.8) TW (Table VII), where the highest values are obtained for a Fully Radiogenic Earth model.
The total amount of HPEs, as well as their distribution in the deep Earth, affect the geoneutrino flux. We will express the dependence of the expected mantle geoneutrino signal S mantle (U+Th) on the mantle radiogenic heat H mantle rad (U+Th). Since the radiogenic heat is unequivocally related with the masses of HPEs, this relation is first expressed through M mantle (U), the U mass in the mantle (Table VII), assuming the mantle Th/U mass ratio of 3.7. The dependence on the HPEs' distribution is taken into account following the approach presented in Sec. 5.2. Assuming the homogeneous mantle (Fig. 16c), we express, in units of TNU, the high scenario signal S HS mantle : S HS mantle (U + Th) = 0.98·[h(U) + 3.7 · h(Th)]· M mantle (U), (42) while assuming a unique HPE-rich layer just above the CMB (Fig. 16a), we express the low scenario signal FIG. 53. Mantle geoneutrino signal expected in Borexino as a function of U and Th mantle radiogenic heat: the area between the red (Eq. 44) and blue (Eq. 45) lines denotes the full range allowed between a homogeneous mantle (high scenario - Fig. 16c) and a unique rich layer just above the CMB (low scenario - Fig. 16a). The average slope of this dependence (black inclined line) corresponds to 1/1.16, see Eq. 46. The blue, green, red, and yellow ellipses represent different BSE models (CC, GC, GD, and FR; Table VII, respectively), while 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).  (Table VII), and the specific heat h(U) = 98.5 µW/kg and h(Th) = 26.3 µW/kg are taken from [26]. From the above equations follow direct relations between the signals in TNU and radiogenic power in TW, according to the two mantle scenarios: S HS mantle (U + Th) = 0.98 · H mantle rad (U + Th), S LS mantle (U + Th) = 0.75 · H mantle rad (U + Th).
These relations are depicted by the red and the blue lines of Fig. 53, respectively, and the area between them denotes the region allowed by all possible U and Th distributions in the mantle, assuming that the abundances in this reservoir are radial and non-decreasing function of the depth. The maximal and minimal excursions of mantle geoneutrino signal is taken as a proxy for the ±3σ error range.
Since the radiogenic heat power of the lithosphere is independent from the BSE model, the discrimination capability of Borexino geoneutrino measurement  (Table VII), while for Borexino the value of 30.0 +13.5 −12.7 TW is inferred from the extracted mantle signal. The difference between H tot and the respective total radiogenic heat is assigned to the heat from secular cooling of the Earth. among the different BSE models can be studied in the space S mantle (U + Th) vs H mantle rad (U + Th). In Figure 53, the solid black horizontal line represents the Borexino measurement, the median S med mantle , which falls within prediction of the Geodynamical model (GD). The 68% coverage interval I 68full S mantle , also represented in this Figure by horizontal black dashed lines, covers the area of prediction of the GD and the Fully Radiogenic (FR) models. We are least compatible with the Cosmochemical model (CC), at 2.4σ with its central value.
The mantle signal measured by Borexino can be converted to the corresponding radiogenic heat by inverting the Eqs. 44 and 45. Since the experimental error on the mantle signal is much larger than the systematic variability associated to the U and Th distribution in the mantle, the average of the respective slopes, in particular its inverse β = 1.16 is used: H mantle rad (U + Th) = β · S mantle (U + Th).
We obtain: Assuming the contribution from 40 K to be 18% of the total mantle radiogenic heat (Sec. 2), the total radiogenic  mantle signal can be expressed as H mantle rad (U + Th + K) = 30.0 +13.5 −12.7 TW, where we have expressed the 1σ errors with respect to the median. If we further add the lithospheric contribution H LS rad (U+Th+K) = 8.1 +1.9 −1.4 TW, we get the 68% coverage interval for the Earth's radiogenice heat H rad (U + Th + K) = 38.2 +13.6 −12.7 TW, as shown in Fig. 54.
The experimental error on the Earth's radiogenic heat power estimated by Borexino is comparable with the spread of power predictions derived from the eight BSE models reported in Table II. This comparison is represented in Fig. 54. Among these, a preference is found for models with relatively high radiogenic power, which correspond to a cool initial environment at early Earth's formation stages and small values of the current heat coming from the secular cooling. 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 Eq. 6. The resulting value of UR CV = 0.78 +0. 41 −0.28 is compared to the UR CV predicted by different BSE models in Fig. 55. The Borexino geoneutrino measurement constrains at 90(95)% C.L. a mantle radiogenic heat power 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 UR CV > 0.13(0.04).

Testing the georeactor hypothesis
The georeactor hypothesis described in Sec. 5.5 was tested by performing the spectral fit after constraining the expected number of reactor antineutrino events (Table XIV) to 97.6 ± 1.7 (stat) ± 5.2 (syst). The geoneu-  (Fig. 34), corresponding to different source positions in the deep Earth ( Fig. 21(a)), in the spectral fit after constraining the reactor antineutrinos. The vertical dashed lines correspond to the upper limit at 95% C.L. for different positions. trino (Th/U fixed to chondritic mass ratio of 3.9) and georeactor contributions were left free in the fit. For each georeactor location ( Fig. 21(a)), we have used the respective georeactor PDF as in Fig. 34. However, their shapes are practically identical and Borexino does not have any sensitivity to distinguish them. The different likelihood profiles obtained using different georeactor PDFs are very similar (including the PDF constructed assuming no neutrino oscillations), as shown in Fig. 56. The vertical lines represent the 95% C.L. limits for the number of georeactor events obtained with different PDFs. Conservatively, the highest limit is transformed to the signal of 18.7 TNU, as the 95% C.L. upper limit on the signal coming from a hypothetical georeactor. Considering the values from Table VIII, that is, the predicted georeactor signal expressed in TNU for a 1 TW georeactor in different locations, these upper limits on the georeactor power are set: 2.4 TW for the location in the Earth center (GR2) and 0.5 TW and 5.7 TW for the georeactor placed at the CMB at 2900 km (GR1) and 9842 km (GR3), respectively.

CONCLUSIONS
Borexino is 280-ton liquid scintillator neutrino detector located at Laboratori Nazionali del Gran Sasso (LNGS) in Italy, and has been acquiring data since 2007. It has proven to be a successful neutrino observatory, which went well beyond its original proposal to observe 7 Be solar neutrinos. Radiopurity of the detector, its calibration, stable performance, its relatively large distance to nuclear reactors, as well as depth of the LNGS laboratory to guarantee smallness of cosmogenic background, are the main building blocks for a geoneutrino measurement with systematic uncertainty below 5%.
The paper provides an in-depth motivation and description of geoneutrino measurement, as well as the geological interpretations of the result. It presents in detail the analysis of 3262.74 days of Borexino data taken between December 2007 and April 2019 and provides, with some assumptions, a measurement of the Uranium and Thorium content of the Earth's mantle and its radiogenic heat.
Borexino detects geoneutrinos from 238 U and 232 Th through inverse beta decay, in which electron flavour antineutrinos with energies above 1.8 MeV interact with protons in the detector. The detection efficiency for optimized data selection cuts is (86.98 ± 1.50)%. This interaction is the only channel presently available for detection of MeV-scale electron antineutrinos. Optimized data selection including an enlarged fiducial volume and a sophisticated cosmogenic veto resulted in an exposure of (1.29 ± 0.05) × 10 32 protons × year. This represents an increase by a factor of two over the previous Borexino analysis reported in 2015. The focus of the paper is to provide the scientific community with a comprehensive study that combines the expertise of neutrino physicists and geoscientists.
The paper documents improved techniques in the in-depth analysis of the Borexino data, and provides future experiments with a description of the substantial effort required to extract geoneutrino signals. We have underlined the importance of muon detection (in particular special categories of muon events that become crucial in low-rate measurements), as well as the α/β pulse shape discrimination techniques. The optimization of data selection cuts, chosen to maximize Borexino's sensitivity to measure geoneutrinos, has been described.
All kinds of background types considered important for geoneutrino measurement have also been discussed, including approaches of their estimation either through theoretical calculation and Monte Carlo simulation, or by analysis of independent data.
Borexino ability to measure electron antineutrinos is calibrated via reactor antineutrino background, that is not constrained in geoneutrino analysis and has been found to be in agreement with the expectations. By observing 52.6 +9.4 −8.6 (stat) +2.7 −2.1 (sys) geoneutrinos (68% interval) from 238 U and 232 Th, a geoneutrino signal of 47.0 +8.4 −7.7 (stat) +2.4 −1.9 (sys) TNU has been obtained. The total precision of +18.3 −17.2 % is found to be in agreement with the expected sensitivity. This result assumes a Th/U mass ratio of 3.9, as found in chondritic CI meteorites, and is compatible with result when contributions from 238 U and 232 Th were both fit as free parameters.
Importance of the knowledge of abundances and distributions of U and Th in the Earth, and in particular around the detector, for both the signal prediction as well as interpretation of results, have been discussed. The measured geoneutrino signal is found to be in agreement with the predictions of different geological models with a preference for those predicting higher signals. The hypothesis of observing a null mantle signal has been excluded at 99% C.L. when exploiting detailed knowledge of the local crust near the LNGS. The latter is characterized by the presence of U and Th depleted sediments. We note that geophysical and geochemical observations constrain the Th/U mass ratio for the bulk lihtosphere to a value of 4.3.
Maintaining the global chondritic ratio of 3.9 for the bulk Earth, the inferred Th/U mass ratio for the mantle is 3.7. Assuming the latter value, we have observed mantle signal of 21.2 +9.6 −9.0 (stat) +1.1 −0.9 (sys) TNU. Considering different scenarios about the U and Th distribution in the mantle, the measured mantle geoneutrino signal has been converted to radiogenic heat from U and Th in the mantle of 24.6 +11.1 −10.4 TW (68% interval). Assuming the contribution of 18% from 40 K in the mantle and adding the relatively-well known lithospheric radiogenic heat of 8.1 +1.9 −1.4 TW, Borexino has estimated the total radiogenic heat of the Earth to be 38.2 +13.6 −12.7 TW. The latter is found to be compatible with different geological predictions. However, there is a ∼2.4σ tension with Earth models predicting the lowest concentration of heat-producing elements. The total radiogenic heat estimated by Borexino can be used to extract a convective Urey ratio of 0.78 +0. 41 −0.28 . In conclusion, Borexino geoneutrino measurement has constrained at 90% C.L. the mantle composition to a mantle (U) > 13 ppb and a mantle (Th) > 48 ppb, the mantle radiogenic heat power to H mantle rad (U+Th) > 10 TW and H mantle rad (U+Th+K) > 12.2 TW, as well as the convective Urey ratio to UR CV > 0. 13. With the application of a constraint on the number of expected reactor antineutrino events, Borexino has placed an upper limit on the number of events from a hypothetical georeactor inside the Earth. Assuming the georeactor located at the center of the Earth, its existence with a power greater than 2.4 TW has been excluded at 95% C.L.
In conclusion, Borexino confirms the feasibility of geoneutrino measurements as well as the validity of different geological models predicting the U and Th abundances in the Earth. This is an enormous success of both neutrino physics and geosciences. However, in spite of some preference of Borexino results for the models predicting high U and Th abundances, additional and more precise measurements are needed in order to extract firm geological results. The next generation of large volume liquid scintillator detectors has a strong potential to provide fundamental information about our planet.