Comprehensive geoneutrino analysis with Borexino M. Agostini,16 K. Altenmler,16 S. Appel,16 V. Atroshchenko,6 Z. Bagdasarian,22 D. Basilico,9 G. Bellini,9 J. Benziger,13 D. Bick,26 G. Bonfini,8 D. Bravo,9,c B. Caccianiga,9 F. Calaprice,12 A. Caminata,3 L. Cappelli,8 P. Cavalcante,15,e F. Cavanna,3 A. Chepurnov,17 K. Choi,21 D. D’Angelo,9 S. Davini,3 A. Derbin,11 A. Di Giacinto,8 V. Di Marcello,8 X. F. Ding,18,8,12 A. Di Ludovico,12 L. Di Noto,3 I. Drachnev,11 G. Fiorentini,28,29 A. Formozov,2,9,17 D. Franco,1 F. Gabriele,8 C. Galbiati,12 M. Gschwender,24 C. Ghiano,8 M. Giammarchi,9 A. Goretti,12,e M. Gromov,17,2 D. Guffanti,18,8,f C. Hagner,26 E. Hungerford,25 Aldo Ianni,8 Andrea Ianni,12 A. Jany,4 D. Jeschke,16 S. Kumaran,22,23 V. Kobychev,5 G. Korga,25,g T. Lachenmaier,24 T. Lasserre,27 M. Laubenstein,8 E. Litvinovich,6,7 P. Lombardi,9 I. Lomskaya,11 L. Ludhova ,22,23 G. Lukyanchenko,6 L. Lukyanchenko,6 I. Machulin,6,7 F. Mantovani,28,29 G. Manuzio,3 S. Marcocci,18,†,d J. Maricic,21 J. Martyn,20,f E. Meroni,9 M. Meyer,19 L. Miramonti,9 M. Misiaszek,4 M. Montuschi,28,29 V. Muratova,11 B. Neumair,16 M. Nieslony,20,f L. Oberauer,16 A. Onillon,27 V. Orekhov,20,f F. Ortica,10 M. Pallavicini,3 L. Papp,16 Ö. Penek,22,23 L. Pietrofaccia,12 N. Pilipenko,11 A. Pocar,14 G. Raikov,6 M. T. Ranalli,8 G. Ranucci,9 A. Razeto,8 A. Re,9 M. Redchuk,22,23 B. Ricci,28,29 A. Romani,10 N. Rossi,8,a S. Rottenanger,24 S. Schert,16 D. Semenov,11 M. Skorokhvatov,6,7 O. Smirnov,2 A. Sotnikov,2 V. Strati,28,29 Y. Suvorov,8,6,b R. Tartaglia,8 G. Testera,3 J. Thurn,19 E. Unzhakov,11 A. Vishneva,2 M. Vivier,27 R. B. Vogelaar,15 F. von Feilitzsch,16 M. Wojcik,4 M. Wurm,20,f O. Zaimidoroga,2,† S. Zavatarelli,3 K. Zuber,19 and G. Zuzel4 (Borexino Collaboration)* 1AstroParticule et Cosmologie, Universit´e Paris Diderot, CNRS/IN2P3, CEA/IRFU, Observatoire de Paris, Sorbonne Paris Cit´e, 75205 Paris Cedex 13, France 2Joint Institute for Nuclear Research, 141980 Dubna, Russia 3Dipartimento di Fisica, Universit`a degli Studi e INFN, 16146 Genova, Italy 4M. Smoluchowski Institute of Physics, Jagiellonian University, 30348 Krakow, Poland 5Kiev Institute for Nuclear Research, 03680 Kiev, Ukraine 6National Research Centre Kurchatov Institute, 123182 Moscow, Russia 7National Research Nuclear University MEPhI (Moscow Engineering Physics Institute), 115409 Moscow, Russia 8INFN Laboratori Nazionali del Gran Sasso, 67010 Assergi (AQ), Italy 9Dipartimento di Fisica, Universit`a degli Studi e INFN, 20133 Milano, Italy 10Dipartimento di Chimica, Biologia e Biotecnologie, Universit`a degli Studi e INFN, 06123 Perugia, Italy 11St. Petersburg Nuclear Physics Institute NRC Kurchatov Institute, 188350 Gatchina, Russia 12Physics Department, Princeton University, Princeton, New Jersey 08544, USA 13Chemical Engineering Department, Princeton University, Princeton, New Jersey 08544, USA 14Amherst Center for Fundamental Interactions and Physics Department, University of Massachusetts, Amherst, Massachusetts 01003, USA 15Physics Department, Virginia Polytechnic Institute and State University, Blacksburg, Virginia 24061, USA 16Physik-Department and Excellence Cluster Universe, Technische Universität Mchen, 85748 Garching, Germany 17Lomonosov Moscow State University Skobeltsyn Institute of Nuclear Physics, 119234 Moscow, Russia 18Gran Sasso Science Institute, 67100 L’Aquila, Italy 19Department of Physics, Technische Universität Dresden, 01062 Dresden, Germany 20Institute of Physics and Excellence Cluster PRISMA+, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany 21Department of Physics and Astronomy, University of Hawaii, Honolulu, Hawaii 96822, USA 22Institut fr Kernphysik, Forschungszentrum Jich, 52425 Jich, Germany 23RWTH Aachen University, 52062 Aachen, Germany 24Kepler Center for Astro and Particle Physics, Universität Tingen, 72076 Tingen, Germany 25Department of Physics, University of Houston, Houston, Texas 77204, USA 26Institut fr Experimentalphysik, Universität Hamburg, 22761 Hamburg, Germany 27IRFU, CEA, Universit´e Paris-Saclay, F-91191 Gif-sur-Yvette, France M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 28Dipartimento di Fisica e Scienze della Terra, Universit`a di Ferrara, Via Saragat 1, I-44122 Ferrara, Italy 29INFN—Sezione di Ferrara, Via Saragat 1, I-44122 Ferrara, Italy (Received 5 September 2019; published 21 January 2020; corrected 25 February 2020) 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, 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 .stat..2.7 .sys. geoneutrinos (68% interval) from 238Uand -8.6 -2.1 232Th, a geoneutrino signal of 47.0.8.4 .stat..2.4 .sys. TNU with .18.3% total precision was obtained. This -7.7 -1.9 -17.2 result assumes the same Th/U mass ratio as found inchondritic CI meteorites but compatible results were found when contributions from 238Uand 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 isexcluded ata 99.0% C.L. when exploiting detailed knowledgeof the local crust near the experimental site. Measured mantle signal of 21.2.9.5 .stat..1.1 .sys. TNU corresponds to the -9.0 -0.9 production of a radiogenic heat of 24.6.11.1 TW (68% interval) from 238Uand 232Th in the mantle. Assuming -10.4 18% contribution of 40Kin themantleand 8.1.1.9 TW of total radiogenic heat of the lithosphere, the Borexino -1.4 estimate of the total radiogenic heat of the Earth is 38.2.13.6 TW, which corresponds to the convective Urey -12.7 ratio of 0.78.0.41 . These values are compatible with different geological predictions, however there is a ~2.4. -0.28 tension with those Earth models which predict the lowest concentration of heat-producing elements in the mantle. In addition, by constraining the number of expected reactor antineutrino events, the existence of a hypothetical georeactor at the center of the Earth having power greater than 2.4 TW is excluded at 95% C.L. Particular attention is given to the description of all analysis details which should be of interest for the next generation of geoneutrino measurements using liquid scintillator detectors. DOI: 10.1103/PhysRevD.101.012009 I. INTRODUCTION Neutrinos, the most abundant massive particles in the universe, are produced by a multitude of different *spokesperson-borex@lngs.infn.it †Deceased. aPresent address: Dipartimento di Fisica, Sapienza Universit`a di Roma e INFN, 00185 Roma, Italy.bPresent address: Dipartimento di Fisica, Universit`a degli Studi Federico II e INFN, 80126 Napoli, Italy. cPresent address: Universidad Automa de Madrid, Ciudad Universitaria de Cantoblanco, 28049 Madrid, Spain.dPresent address: Fermilab National Accelerator Laboratory (FNAL), Batavia, Illinois 60510, USA. ePresent address: INFN Laboratori Nazionali del Gran Sasso, 67010 Assergi (AQ), Italy.fInstitute of Physics and Excellence Cluster PRISMA+, Jo­hannes Gutenberg-Universität Mainz, 55099 Mainz, Germany. gAlso at MTA-Wigner Research Centre for Physics, Depart­ment of Space Physics and Space Technology, Konkoly-Thege Mikls t 29-33, 1121 Budapest, Hungary. Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. processes. They interact only by the weak and gravitational interactions, and so are able to penetrate enormous dis­tances 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 atmospheric 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 us to assess the Earth’s heat budget, specifically the heat emitted in the radioactive decays. The latter, the so-called radiogenic heat of the present Earth, arises mainly from the decays of isotopes with half-lives comparable to, or longer than Earth’s age (4.543 × 109 years): 232Th (T1=2 1 1.40×1010 years), 238U (T1=2 1 4.468×109 years), 235U(T1=2 1 7.040×108 years), and 40K(T1=2 1 1.248 × 109 years) [4]. All these isotopes are labeled as heat-producing elements (HPEs). The natural Thorium is fully composed of 232Th, while the natural isotopic abundances of 238U, 235U, and 40K are 0.992742, 0.007204, and 1.17 × 10-4, respectively. In each decay, the emitted radiogenic heat is in a well-known ratio1 to the number of emitted geoneutrinos [5]: 238U › 206Pb . 8. . 8e- . 6—.e . 51.7 MeV .1. 235U › 207Pb . 7. . 4e- . 4.—e . 46.4 MeV .2. 232Th › 208Pb . 6. . 4e- . 4—.e . 42.7 MeV .3. 40K › 40Ca . e- . .—e . 1.31 MeV .89.3%..4. 40K . e- › 40Ar . .e . 1.505 MeV .10.7%..5. Obviously, the total amount of emitted geoneutrinos scales with the total mass of HPEs inside the Earth. Hence, geoneutrinos’ detection provides us a way of measuring this radiogenic heat. This idea was first discussed by G. Marx and N. Menyhárd [6], G. Eder [7], and G. Marx [8] in the 1960s. It was further developed by M. L. Krauss, S. L. Glashow, and D. N. Schramm [9] in 1984. Finally, the potential to measure geoneutrinos with liquid scintillator detectors was suggested in the 1990s by C. G. Rothschild, M. C. Chen, and F. P. Calaprice [10] and independently by R. Raghavan et al. [11]. It took several decades to prove these ideas feasible. Currently, large-volume liquid-scintillator neutrino experi­ments KamLAND [12–15] and Borexino [16–18] have demonstrated the capability to efficiently detect a geo­neutrino 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 7Be solar neutrinos. However thanks 1The energy expressed in the following equations is the total energy, from which the released geoneutrinos take away about 5% as their kinetic energy. to the unprecedented levels of radiopurity, Borexino has surpassed its original goal and has now measured all2 the pp-chain neutrinos [20–22]. We report here a comprehen­sive 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. × 1032 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 continuous develop­ment of analysis techniques. This paper is structured as follows: Section II introduces the fundamental insights on what the geoneutrino studies can bring to the comprehension of the Earth’s inner structure and thermal budget. Section III details a descrip­tion of the Borexino detector and the structure of its data. In Sec. IV, the .—e detection reaction—the inverse beta decay on 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 40K and 235U geoneutrinos completely unreachable with present-day detection techniques. Section V deals with the estimation of the expected antineutrino signal from geoneutrinos, through background from reactor and atmos­pheric neutrinos, up to a hypothetical natural georeactor in the deep Earth. Section VI describes the nonantineutrino 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. VII, which involves the optimization of the signal-to­background ratio. Section VIII 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 IX introduces the analyzed dataset and discusses the number of expected signal and background events passing the optimized cuts, based on Secs. V and VI.In Sec. X, the Borexino sensitivity to extract geoneutrino 2The upper limit was placed for hep solar neutrinos, the flux of which is expected to be about 3 orders of magnitude smaller than that of 8B solar neutrinos [19]. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) signals is illustrated. Finally, Sec. XI discusses our results. The golden IBD candidate sample is presented (Sec. XI A) together with the spectral analysis (Sec. XI B) and sources of systematic uncertainty (Sec. XI C). The measured geo­neutrino signal at LNGS is compared to the expectations of different geological models in Sec. XI D. The extraction of the mantle signal using knowledge of the signal from the bulk lithosphere is discussed in Sec. XI E. The conse­quences of the new geoneutrino measurement with respect to the Earth’s radiogenic heat are discussed in Sec. XI F. Placing limits on the power of a hypothetical natural georeactor, located at different positions inside the Earth, is discussed in Sec. XI G. Final summary and conclusions are reported in Sec. XI. The acronyms used within the text are listed in alphabetical order in the Appendix. II. WHY STUDY GEONEUTRINOS? Our Earth is unique among the terrestrial planets3 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, geodynam­ical, 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 of differentiation, i.e., a transformation from a homo­geneous 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 ME 1 5.97 × 1024 kg, is relatively well established from its density profile, which is obtained by precise measure­ments of seismic waves on its surface. During the first differentiation, metallic segregation occurred and the core (~32.3% of ME) separated from the silicate primitive mantle or bulk silicate earth (BSE). The latter further differentiated into the present mantle (~67.2% of ME) and crust (~0.5% of ME). 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 3Mercury, Venus, Earth, and Mars. FIG. 1. Schematic cross section of the Earth. 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 tomogra­phy 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 midoceanic ridges where the oceanic crust is formed. The continental crust is more complex and thicker than the oceanic crust. 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 discontinuity4 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 4The Moho (Mohorovičić) discontinuity is the boundary between the crust and the mantle, characterized by a jump in seismic compressional waves velocities from ~7 to ~8 km=s occurring beneath the CC at typical depth of ~35 km. history being the most differentiated and heterogeneous layer. It consists of igneous, metamorphic, and sedimentary rocks. The OC with .83.km thickness [35] is created along the midoceanic ridges, where the basaltic magma differentiates from the partially melting mantle up-welling toward 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 cm2 at MeV energy for the IBD, Sec. IV), 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 composition of the Earth’s interior [5]. Consequently, it also provides an insight into the radiogenic heat contribution to the mea­sured Earth’s surface heat flux. The heat flow from the Earth’s surface to space results from a large temperature gradient across the Earth [37]. Table I shows estimations of this heat flow, Htot, integrated over the whole Earth’s surface. Different studies of this flux are based on several thousands of inhomogeneously dis­tributed 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 midocean ridges (where the data are sparse), requires the application of energy-loss models [38]. Except for Ref. [39], the papers account for the hydrothermal circu­lation 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 midocean TABLE I. Integrated terrestrial surface heat fluxes Htot esti­mated by different authors. The lower limit estimation [39] 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) [43] Davies (1980) [44] Sclater et al. (1980) [45] Pollack et al. (1993) [46] Hofmeister et al. (2005) [39] Jaupart et al. (2007) [41] Davies & Davies (2010) [42] 43 41 42 44 1 31 1 46 3 47 2 ridges, where the oceanic crust is created [40] and from where the older crust is pushed away by a newly created crust. This approach leads to an Htot estimation between 41 and 47 TW, with the oceans releasing ~70% of the total escaping heat. The most recent models [41,42] are in excellent agreement and provide a value of (46–47) TW with (2–3) TW error. However, Ref. [39], based only on direct measurements and not applying the half-space cool­ing model, provides a much lower Htot of .31 1.TW. We assume a Htot 1.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 Htot is typically expected to originate from two main processes: (i) secular cooling HSC of the Earth, i.e., cooling from the time of the Earth’s formation when gravitational binding energy was released due to matter accretion, and (ii) radiogenic heat Hrad from HPEs’ radioactive decays in the Earth. The relative contribution of radiogenic heat to the Htot 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 (URCV) quantifies the ratio of internal heat generation in the mantle over the mantle heat flux, as the following ratio [37]: Hrad - HCC URCV 1 rad ; .6. Htot - HCC rad where HCC 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 [38], 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 URCV must be in the range between 0.12 and 0.49 [38]. Additionally, HPEs’ abundances, and thus Hrad 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 abun­dances (e.g., O, Si, Mg, Fe) within 10% [47]. Uranium and Thorium are refractory (condensate at high temperatures) and lithophile (preferring to bind with silicates over metals) elements. The relative abundances of the refractory lith­ophile elements are expected to be stable to volatile loss or core formation during the early stage of the Earth [48]. The content of refractory lithophile elements (e.g., U and Th), M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) TABLE II. Masses M and abundances a of HPEs in the bulk silicate earth (MBSE 14.04 × 1024 kg [32]) predicted by different models: J: Javoy et al., 2010 [34], L & K: Lyubetskaya & Korenaga, 2007 [52], T: Taylor, 1980 [53], M & S: McDonough & Sun, 1995 [47], A: Anderson, 2007 [54], W: Wang et al., 2018 [55], P & O: Palme and O’Neil, 2003 [56], T & S: Turcotte & Schubert, 2002 [57]. The cosmochemical (CC), geochemical (GC), and geodynamical (GD) BSE models correspond to the estimates reported in [58]; the fully radiogenic (FR) model is defined adopting the approach of [59], assuming that the total heat Htot 1.47 2.TW is due to only the radiogenic heat production Hrad.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 [57], respectively. The GC is based on estimates reported by McDonough & Sun, 1995 [47] with K abundances corrected following [60]. The radiogenic heat Hrad released in the radioactive decays of HPEs is calculated adopting the element specific heat generation h (Hrad 1h × M with h.U.198.5 µW=kg, h.Th.126.3 µW=kg, and h.K.1 3.33 × 10-3 µW=kg taken from [59]). It is assumed that the uncertainties on U, Th, and K abundances are fully correlated. Model a.U. 1ng=g a.Th. 1ng=g a.K. 1µg=g M.U. 11016 kg M.Th. 11016 kg M.K. 11019 kg Hrad.U. [TW] Hrad.Th. [TW] Hrad.K. [TW] Hrad.U .Th .K. [TW] J 12 43 146 4.85 17.4 59 4.8 4.6 2.0 11.3 L & K 17 63 190 6.87 25.5 76.8 6.8 6.7 2.6 16 T 18 70 180 7.28 28.3 72.8 7.2 7.5 2.4 17 M & S 20 80 240 8.09 32.4 97.1 8.0 8.5 3.2 19.7 A 20 77 151 8.09 31.1 61.1 8.0 8.2 2.0 18.2 W 20 75 237 8.09 30.3 95.8 8.0 8.0 3.2 19.1 P & O 22 83 260 8.9 33.6 105.1 8.8 8.9 3.5 21.1 T & S 35 140 350 14.2 56.6 141.5 13.9 14.9 4.7 33.5 CC 12 2 43 4 146 29 5 1 17 2 59 12 4.8 0.84.6 0.42.0 0.4 11.3 1.6 GC 20 480 13 280 60 8 2 32 5 113 24 8.0 1.68.5 1.43.8 0.8 20.2 3.8 GD 35 4 140 14 350 35 14 2 57 6 142 14 13.9 1.614.9 1.54.7 0.5 33.5 3.6 FR 49 2 189 8 554 24 20 1 77 3 224 10 19.4 0.820.2 0.87.5 0.3 47 2 which are excluded from the core,5 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 compositional models. The contributions to the radiogenic heat of U, Th, and K vary in the range of [39–44]%, [40–45]%, and [11–17]%, respectively. Three classes of BSE models are adopted in this work: the cosmochemical, geochemical, and geodynamical mod­els, as defined in [32,58]. The cosmochemical (CC) model [34] is characterized by a relatively low amount of U and Th producing a total Hrad 1.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 con­strained by terrestrial samples [47,60]. 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 [57]. Additionally, an extreme 5Recent speculations [49] about possible partitioning of some lithophile elements (including U and Th) into the metallic core are still debated [50,51]. 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. model can be obtained following the approach described in [59], where the terrestrial heat Htot of 47 TW is assumed to be fully accounted for by radiogenic production Hrad. 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 MTh=MU 1.3.876 0.016.[61]. Recent stud­ies [62], based on measured molar 232Th=238U values and their time integrated Pb isotopic values, are in agreement estimating MTh=MU 13.90.0.13 Significant deviations -0.08 . from this average value can be found locally, especially in the heterogeneous continental crust. This fact is attrib­utable 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. VB), the reservoirs of the sedimentary cover, which account for 30% of the geo­neutrino 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 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 lithospheric contribution, as we show in Sec. XI E. 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. XI G) can be set on the power of a hypothetical Uranium natural georeactor suggested in [66–69] and discussed in Sec. VE. In future, by combining measurements from several experi­ments placed in distant locations and in distinct geological environments, one could test whether the mantle is laterally homogeneous or not [58], as suggested, for example, by the large shear velocity provinces observed at the mantle base below Africa and Pacific ocean [70]. 40K In future, detection of geoneutrinos might be possible [71,72]. This would be extremely important, since Potassium is the only semivolatile HPE. Our planet seems to show ~1=3 [47] 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 [47] 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 16 000 [58]. 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. XI F). III. 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. (meter water equivalent). The rock above the detector provides shielding against cosmogenic backgrounds such that the muon flux is decreased to -2 .3.432 0.003. × 10-4 ms-1 [75]. The general scheme of the Borexino detector is shown in Fig. 2. The detector has a concentric multilayer structure. The outer layer [outer detector (OD)] serves as a passive shield against external radiation as well as an active Cherenkov veto of cosmo­genic 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 ultrapure 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, C6H3.CH3.3. solvent doped with a fluorescent dye PPO (2,5­diphenyloxazole, C15H11NO) in concentration of 1.5 g=l. The scintillator density is .0.878 0.004. gcm-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. × 1028 per 1 ton. A careful selection of detector materials, accurate assem­bling, and a complex radio-purification of the liquid scintillator guaranteed extremely low contamination levels of 238U and 232Th. After the additional LS purifica­tion 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 dime­thylphthalate (DMP, C6H4.COOCH3.2) light quencher in PC, shields the core of the detector against external .s M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 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=lto 3 g=land 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. III C. Borexino has a main data acquisition system (the main DAQ) and a semi-independent fast wave form digitizer (FWFD) or flash analog-to-digital converter (FADC) subsystem 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).6 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. III A. The effective light yield in Borexino is approximately 500 detected photoelectrons (p.e.) per 1 MeV of deposited p................. energy. This results in the 5%=E.MeV. energy reso­lution. Borexino is a position sensitive detector. For pointlike 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 per­formed 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 inter-action 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 and 1.46 MeV were used for studying the energy scale. 222Rn source, emitting alpha particles characterized by pointlike interactions, was applied to study the homogeneity of the detector’s response 6The boards were designed and built in collaboration with Laben s.p.a. as well as position reconstruction. For geoneutrino studies, 241Am–9Be employment of the source is of particular interest, since the emitted neutrons closely represent the delayed signal of an inverse beta decay (Sec. IV), which is the interaction used to detect geoneutrinos. In addition, a 228Th 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 per­formed once a week with a special laser (. 1 394 nm, 50 ps wide peak) calibration run. This procedure is of utmost importance for position and muon track reconstruc­tions, as well as for the .=ß discrimination (Sec. III D), based on their different fluorescence time profiles. The charge calibration of the single photoelectron response of each PMT is performed typically 4 times a day. It is based on plentiful ß- decays of 14C(Q 1 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. Anew 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 perfor­mance and further improved stability of the detector. A. Borexino data structure Borexino electronics must handle about 106 events per day, which are dominated by 14C decays. The residual flux of cosmogenic muons results in the detection of approx­imately 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. 1. 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 subdetectors 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 22 bit of the trigger word, and is referenced as the BTB4 condition,or muon trigger flag (MTF). This is described in Sec. III B. 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 (TT), 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. (i) TT1 & BTB0—pointlike events—the main physics trigger type for ID events, 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. (ii) TT1 & BTB4—internal muons—the main category of internal muons which did trigger both the OD as well as the ID. (iii) TT128—neutron trigger—special 1.6 ms long trig­ger 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]. (iv) TT2—external muons—the main category of exter­nal muons detected by the OD only. (v) TT8—laser—calibration trigger of the timing laser used to monitor the quality of the laser pulse. (vi) TT32—pulser—calibration trigger with the elec­tronics pulse, used to monitor the number of work­ing electronic channels. (vii) TT64—random—forced trigger, used to monitor the dark noise and the 14C-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, i.e., the hits that are not accompanied by possible error messages from the Laben digital boards. 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: (i) NP —number of triggered PMTs; (ii) Nh —number of hits within the cluster, including possible multiple hits from the same PMT; (iii) Npe —number of photoelectrons, calculated as a sum of charges of all individual hits contributing to Nh. 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. 2. The FADC DAQ subsystem There is an auxiliary data acquisition system based on 34 FWFD (also known as 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. III B), 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 pointlike event, a muon, and a calibration pulser signal. As it will be discussed in Sec. III B, 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 wave­forms are extended up to 16 µs (or 1.6 ms for TT128 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. B. 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 methods in Borexino, and also provide updated analysis of the overall muon tagging efficiency using the main DAQ and evaluate its stability M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 104 3 10 102 10 1 104 3 10 102 10 1 104 103 102 10 1 3 10 102 10 1 106 5 10 104 3 10 102 10 1 102 10 1 -15000 -10000 -5000 0 ID hit time [ns] OD hit time [ns] FIG. 3. Event structure, i.e., hit time distribution in a 16 µs DAQ gate with respect to the trigger time, for different event types as acquired during a typical 6 hour run by the ID (left column) and the OD (right column). The hit times are negative, since the hits are detected before the trigger decision is made. Different rows represent, from top to bottom, the different trigger types: pointlike events in the ID (TT1 & BTB0); internal muons (TT1 & BTB4); external muons (TT2); and the three types of calibration triggers: laser (TT8), pulser (TT32), and random (TT64). The total number of decoded hits detected for each trigger type is shown on the top right corner of the respective pads. ID hit time [ns] OD hit time [ns] FIG. 4. Analogy of the Fig. 3, but the event structure is shown for individual events of each trigger type. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) Neutron-muon pulse shape ID hit time [ns] Single muon event followed by neutrons ID hit time [ns] FIG. 5. Top: Event structure of TT1 & BTB4 muons (blue) and the start of the follow-up 1.6 ms-long TT128 neutron gates (red) collected in a typical 6 hour run. Due to the large amount of light for muon events, the detected hits at the end of TT1 & BTB4 events are prescaled, causing the step between the two trigger types. Bottom: example of a TT1 & BTB4 muon (blue) followed by a TT128 event (red) with high neutron multiplicity. 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 antineu­trino candidates per year. Below is a summary of different categories of muon detection in Borexino. 1. Muon trigger flag (MTF) MTF muons trigger the OD and set the bit BTB4 of the trigger board, as described above in Sec. III A. 2. 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 FADC DAQ gate [ns] FIG. 6. Examples of the FADC 1.28 µslong waveforms with 2.5 ns binning: 3.9 MeV pointlike event (top), muon with 90 MeV deposited energy (middle), calibration pulser event (bottom). 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. 3. Inner detector flag (IDF) The IDF identifies muons based on different time profiles of hits originating from muon tracks with respect to those from pointlike 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 optimized in three different energy ranges. For muons, the Nh energy esti­mator is proportional to the track length across the ID rather than to the muon energy. Events with Nh >2100 and with mean time greater than 100 ns are considered as muons. In the lower energy interval 80 30.40. ns are tagged as muons above (below) FIG. 7. Mutual efficiencies of the three strict muon flags MTF (top), MCF (middle), and IDF (bottom) for events with Nh > 80 calculated for calendar years (2019 contains data up to end of April). Nh 1900, respectively. In addition, the Gatti .=ß dis­crimination parameter (see Sec. VII D) is used in IDF to reject electronics noise, retriggering of muons, as well as scintillation pulses coming from the buffer. In the very low energy interval Nh < 80 dominated by 14C background, IDF is not applied due to the limited performance of pulse­shape identification. A detailed explanation of the IDF flag can be found in [77]. 4. 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 TT2 events, very few also TT1 triggers, that satisfy MTF or MCF (slightly modified) conditions. 5. 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 TT1 events. However, we include in this category also the rare TT2 events, which did not trigger the ID, but have a cluster with more than 80 hits in the ID-data. In the lack of a pure muon sample, we introduce the concept of mutual efficiencies. We first define the muon reference sample with one flag, with respect to which we express the efficiency of the other two flags. Such mutual efficiencies of the three strict muon flags for muons with Nh > 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. IX C). 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 × .1–0.9989.1.1.31 0.5.× 10-5 . 6. 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 TABLE III. Mutual efficiencies for the three muon flags MTF, MCF, and IDF are given for the period December 2007—April 2019 used in the geoneutrino analysis. In the three last columns, the muon reference sample was defined, from left to right, by IDF, MTF, and MCF flags, respectively. In the last 4 rows, the IDF efficiency is shown for different energy ranges as well. Visible energy (Nh) Mutual efficiency . vs IDF vs MTF vs MCF .80 .80 .80 .MTF .MCF .IDF 0.9957 0.0004 0.9894 0.0004 Not applicable Not applicable 0.9927 0.0004 0.9882 0.0004 0.9989 0.0004 Not applicable 0.9881 0.0004 80–900 900–2100 .2100 .IDF .IDF .IDF Not applicable Not applicable Not applicable 0.7941 0.9988 1.0000 0.0014 0.0008 0.0005 0.7921 0.9984 1.0000 0.0014 0.0008 0.0005 M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) decreased purity of the muon sample. In fact, special tags 500 mark as muons also noise events. In the antineutrino analysis, all special muons are conservatively treated as internal muons (see Sec. VI A), regardless of the fact whether they have a cluster of hits in the ID-data. 400 Decoded hits/ 100 ns (i) 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 300 200 NR > 200 when only less than 5% of them are decoded (ND=NR < 0.05). These events can have 100 trigger type TT1 or TT2, without any condition 0 -15000 -10000 -5000 regarding the BTB trigger word. ID hit time [ns] (ii) Special flag 2 This flag tags events of trigger type Special flag 6 TT1 or TT2 with ND > 100 that have the 22 bit of the trigger word raised, but the OD triggered in coincidence with the service interrupt which gen-1200 erated service triggers. In this case, additional bits of 1000 800 600 400 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. (iii) Special flag 3 Events with a cluster start time out of the DAQ gate, as shown in Fig. 8, are tagged as special muons. These events typically tag muons Decoded hits/ 100 ns 200 that pass through the detector during the 2–3 µs dead time after TT1 & BTB0 events. The muon itself does not generate a trigger, but the hits from the PMT after-pulses do. Even if the detection efficiency for out-of-gate hits is < 100%, the main muon pulse is detected and positioned before the start of the DAQ gate. (iv) Special flag 4 The pointlike events are expected to have cluster mean time smaller than 200 ns. Conservatively, all events with a cluster mean time greater than 200 ns are tagged as special muons. (v) Special flag 5 There is an extremely small number of TT1 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. (vi) Special flag 6 The events of trigger type TT8, TT32, or TT64 (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. 7. 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 0 -15000 -10000 -5000 0 ID hit time [ns] FIG. 8. Top: Muon type special flag 3 showing the muon which crossed the detector during the dead-time after the previous trigger. While the main peak is out of gate, the muon was detected anyway, triggering on the after-pulse. The main muon peak is out of the gate. Bottom: Muon special flag 6 type showing the pulser trigger TT32 with a muon occurring during the same gate. The vertical dashed lines represent the start of the DAQ gate at -16 µs. 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 multilayer 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 6 of reliability is suitable for the antineutrino analysis. In this case, the maximum tagging efficiency is achieved with the 4 price of the greatest over-efficiency. The combined muon tagging efficiency of the main and the FADC DAQ systems is 99.9969% [81]. 2 5 4.5 4 3.5 3 0 z [m] 8. Internal large muon flag The term internal large muon flag is used for the muon category that considers as internal muons not only the strict internal flag muons, but includes also special and FADC muon flags. This approach is conservative, since internal muons, contrary to external muons, are able to create cosmogenic background other than neutrons and thus, generally, require longer veto. This will be discussed in Sec. VII A. C. Inner vessel shape reconstruction The shape of the thin nylon 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 estimated 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 behavior. 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. VII F) 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 (Npe 1 290–350 p:e:) recon­structed on the IV surface (Fig. 9). These events originate in the radioactive contamination of the nylon and are domi­nated by 210Bi, 40K, and 208Tl. The reconstructed position of selected events is fit assuming uniform azimuthal symmetry (x-y 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 par­ticular the reconstructed radius R as a function of .. Fixed parameters of the fit are the two endpoints (. 1 0, . 1 .) 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 2 -2 1.5 1 -4 0.5 -6 0 -6 -4 -20 2 4 6 x [m] FIG. 9. A cross section (z-x plane, jyj < 0.5 m) view of the distribution of 2478 events, acquired during a 3 week period, and selected for the IV shape reconstruction. The color axis represents the number of events per 0.0016 m3, in a pixel of 0.04 m× 1.00 m× 0.04 m(x × y × z). This distribution reveals the IV shift and deformation with respect to its nominal spherical position shown in solid black line. 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 IV. 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. XI C. D. .=ß 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 discrimina­tion parameters which aid in the distinction of . 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 Radon-correlated 214Bi.ß-.–214Po... coincidence sample that was present in the detector during M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 4.6 4.4 4.2 4 3.8 . [rad] (a) 1. Gatti optimal filter The Gatti optimal filter (G) is a linear discrimination technique, which allows us to separate two classes of events with different time distributions [20,82]. First, using the typical ß and . time profiles after the time-of-flight subtraction, the so-called weights w.t. are defined for n time bins t: n P..t.- Pß.t. nn w.t.. ; .7. n P..tn..Pß.t. n where P..t. and Pß.t. are the probabilities that a nn photoelectron is detected at the time bin tfor . and ß n events, respectively. The Gatti parameter G for an event IV volume [m3] Radius [m] 315 310 305 300 295 290 2008 2010 2012 2014 2016 2018 Year (b) FIG. 10. (a) Examples of the reconstructed IV shapes, i.e., reconstructed IV radius Ras a function of ., resulting from fitting of the (R, .) distributions of the selected events originating in the IV contamination. The solid blue line shows one recent vessel shape. The dashed lines represent the least (green) and the most (red) deformed IV shapes registered. The dashed black line shows the ideal sphere with 4.25 m radius. (b) Trend of the reconstructed IV volumes as a function of time (1 week bin), starting from 2007, December 09 up to 2019, April 28. Each point represents 3 weeks of data. The error bars reflect the goodness of the fit of (R, .) distributions. The three dashed vertical lines represent the LS refill into the IV: we observe that the increase in the reconstructed volume corresponds to the amount of inserted LS. During the two periods (small red arrows in the lower part of the plot), the concentration of the DMP in the buffer was decreased from the original 5.0 g=l to 3.0 and then to 2.0 g=l. Thanks to the resulting better match between the densities of the LS and the buffer liquid, the rate of the leak was minimized, but not fully stopped. the 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. VII D. with the hit time profile f.t., after the time-of-flight n subtraction, is then defined as: X G1 f.tn.w.tn.: .8. n Figure 11(a) shows the distributions of the Gatti parameter for ß particles from 214Bi (G<0) and . particles from 214Po (G>0). 2. Multilayer perceptron The multilayer perceptron (MLP) is a nonlinear 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 pulse­shape variables, parametrizing the event hit-time profile, as input. Among these variables are, for example, tail-to­total ratio for different time bins t, mean time of the hits in n 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). IV. ANTINEUTRINO DETECTION Antineutrinos are detected in liquid scintillator detectors through the inverse beta decay (IBD) reaction illustrated in Fig. 12: .—e.p› e..n; .9. in which the free protons in hydrogen nuclei, that are copiously present in hydrocarbon (CnH2n) molecules of organic liquid scintillators, act as target. IBD is a charge­current interaction which proceeds only for electron fla­vored 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 Gatti parameter (a) MLP parameter (b) 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 Ep, which is directly correlated with the incident antineu­trino energy E—: .e Ep ~ E.—e - 0.784 MeV: .10. The offset results mostly from the difference between the 1.806 MeV, absorbed from E—in order to make the IBD .e 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 scintilla­tion light that is detected in a single coincident delayed FIG. 12. Schematic of the proton inverse beta decay interaction, used to detect geoneutrinos, showing the origin of the prompt (violet area) and the delayed (blue area) signals. The visible energy of the prompt signal includes the contribution from the kinetic energy of the positron as well as from its annihilation. The neutron thermalizes and scatters until it is captured on a free proton. The 2.2 MeV deexcitation gamma of the deuteron represents the delayed signal. signal. In Borexino, the neutron capture time was measured with the 241Am–9Be 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 Npe spectrum of delayed signals due to the gammas from captures of neutrons emitted from Q [p.e.] d FIG. 13. The Npe charge spectrum of delayed signals, ex­pressed in the number of detected photoelectrons, due to the gammas from captures of neutrons emitted from the 241Am–9Be calibration source placed in the center of the detector. The clearly visible peaks of 2.22 MeV and 4.95 MeV gammas from the neutron captures on proton and 12C are positioned at 1090 p.e. and 2400 p.e., respectively. The other peaks are from neutron captures on nuclei of stainless steel used in the source and its insertion system construction: at >7 MeV energies due to captures on (Fe, Ni, Cr) and at 477.6 keV on 10B. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) the 241Am–9Be 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 .soriginate from the neutron captures on 12C 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 signature of .—e interactions, which strongly suppresses backgrounds. In the following sections, we refer to these signals as prompt and delayed, respectively. V. EXPECTED ANTINEUTRINO SIGNAL This section describes the expected antineutrino signals at the LNGS location (r .142.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. IV) over 1 year by a detector with 100% detection efficiency containing 1032 free target protons (roughly corresponds to 1 kton of LS). Since we detect only electron flavor of the total anti­neutrino flux (Sec. IV), neutrino oscillations affect the expected signal expressed in TNU. Thus, the neutrino oscillations and the adopted parameters are discussed in Sec. VA. The evaluation of the expected geoneutrino signal from the Earth’s crust and mantle is described in Sec. VB. Section VC details the estimation of the signal from antineutrinos from the world reactors, the most important background for geoneutrino measurements. Atmospheric neutrinos, discussed in Sec. VD, 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. VE. The final number of the expected events from each of these sources, expected in our dataset (Sec. III A) and passing all optimized selection cuts (Sec. VII H) will be presented in Sec. IX B. A. Neutrino oscillations The presently accepted Standard Model of elementary particles describes neutrinos existing in three flavors (electron, muon, and tau) with masses smaller than 1=2 of the Z0 boson mass. The experiments with solar, atmospheric, as well as reactor antineutrinos observed that the neutrino flavor can change during the travel between the source and the detector. The process of neutrino oscilla­tions has been established and confirmed that neutrinos have a nonzero 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 antineu­trino flux at a given site, it is necessary to consider the survival probability Pee of the electron flavored 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 Pee for MeV antineu­trinos, 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: Pee.L; E.—e .11 - 4c4 12c2 sin2. 13s2 12 - 4s2 13c2 sin2.. ..=2. 13c2 12 - s2 c2 s2 sin2.-. ..=2.; .11. 131312 where .m2L .m2L . 1 . 1 4E—4E— .e .e 1 .m2 1.m2 .m2 1j.m2 ..m2 j 12 3132 2 .mi;j 1m2 i - m2 j c12 1cos .12 s12 1sin .12 c13 1cos .13 s13 1sin .13 .12. and L and E—are expressed in natural units (. 1c 11). .e We assume normal hierarchy for neutrino mass eigen­states .m1 0.8 was applied to the delayed using the better .=ß discrimination power of the MLP (Sec. III D) when compared to the Gatti. The cut threshold was chosen based on the 241Am - 9Be calibration data as shown in Fig. 29. It was found that only a 5.4.8.1.× 10-3 fraction of the neutron-capture .s remains below the MLP threshold of 0.8 for the source position in the detector’s center (close to IV border).” MLP FIG. 29. Distribution of the MLP .=ß discriminator (Sec. III D) for the gammas from the capture of neutrons from 241Am - 9Be calibration source placed at the detector center .x; y; z.1 .0; 0; 0.m, as well as at .0; 0; -4.m, where the observed distribution is the broadest. The dashed line shows the threshold set at 0.8 for the delayed in IBD search. E. Energy cuts Energy cuts were applied to the prompt and delayed to aid the identification of IBD signals. The analysis was performed using the charge energy estimator Npe (Sec. III A). The intervals in the respective charges Qp and Qd are optimized as explained in this section. 1. Charge of prompt signal The energy spectrum of the prompt Ep 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 Qmin 1408 p:e: This threshold value p corresponds to approximately 0.8 MeV and remains unchanged from previous analyses. No upper limit is set on the charge of the prompt candidate. 2. Charge of delayed The delayed signal can be either due to a 2.22 MeV (n-capture on 1H), or a 4.95 MeV gamma (n-capture on 12C) with about 1.1% probability, as described in Sec. IV. The corresponding values in the Npe 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). Q [p.e.] d FIG. 30. The Npe spectrum of delayed candidates from the 241Am - 9Be calibration source placed at .x; y; z.1.0; 0; -4.m inside the detector (dashed line) compared to the spectrum for the source placed at the detector’s center (solid line). The two peaks from the left correspond to 2.2 MeVand 4.95 MeV gammas from neutron captures on 1H and 12C nuclei, respectively, while the higher energy peaks are from neutron captures on stainless steel nuclei (Fe, Ni, Cr) used in the source construction. All gamma peaks in the off-center spectrum are shifted to lower energies and develop tails due to the partial energy deposits in the buffer. The spectra are normalized to one. In the last geoneutrino analysis [18], Qmind was set to 860 p.e., a conservatively large value because of the radon­214Po.. ... correlated background, particularly due to decays. This was discussed in Sec. VI E. In this analysis, Qmin was decreased to 700 p.e. based on the improved d performance of .=ß separation with MLP (Sec. VII D), 214Po.. ... which improved the ability to suppress 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 Qmin 860 p.e. was d retained. The Qmin cut was not decreased below 700 p.e. for d the following reasons: (1) The Npe spectrum of accidental background in­creases at lower energies, as shown in Fig. 37. (2) The endpoint of the . peak from the main 214Po decay in the radon-correlated background is ~600 p:e:,as described in Sec. IX G. In this analysis, we include the neutron captures on 12C (4.95 MeV) as well. Consequently, Qmax was increased d from 1300 p.e. (2.6 MeV) to 3000 p.e. (.5.5 MeV). F. 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. III C). A DFV cut, i.e., a require­ment of some minimal distance of the prompt from the IV, dIV, is applied along the reconstructed IV shape. In the previous analysis [18] a conservative cut of dIV 130 cm has been applied to account for the uncertainty of the IV shape reconstruction and the potential background coinci­dences near the IV. In this analysis, the DFV was increased by using dIV 110 cm, which leads to a 15.8% relative increase in exposure. This choice is justified below. The geoneutrino sensitivity studies (Sec. X)were performed for different combinations of the DFV cut and the Qd d does min, as shown in Fig. 31. The choice of Qmin not have a big impact on the expected precision, for a given DFV cut. Order of 5% improvement in the statistical uncertainty of the geoneutrino measurement is expected when the cut is lowered to dIV 110 cm, while there is no further improvement when no DFV cut is applied. A dIV 110 cm DFV cut is also sufficient to account for the precision of the IV shape reconstruction (Sec. III C). In addition, we have verified that no excess of the IBD candidates is observed at large radii (close to the IV), as it will be shown in Sec. XI A. G. Multiplicity cut The multiplicity cut requires that no additional “high­energy” (Npe > 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 M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 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 Npe > 400 p:e: events within 2 ms is of the order of 0.01%, which is negligible. H. Summary of the selection cuts The summary of all the optimized selection cuts is listed in Table XII. TABLE XII. Summary of the optimized selection cuts for IBD candidates: charge cut on the prompt, Qp and delayed, Qd, time and space correlation dt and dR, respectively, distance to the IV, dIV, MLP .=ß particle identification parameter cut on delayed, and multiplicity cut. The scheme indicating the duration and geometry of the different muon vetoes is shown in Fig. 23. Cut Condition Qp >408 p:e: Qd (700–3000) p.e. (860–3000) p.e. (WE period) dt Double cluster: .2.5–12.5. µs Single cluster: .20–1280. µs dR 1.3 m Muon veto 2 s or 1.6 s or 2 ms (internal µ) 2 ms (external µ) dIV 10 cm (prompt) PID (.=ß) MLPd > 0.8 Multiplicity No Npe > 400 p:e: event 2 ms around prompt/delayed VIII. MONTE CARLO OF SIGNAL AND BACKGROUNDS The spectral fit of the prompt charge Qp (Sec. XA), relies extensively on the use of MC-constructed prob­ability distribution functions (PDFs) which represent the shapes of geoneutrino signal and, with the exception of accidental coincidences (Sec. IX D), all backgrounds. The construction of these PDFs is described in Sec. VIII A.The GEANT4 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 241Am - 9Be neutron source is of particular impor­tance, since the delayed IBD (Sec. IV) signal is repre­sented by a neutron. The comparison of the neutron spectra from the 241Am - 9Be calibration source at the detector center and at .0; 0; -4. m position inside the detector is shown in Fig. 30. A. Monte Carlo spectral shapes Once the full G4Bx2 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. XA). 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 dataset from December 2007 to April 2019 is simulated individually. After the simulation, the optimized geoneutrino selection cuts (Sec. VII H) 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. V are used. The antineutrino energy spectra are transformed to positron energy spectra follow­ing Eq. (10). In particular, for geoneutrinos, the energy spectra as in Fig. 18 are used. Individual spectra from 232Th and 238U chains were also simulated, so that they can be weighted according to the expected Rs ratio (Eq. (18) for different geological contributions. For reactor antineutri­nos, 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 nonantineutrino backgrounds are shown in Fig. 33. A dedicated code is developed within the G4Bx2 simulation framework for the generation of 9Li events, based on Nuclear Data Tables and literature data [125]. The input for (.;n) background simulation is discussed in Sec. VI C, while for atmospheric neutrinos in Sec. VD. 0.05 0.03 0.025 0.04 0.02 0.03 0.01 0 0 500 1000 1500 2000 2500 3000 3500 4000 Qp [p.e.] Qp [p.e.] 500 1000 1500 2000 2500 3000 3500 4000 (a)(b) 1.2 Events/ 10 p.e. Events/ 10 p.e. Events/ 10 p.e. 0.015 0.02 0.01 0.006 0.005 0.004 0.003 0.002 0.001 1.15 1.1 Ratio/ 100 p.e. 1.05 1 0.95 0 500 1000 1500 2000 2500 3000 3500 4000 Q [p.e.] Q [p.e.] p p (c)(d) FIG. 32. The MC-based PDFs normalized to one, i.e., the expected shapes of prompts for geoneutrinos and reactor antineutrinos, after optimized geoneutrino selection cuts, that include the detector response. Top-left: geoneutrinos with Th/U ratio fixed to the chondritic value (RS 1 0.27). Top-right: 238U and 232Th PDFs shown separately. Bottom-left: reactor antineutrinos “without 5 MeV excess.” Bottom-right: ratio of reactor antineutrino spectra “with/without 5 MeV excess” (Sec. VC), normalized to the same number of events each, in order to demonstrate the difference in shape only. The PDFs of prompts due to antineutrinos from a hypothetical georeactor (Sec. VE) 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. B. 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. VII H) 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. XI C. The error on the detection efficiency was estimated based on the comparison of the calibration data (most importantly 241Am - 9Be neutron source data) with MC simulation. The major contribution comes from the uncertainties of the detector response close to the edge of IV. IX. 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 opti­mized selection cuts (Table XII) is described. The dataset and the total exposure are presented in Sec. IX A. The number of expected antineutrino events from different sources, based on the estimated antineutrino signals as in Table VIII, is presented in Sec. IX B. The next sections treat the nonantineutrino background, following the struc­ture of Sec. VI, where the physics of each of these M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) Q [p.e.] (a)p Events/ 10 p.e. Events/ 100 p.e. Qp [p.e.] (b) after optimized IBD selection cuts, normalized to one. We show both the case of nonoscillated spectrum, as well as the oscillated spectra for the georeactor placed in three positions at different depths: GR1 (d 12900 km), GR2 (d 16371 km), and GR3 (d 19842 km), defined in Fig. 21(a). A. Data set and exposure In this analysis the data taken between December 9, 2007 and April 28, 2019, corresponding to tDAQ 13262.74 days of data taking, are considered. The average life-time weighted IV volume and the FV used in this analysis — (DFV cut dIV 110 cm) are VIV 1.301.3 10.9.m3 and — VFV 1.280.1 10.1.m3 , respectively, after taking the Qp [p.e.] (c) changing shape of the IV into account (Sec. III C). These correspond to the average IV and FV mass of m—IV 1 .264.59.6.ton and m—FV 1.245.88.7.ton, respec­tively, considering the scintillator density of .LS 1.0.878 0.004.gcm-3 . The total exposure after the cosmogenic veto (Sec. VII A) and in the FV is E 1.2145.8 82.1.ton × yr, considering the systematic uncertainty on position reconstruction (Sec. XI C). This can be expressed as Ep 1.1.29 0.05.× 1032 protons × yr, using the proton density in Borexino FIG. 33. MC-based PDFs of prompts for different backgrounds after optimized geoneutrino selection cuts, normalized to one. Top: cosmogenic 9Li background. Middle: (., n) background. Bottom: atmospheric neutrino background. background categories is described. In particular, cosmo­genic background is discussed in Sec. IX C, accidental background in Sec. IX D,(alpha, n) interactions in Sec. IXE,(., n) and fission in PMTs in Sec. IX F, radon correlated background in Sec. IX G, and finally 212Bi - 212Po background in Sec. IX H. Section IX I summarizes the expected total number of nonantineutrino background events. TABLE XIII. Detection efficiencies after the optimized selec­tion cuts for geoneutrinos from 238U and 232Th chains individually and for their summed contribution (according to the chondritic Th/U mass ratio), as well as for antineutrinos from reactors and from a hypothetical georeactor. The efficiencies were determined based on the MC simulation of each component. The error is estimated using the calibration data. Source Efficiency [%] 238U geoneutrinos 232Th geoneutrinos Geoneutrinos (Rs 10.27) Reactor antineutrinos 87.6 84.8 87.0 89.5 1.5 1.5 1.5 1.5 Georeactor 89.6 1.5 LS of Np 1.6.007 0.001.×1028 protonston-1. Applying the geoneutrino detection efficiency described in Sec. VIII B, the effective exposure for the geoneutrino detection reduces to E01.1866.4 78.4.ton×yr and E0 p 1.1.12 0.05.× 1032 protons × yr. B. Antineutrino events This section summarizes the number of expected anti­neutrino events detected with the optimized selection cuts, assuming the expected antineutrino signals as given in Table VIII. An overview is given in Table XIV. C. Cosmogenic background The various cosmogenic backgrounds in Borexino which affect the geoneutrino analysis were explained in Sec. VI A, while the analysis and technical evaluation of these back­grounds is explained in this section. TABLE XIV. Summary of the expected number of antineutrino events with optimized selection cuts in the geoneutrino (408– 1500) p.e and reactor antineutrino (408–4000) p.e. energy ranges. The errors include uncertainties on the predicted signal only. The reference exposure Ep 1.1.29 0.05.× 1032 protons × yr cor­responds to the analyzed period. Energy range Signal Model [p.e.] [Events] Geoneutrinos Bulk lithosphere 408–1500 28.8.5.5 -4.6 CC BSE (total) 408–1500 31.9.6.2 -5.4 CC BSE (mantle) 408–1500 2.80.6 GC BSE (total) 408–1500 38.8.6.2 -5.4 GC BSE (mantle) 408–1500 9.80.9 GD BSE (total) 408–1500 51.1.6.3 -5.5 GD BSE (mantle) 408–1500 22.01.2 FR (total) 408–1500 62.0.6.4 -5.6 FR (mantle) 408–1500 33.01.7 Reactor antineutrinos without 408–1500 42.60.7 “5 MeV excess” 408–4000 97.6.1.7 -1.6 with 408–1500 39.50.7 “5 MeV excess” 408–4000 91.9.1.6 -1.5 Atmospheric neutrinos 408–1500 2.21.1 408–4000 3.31.6 408–8000 9.24.6 1 TW Georeactor GR2: Earth’s center 408–1500 3.60.1 408–4000 8.90.3 GR1: CMB at 2900 km 408–1500 17.60.5 408–4000 43.11.3 GR3: CMB at 9842 km 408–1500 1.50.04 408–4000 3.70.1 1. Hadronic background The hadronic background, expected to be dominated by 9Li (Sec. VII A), 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 9Li 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. VII A and summa­rized in Fig. 23. a.Time distribution dtLi-µ.—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 QLi p charge energy spectrum of the prompts is compatible with the expected MC spectrum of 9Li, as it is shown in Fig. 35(a). The distribution of the time differences between the prompt and the preceding muon, dtLi-µ, is shown in Fig. 35(b). As it can be seen, no events are observed after dtLi-µ > 1.6 s. The decay time . is extracted by performing an exponential fit to the dtLi-µ distribution and is found to be .0.260 0.021.s. This is compatible with .9Li 10.257 s decay time of 9Li. b.Spatial distribution dRLi-µ.—The distance of the 9Li prompt from the muon track, dRLi-µ, 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 normalization factor n: n 2µ ..2.-1 - 2dRLi-µ f.dRLi-µ; .; .; µ;n.1 × exp 2. 2. µ ..2.-1 - dRLi-µ × erfc p ... : 2. .30. The fit results in . 1.0.35 0.11.m, which represents well the combined position reconstruction of the muon track and the prompt, and in . 1.0.68 0.24.m. Considering these values and the fit function in Eq. (30), 2.75% of 9Li 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 catego­ries (Fig. 23) is evaluated, considering the dtLi-µ and dRLi-µ parametrizations described above. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) Charge distribution of 9Li (i) 1.6 s and 2.0 s vetoes of the whole detector: For 8.1(7.1)% of the total muons, we veto the whole detector for 1.6(2.0) s. In the time window [2 ms, 1.6(2.0) s] after these muons, whereIBD-likecandidatesweresearched, falls 20 15 10 5 0 0 1000 2000 3000 4000 5000 6000 Q [p.e.] p (a) Time distribution of 9Li 99.02(99.19)% of the corresponding 9Li back­ground. In this window, 282(23) candidates were found. Finally, in the time window after the time veto, the remaining background is exp.-1.6.2.0.=..10.21.0.05.% of the total back­ground. Adding the two components and including the statistical error and the error on ., the total amount of hadronic background in the antineutrino 0.09 candidate sample is 0.18.events. -0.06 (ii) 1.6 s cylindrical veto: For 27.0% of the muons with a reliably recon­structed 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 ob­served in the whole detector. Thus, after 1.6 s, we expect 0.015.0.013 events distributed in the -0.011 whole detector. Of these, 2.75% would be recon­structed 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 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 (2.20.2) expected IBD events. We observe one dtLi-µ [s] additional candidate. (b) (iii) 2 ms veto of the whole detector: For 57.8% of all muons which have a lower Spatial distribution of 9Li probability to produce detectable hadronic back- ground [.µ - n.], we restrict the time veto of <8000 the whole detector to 2 ms, to veto only the cosmogenic neutrons. Consequently, the exposure increases by 6.6%, corresponding to (9.70.8) expected IBD candidates. Seven candidates were observed, well within the expectation. However, 15 10 this does not guarantee that we did not introduce additional 9Li 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 analy­sis. It is reasonable to assume that the production ratios for 9Li 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 9Li decays are the same. Thus, 0 0123456 dRLi-µ [m] (c) FIG. 35. Distributions of the observed 9Li cosmogenic back­ground, i.e., the IBD-like candidates after optimized selection cuts following muons after at least 2 ms. (a) Comparison of the prompt QLi p charge data spectra (points) with the MC spectrum (solid line). (b) Distribution of dtLi-µ, the time difference between the prompt and the preceding muon, fit with the exponential function. (c) Dis­tribution of dRLi-µ, the distance of the prompt from the well reconstructed muon track, fit with the function in Eq. (30). the equality of the ratios of the number of observed 9Li candidates (with decay neutron, NLi) and cosmogenic neutrons (Nn) should hold: .NLi 1 305.>8000 .NLi.<8000 1 : .31. .Nn 1 8.6×105.>8000 .Nn 1 9181.<8000 From this equality, the expected number is 3.2 events. This is the number of .NLi.<8000 expected 9Li events produced by muons with less than 8000 hits, independent of whether this muon wasfollowedbyaneutronornot.Thisisa conservative number for the background estima­tion 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 9Li than .µ - n.<8000 muons. However, we do not observe any 9Li candidate for .µ . n.<8000 muons. Thus to summarize, our expected 9Li background for .µ - n.<8000 muons is 3.21.0 events, which includes larger systematic error. In total, after summing all the contributions, the expected 9Li background within our golden IBD candidates is (3.61.0) events. 2. Untagged muons The .0.0013 0.0005.% mutual inefficiency of the strict muon flags shown in Sec. III B corresponds to (195 75) undetected muons in the entire dataset. Following the discussion in Sec. VI A, these muons could eventually cause background of three types: (i) µ . µ: Considering the small amount of undetected muons in the entire dataset, the probability that two undetected muons would fall in a 1260 µs time window of the delayed coincidence is completely negligible. (ii) µ . n: The most dangerous are pairs of buffer muons (possibly fulfilling the Qp cut) followed by a single neutron (multiple neutrons are removed by the multiplicity 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.70.003. × 10-5 . Hence, there will be (0.019 0.007) events of this kind in the IBD sample due to the untagged muons. (iii) Muon daughters: After 1.497 × 107 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 back­ground originating from untagged muons is (0.023 0.007) events in the IBD sample. We note that this is a very small number. 3. Fast neutrons As described in Sec. VI A, 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 [126]. 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. III A, 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 simu­lations are used to obtain an estimation of fast neutron background due to the muons passing through the sur­rounding rock and not the detector. Both estimations are given below. a.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: (i) 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. (ii) 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 pointlike 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. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) b.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 [126] 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 pen­etrating 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. D. Accidental coincidences In order to evaluate the amount of accidental coinci­dences in the antineutrino sample, coincidence events were searched for in the off-time interval dt 112 s; 20 s and were then scaled to the 1270 µs duration of the geoneutrino selection time window (dt 112.5 µs; 12.5 µs .120 µs; 1280 µs , Sec. VII H). 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 112 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 112 ms; 2 s is as follows: racc01racc0 · exp.-rµ · dt.; .32. 0 where racc0 would be the rate of accidental background with 0 relaxed cuts and without the muon veto suppression factor. After 2 s, the exponential suppression factor becomes constant and consequently the fit function for dt > 2 s acquires a constant form: dt(delayed-prompt) [s] (a) dt(delayed-prompt) [s] (b) FIG. 36. Distribution of dt (delayed-prompt) for accidental coincidences (a) with relaxed selection cuts to show a decreasing trend until 2 s and a constant trend after 2s and (b) with geoneutrino selection cuts in the time window [2 s, 20 s]. In both cases the search is performed by applying a 2 s veto for all internal muons. racc01racc0 0 · exp.-rµ · 2s.: .33. When fitting the spectrum with relaxed cuts, rµ 1.0.0501 0.0016.s-1 was obtained. This is compatible with the measured rate of internal muons of .0.05311 0.000 01.s-1. The validity of this behavior 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.999 93 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 racc for the geoneutrino measurement, the dt 112 s; 20 s 0 distribution of 49 004 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: Q [p.e.] d 400 p.e. and scaled to the number of prompts in the solid-line spectrum. The dashed vertical line shows the chosen Qmin d charge threshold of 700 p.e. 1racc racc 0 · exp.-rµ · 2s.: .34. 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 racc is .3029.0 12.7.s-1. This means 0 that the number of accidental coincidences among our IBD racc candidates can be estimated as × 1270 µs, that is 0 (3.846 0.017) events. The Npe spectra of the prompt and delayed signals of the accidental coincidences, selected with optimized geo­neutrino cuts in dt 112 s; 20 s time window, are shown in Fig. 37. E. (., 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 210Po that is found out of equilibrium with the rest of 238U chain [20]. In the energy region of 210Po (Npe 1 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.210Po.1.12.75 0.08.events=.day · ton. is used to evaluate the (., n) background from the 210Po contami­nation of the LS. In the second stage, the neutron yield, i.e., the probability that 210Po . would trigger an (., n) reaction in the LS, was calculated with the NeuCBOT program [127–130], which is based on the TALYS software for simulation of nuclear reactions [122,131,132]. Only PC was considered as a target material. The contribution from PPO is negligible, as 70 60 50 40 30 20 10 0 200 400 Time [weeks] FIG. 38. The evolution of the weekly 210Po .-decay rates in the DFVof geoneutrino analysis, in the period from December 2007 to April 2019. The horizontal line shows the mean rate — RDFV.210Po.1.12.75 0.08.events=.day · ton.. 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 experi­ .-decay rate [events/ (day·ton)] mental data for energies of the 210Po . 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 Yn of the (., n) reaction in the PC is found to be .1.45 0.22.× 10-7 neutrons per a single 210Po decay. Note that the corresponding value used in the previous Borexino geoneutrino analysis, based on [133], was approximately three times smaller. Even if the (.;n) background is directly proportional to Yn [see Eq. (35) below], this has a negligible impact on the final geoneutrino result, thanks to a high radio-purity of the Borexino scintillator. The final calculation of the number of IBD-like coinci­dences N..;n. triggered by the neutrons from the (., n) reaction over the whole analysis period can be computed using the following formula: 210Po.· E · Yn · ...;n.; N..;n.1 R—DFV..35. where E 1.2145.8 82.1.ton × yr (Sec. IX A) is the exposure and ...;n.156% 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 210Po contamination of the LS is (0.81 0.13) events. Another potential source of background are (.,n) reac­tions due to 210Po 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 occur­ring in the inner buffer, this probability was estimated to M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 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 210Po 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 [134], that not any more opera­tional. 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 210Po contamination in the buffer. In January 2010, the DMP concentration in the buffer was reduced to 2 g=l (the original concentration was 5 g=l), as discussed in Sec. III. Since then, no further operations have been performed with the buffer and the 210Po contamination is expected only to decay (. 1 199.6 day) and to be suppressed 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 210Po 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. XI A, the golden IBD candidates are evenly distributed in time and no excess close to the IV is observed. F. (., n) interactions and fission in PMTs In order to obtain an upper limit to the possible back­ground 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 Npe. 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 2His .D 1 1.6 mb. Since the deuteron density is .D 1 7.8 × 1018 atoms=cm3, the upper Events/ 5 p.e. 2000 1500 1000 500 0 Npe FIG. 39. The Npe distribution of 3 MeV .s generated in the DFV and in the entire analyzed period of geoneutrino analysis. The vertical line indicates the 1200 p.e. threshold used in the evaluation of the (., n) background. limit on the number of events N..;n. due to this background, taking into account the estimated detection efficiency ...;n.1 50%,is N..;n. 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 Qmin to 700 p.e., safely above the 214Po .-peak d [Fig. 41(b)]. The total number of background events correlated with radon contamination is expected to be (0.003 0.001), which is completely negligible. H. 212Bi - 212Po background A MC study proves that the cut on Qmin 1 860 p:e:, d adopted to reject the radon contamination during the WE 214Po charge spectra for 80 ton fiducial volume Q [p.e.] d (a) 214Po decays (MC) Q [p.e.] d (b) FIG. 41. Top: Comparison of the charge spectra of . decays from 214Po from data (circles with error bars) and MC (solid line) in the fiducial volume of ~80 tons around the detector center chosen to tune the alpha particle quenching factor. Bottom: Charge distributions for the three 214Po decays (Table XI) obtained using MC with the . energy scale tuned on pure .­decays [Fig. 41(a)]. The main . decay (red line) and the two subdominant (. . .) branches are shown: 10-4 (green line) and 10-7 (blue line) branch. The three spectra are normalized to have the same area. periods, is effective in removing the 212Bi - 212Po fast coincidences, as shown in Fig. 42. It can be seen that 212Po . the endpoint of the peak is around 700 p.e. Therefore, a Qmin of 700 p.e., combined with the MLP d pulse shape cut on the delayed, makes the overall 212Bi - 212Po background fully negligible in geoneutrino analysis. I. Summary of the estimated nonantineutrino background events Table XV summarizes the expected number of events from all nonantineutrino backgrounds passing the opti­mized selection cuts listed in Table XII. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) Events/ 5 p.e. 350 300 250 200 150 100 50 0 350 400 450 500 550 600 650 700 750 800 Qd [p.e.] FIG. 42. Charge spectrum of MC generated 212Po . peak. The vertical line at 700 p.e. shows the optimized Qmin for IBD d selection. X. SENSITIVITY TO GEONEUTRINOS This section describes the Borexino sensitivity to geo­neutrinos and the MC based procedure with which it was evaluated. In Sec. XA the description of the basic ingre­dients of the analysis focuses on the spectral fit of the Qp spectrum. Section XB describes the sensitivity tool that performs such fits on 10 000 QMC spectra, each corre­ p sponding to a MC generated pseudoexperiment. The expected precision for the Borexino geoneutrino measure­ment as well as its sensitivity to the mantle signal is discussed in Sec. XC. This approach was also used in the optimization of the selection cuts, as mentioned in Sec. VII. The systematic uncertainties given in Sec. XI C are not included in the sensitivity studies and are only considered in the final results of Sec. XI. As it will be shown, the error on the geoneutrino measurement is largely dominated by the statistical error. TABLE XV. Summary of the expected number of events from nonantineutrino backgrounds in the antineutrino candidate sam­ple (exposure Ep 1.1.29 0.05.× 1032 protons × yr). The lim­its are 95% C.L. Background type Events 9Li background Untagged muons Fast n’s(µ in WT) Fast n’s(µ in rock) Accidental coincidences 3.6 1.0 0.023 0.007 <0.013 <1.43 3.846 0.017 (., n) in scintillator (., n) in buffer (.,n) Fission in PMTs 0.81 0.13 <2.6 <0.34 <0.057 214Bi - 214Po 0.003 0.001 Total 8.28 1.01 A. Geoneutrino analysis in a nutshell The geoneutrino signal is extracted from the spectral fit of the charges of the prompts of all selected IBD candi­dates. Since the number NIBD of selected candidates is relatively small (in this analysis, NIBD 1154 candidates, see Sec. XI A), an unbinned likelihood fit is used: NIBD Y ... L 1..; Qp.1 L..; Qi .; .37. pi11 where Q.p is the vector of individual prompt charges Qip, and index i runs from 1 to NIBD. 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 back­ground, 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 pullterms in the like­lihood 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 232Th and 238U contributions are summed and weighted according to the chondritic mass ratio of 3.9, corresponding to Rs signal ratioof0.27(Sec. VB). Alternatively, 232Th and 238U 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. XI E), 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. XI B 1) with the relatively­well known prediction of reactor antineutrino signal (Sec. VC). 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. VE), as it will be discussed in Sec. XI G. Typically, we include the following nonantineutrino backgrounds in the fit: cosmogenic 9Li, accidental coinci­dences, and (., n) interactions. Atmospheric neutrinos are included in the calculation of systematic uncertainties, as itwill be described in Sec. XI C. These background compo­nents 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 nonantineutrino backgrounds and in Table XIV given for atmospheric neutrinos. B. 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. VII). This so-called sensitivity study can be divided in the following four steps: (i) 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. VIII A, 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 nonantineutrino backgrounds in Table XV. (ii) The generated spectra are fit in the same way as the data (Sec. XA), using in the fit the same PDFs that were used for the generation of these pseudoexperi­ments. 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. (iii) The procedure is repeated 10 000 times for each configuration. In each pseudoexperiment, the number of generated events for signal and all backgrounds is varied according to the statistical uncertainty. (iv) The distributions of ratios of the resulting fit value (estimated) over the MC-truth (generated) 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 estimated from the fit over the number of generated geoneutrinos 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. C. Expected sensitivity Using the sensitivity tool as explained in Sec. XB, the expected statistical uncertainty of the Borexino geo­neutrino 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 232Th=238U 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 232Th=238U (considering the statistical fluctuations), the ratio resulting from the fit has nearly a flat distribution for the 10 000 pseudoexperiments. This will be also confirmed by large 232Th versus 238U contours shown in Fig. 48(d) for the fit of the data with free 238U and 232Th components. The sensitivity of Borexino to measure the mantle signal was studied using the log-likelihood ratio method [135] 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 pseu­doexperiments with the mantle geoneutrino component included. In addition, we have generated 1.2 million pseudoexperiments without the mantle contribution. In each dataset, we have included the relatively-well known lithospheric contribution (Table VI), as well as the reactor antineutrino “without 5 MeV excess” (Table XIV) and nonantineutrino backgrounds (Table XV). Each pseudoexperiment from all five datasets (one without the mantle and four with mantle signal according to four geological models), are fit twice: with and without the mantle contribution. The best fit with the mantle contribution fixed to zero corresponds to the likelihood Lf0g. The fit with the mantle component left free results in the likelihood Lfµg. Obviously, for the dataset without the mantle being generated, the two likelihoods tend to be the same. For the datasets with the mantle included, the Lf0g tends to be worse than Lfµg: the bigger this difference, the better the sensitivity to observe the mantle signal. We define the test statistics q (q . 0): q 1-2.ln Lf0g- ln Lfµg.; .38. that we call q0 for the dataset without the mantle generated. The q0 and the four q distributions for different geological models are shown in Fig. 45. The q0 corresponds to the theoretical f.qj0. distribution: 111 f.qj0.1 ..q.. p ........ exp - q: .39. 2 22.q 2 The four q distributions we fit with the f.qjµ. µ f.qjµ.1 1 - . ..q. . 2 11 p ... µ .p ........ exp - q- ; .40. 22.q 2 . where . stands for a cumulative Gaussian distribution with mean µ and standard deviation .. For high statistical M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Estimated/ Generated geoneutrino events Estimated/ Generated geoneutrino events FIG. 43. Results of the sensitivity study for geoneutrinos considering the conditions of the presented analysis: the PDFs for the ratio Estimated/Generated geoneutrino events for the CC, GC, GD, and FR Earth models and correspoding fits. Each PDF is based on 10 000 generated spectra. The Gaussian fits are all centered at one, so no systematic bias is expected; their .’s vary from .13.76 0.10.% to .23.09 0.17.%, depending on the expected geoneutrino signal for different BSE models (Table XIV) and represent the expected statistical uncertainty of the measurement. µ significance, µ=. is very large and ....› 1. In Figure 45 we show also qmed 1.µ=..2, the median value of f.qjµ., 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: Z . p 1 f.qj0.: .41. qmed The differences in qmed values shown in Fig. 45 for the 4 geological models, that correspond to different p-values and different sensitivity to observe the mantle signal, are to be ascribed only to the differences in the central values of the expected signals (Table XIV), which in turn come from measure 232Th=238U geoneutrino signal ratio. Solid line shows the different central values of U and Th masses associated to the different models (Table VII). the distribution of this ratio, assuming its chondritic value, for 10 000 generated pseudoexperiments. Distribution of this ratio, as The qobs from the data fit should be used to obtain the obtained from the fit (dotted line), is nearly flat, with a clear peak final statistical significance of the mantle signal, which will at 0, due to the 232Th contribution railed to 0. be described in Sec. XI E. Cosmochemical model Geochemical model 1 1 10-1 10-1 10-2 10-2 = 0.5594 med q p-value = 2.273 x 10-1 Probability/ 0.1 Probability/ 0.1 Probability/ 0.1Probability/ 0.1 10-3 10-4 10-3 10-4 10-5 0 5 10152025 q = - 2(ln L{0} - ln L{µ}) q = - 2(ln L{0} - ln L{µ}) Geodynamical model Fully radiogenic model 10-5 1 1 10-1 10-2 10-3 10-4 10-1 10-2 10-3 10-4 10-5 10-5 q = - 2(ln L{0} - ln L{µ}) q = - 2(ln L{0} - ln L{µ}) FIG. 45. Distributions of the test statistics q 1 -2.ln Lf0g - ln Lfµg., where Lf0g and Lfµg are likelihoods of the best fits obatined with the mantle contribution fixed to zero and left free, respectively. Brown solid represents the q values obtained from 10 000 pseudoexperiments with the generated mantle geoneutrino signal, based on the predictions of the different geological models (CC, GC, GD, and FR) and fit with f.qjµ. according to Eq. (40). The dark blue points show q 1 q0 test statistics obtained using 1.2M pseudoexperiments without any generated mantle signal and fit with f.qj0. [Eq. (39)]. The vertical dashed lines represent the medians qmed of the q distributions. The corresponding p-values are also shown. From the top left to the bottom right panels, the qmed values are increasing because of increasing expected mantle signal (i.e., increasing predicted U and Th masses in the mantle). XI. RESULTS This section describes the results of our analysis. In Sec. XI A the final IBD candidates selected with the optimized selection cuts are presented. In Sec. XI B the analysis, and in particular the spectral fit with the 238U=232Th ratio fixed to the chondritic value (Sec. XI B 1)or left free (Sec. XI B 2), is described. The systematic uncertainties are discussed in Sec. XI C. A summary of the geoneutrino signal as measured at the LNGS is given in Sec. XI D. Considering the expected signal from the bulk lithosphere (Table VI), we estimate the geoneutrino signal from the mantle in Sec. XI E. The consequences with regard to the Earth radiogenic heat are then presented in Sec. XI F. Finally, in Sec. XI G the constraints on the power of a hypothetical georeactor (Sec. VE) are set. A. Golden candidates In the period between December 9, 2007 and April 28, 2019, corresponding to 3262.74 days of data acquisition, NIBD 1 154 golden IBD candidates were observed to pass the data selection cuts described in Sec. VII. 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 IV of the prompt signal was also studied. This test would be particularly sensitive to a potential background originated from the IV 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. In fact, in a deformed IV, the points characterized by the same distance from the IV (and thus a potential source of background) can correspond to different radii. 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. B. Analysis An unbinned likelihood fit, as described in Sec. XA, was performed with the prompt charge of the 154 golden M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) Time distribution of golden candidates Radial distribution of golden candidates 50 45 40 35 30 15 10 10 5 5 0 (a) (b) Charge distribution of prompt golden candidates Charge distribution of delayed golden candidates Events/ 246 ton × yr 25 20 500 1000 1500 2000 2500 3000 3500 4000 1000 1500 2000 2500 3000 [p.e.] pQ [p.e.] dQ (c) (d) FIG. 46. Distributions for 154 golden IBD candidates (black data points). MC distributions (blue solid lines) are all normalized to the same number of events. (a) Observed IBD-rate (per average FVand one year) as a function of time in one year bins (December 2007 is included in 2008 data point). The dashed line represents the average IBD rate in the entire dataset. (b) Radial distribution of the prompt signals compared to the MC expectation. (c) Charge distribution of the prompts compared to the MC, assuming the geoneutrino and reactor antineutrino events follow the expectations as in Table XIV. (d) Charge distribution of the delayed compared to MC. The two peaks due to the captures on proton and on 12C are clearly visible. candidates shown in Sec. XI A. The three major non­antineutrino backgrounds, namely, the cosmogenic 9Li 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 antineu­trinos were unconstrained in the fit, using the PDF as in Fig. 32(c). The differences in the shape of the reactor antineutrino spectra “with 5 MeV excess” and “without 5 MeV excess” [bottom right in Fig. 32(d)] are included in the systematic uncertainty calculation (Sec. XI C). Obviously, geoneutrinos were also kept unconstrained. The fit was performed in two different ways with respect to the relative ratio of the 232Th and 238U contributions, as detailed in the next two subsections. The presented fit results are obtained following the recommendations given under the statistics chapter of [136] 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 like­lihood 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. 1. 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(a). The resulting spectral fit is shown in All events -0.50 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Prompt distance to IV [m] (a) Events with Q < 1500 p.e. Events with Q > 1500 p.e. pp Prompt distance to IV [m] Prompt distance to IV [m] (b) (c) FIG. 47. (a) Prompt’s distance to the IV for the 154 golden candidates (black data points) compared to MC (solid blue line) scaled to the same number of events. (b) and (c) show the same distribution but split in two energy windows: below and above the endpoint of the geoneutrinos at 1500 p.e. Fig. 48(a) and the numerical results are summarized in Table XVII. The likelihood profile for the number of geoneutrinos Ngeo [Fig. 49(a)], yields the best fit value Nbest 151.9, the median value Nmed 152.6, and the 68% coverage interval I68stat 1144.0–62.0 events. The likeli­ geo geo Ngeo hood profile for the number of reactor antineutrinos is shown in Fig. 49(b): Nbest 192.5 events was obtained with rea the median value Nmed 193.4 and the 68% coverage rea interval I68stat 1182.6–104.7 events. This is compatible Nrea 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. Thus, from the total of 154 golden IBD candidates, the number of detected antineutrinos (geo .reactor) is Nbest antinu 1144.4 events. This leaves the number of back­ground events compatible with the expectation (Table XV). The contour plot for Ngeo versus Nrea is shown in Fig. 48(c). The fit was also performed by constraining the expected number of reactor antineutrino events to (97.61.7) events (Table XIV). Theresult(thebestfitvalue Nbest Nmed geo 151.3, the median geo 152.0, and the 68% coverage interval I68stat 1143.6–61.1 events) is nearly Ngeo 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 endpoint well constrain the reactor antineu­trino contribution also in the geoneutrino window. The fact that without any constraint on Nrea the fit returns a value compatible with expectation is an important con­firmation of the Borexino ability to measure electron antineutrinos. M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 500 1000 1500 2000 2500 3000 3500 Qp [p.e.] Qp [p.e.] (a)(b) Nrea [events] NU [events] (c) (d) FIG. 48. Results of the analysis of 154 golden IBD candidates. (a) Spectral fit of the data (black points with Poissonian errors) assuming the chondritic Th/U ratio. The total fit function containing all signal and background components is shown in brownish-grey. Geoneutrinos (blue) and reactor antineutrinos (yellow) were kept as free fit parameters. Other nonantineutrino backgrounds were constrained in the fit. (b) Similar fit as in (a) but with 238U (dark blue) and 232Th (cyan) contributions as free and independent fit components. (c) The best fit point (black dot) and the contours for the 2D coverage of 68, 99.7, .100–5.7 × 10-5.%, and .100–1.2 × 10-13.%, (corresponding to 1, 3, 5, and 8., respectively), for Ngeo versus Nrea assuming Th/U chondritic ratio. The vertical lines mark the 1. bands of the expected reactor antineutrino signal (solid—without “5 MeVexcess,” dashed—with “5 MeVexcess”). For comparison, the star shows the best fit performed assuming the 238U and 232Th contributions as free and independent fit components. (d) The best fit (black dot) and the 68, 95.5, and 99.7% coverage contours (corresponding to 1., 2., and 3. contours) NTh versus NU. The dashed line represents the chondritic Th/U ratio. 2. Th and U as free fit parameters The second type of fit was performed by treating 238U and 232Th contributions as free and independent fit com­ponents. The corresponding MC PDFs from Fig. 32(b) were used. The spectral fit is shown in Fig. 48(b) and the numerical results are summarized in Table XVII. The 232Th likelihood profiles for the number of 238U and geoneutrinos are shown in Figs. 49(c) and 49(d), respec­ 127.8, Nmed tively. The fit yielded Nbest 129.0, and the UU 68% coverage interval I68NUstat 1116.1–43.1 events for the 121.1, Nmed Uranium contribution and Nbest 121.4, and Th Th the 68% coverage interval I68stat 1112.2–30.8 events for NTh the Th contribution. The best fit leads to 48.9 geoneutrinos in total, which is fully compatible with 51.9 geoneutrino events obtained in the case when Th/U ratio was fixed to the chondritic value. The only difference is significantly larger error in case of the fit with free U and Th contributions. Nbest I68stat For reactor antineutrinos, 195.8 and 1 rea Nrea 185.2–109.0 events were obtained, which is also compat­ible with the expectation. The total number of detected Nbest antineutrinos (geo .reactor) is antinu 1144.7 events. The contour plot for Ngeo versus Nrea is shown in Fig. 48(c). The contour plot for NU versus NTh is shown in Fig. 48(d). The results obtained after constraining the expected Nrea were again fully compatible with the results obtained when leaving the reactor antineutrino component free and without any significant reduction on error. ×10-3 ×10-3 0.45 0.35 0.4 0.3 0.35 0.25 Likelihood L Likelihood L Likelihood L Likelihood L 0.3 0.25 0.2 0.2 0.15 0.15 0.1 0.1 0.05 0.05 0 0 30 40 50 60 70 80 90 60 70 80 90 100110120130140 Ngeo [events] Nrea [events] (a) (b) ×10-3 ×10-3 0.45 0.3 0.4 0.25 0.35 0.1 0.05 0.05 0 0 0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 NU [events] NTh [events] (c) (d) FIG. 49. The likelihood profiles for the number of geoneutrino events Ngeo (a) and reactor antineutrino events Nrea (b) obtained from the fit assuming the chondritic Th/U ratio. The lower row shows the likelihood profiles for the number of 238U (c) and 232Th (d) events obtained from the fit assuming the 238U and 232Th contributions as the two independent free fit components. In each plot, the vertical solid red line indicates the best fit, while the vertical solid black and green lines indicate the median and mean values of the distributions, respectively. The vertical dashed/dotted lines show the 68%/99.7% confidence intervals of the distributions, corresponding to the 0.3 0.2 0.25 0.2 0.15 0.15 0.1 signal values given in Table XVII. C. 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. 1. Atmospheric neutrinos Atmospheric neutrinos as the source of background were discussed in Sec. VD, 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%. In addition, there is an indication of some over-estimation of this background, since above the endpoint of the reactor spectrum, where we would expect (3.31.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 TABLE XVI. Summary of the different sources of systematic uncertainty in the geoneutrino and reactor antineutrino measure­ment. Different contributions are summed up as uncorrelated. Source Geo error [%] Reactor error [%] .0.00 .0.00 Atmospheric neutrinos -0.38 -3.90 .0.00 .0.04 Shape of reactor spectrum -0.57 -0.00 .3.46 .3.25 Vessel shape -0.00 -0.00 Efficiency 1.5 1.5 Position reconstruction 3.6 3.6 .5.2 .5.1 Total -4.0 -5.5 M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 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 events from atmos­pheric neutrino background was constrained according to the values in Table XIV. First, the fit was performed for data up to the endpoint of the reactor antineutrino spectrum at 4000 p.e., which is the interval containing 63% of atmos­pheric neutrino background. The resulting number of atmospheric neutrino events is (4.63.2) and is compat­ible with the expectation of (6.73.4). The geoneutrino signal is almost unchanged, while Nrea decreased to (89.0 11.3) events. Second, we have performed the fit up to the endpoint 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.24.1) events. Fortunately, the resulting Ngeo and Nrea are nearly unchanged. To summarize, we estimate the respective systematic uncer­tainty on geoneutrinos as .0.00 % and on reactor antineu­ -0.38 trinos as .0.00 %. -3.90 2. Shape of the reactor spectrum The likelihood fit, as described in Sec. XI B 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. VC. In order to study the changes that might arise due to the observed “5MeV 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 Nrea and the two spectral shapes are relatively similar, the change in Nrea is very small: we observe an increase of 0.05 events. In case of Ngeo, we observe a decrease of 0.3 events. 3. Inner vessel shape reconstruction We consider a conservative 5 cm error on the IV position (Sec. III C). This means that the function defining our DFV (dIV 1 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 . 1 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 Ngeo and Nrea 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 % and on -0.00 reactor antineutrinos as .3.25 %. -0.00 4. MC efficiency The major source of uncertainty for the MC efficiency arisesfrom theevent lossesclose to theIVedges,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 241Am - 9Be calibration source compared to MC simula­tions for different source positions was studied. The uncertainty in the efficiency was then set to a conservative limit of 1.5%. 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 222Rn and 241Am - 9Be sources placed at 182 and 29 positions in the scintilla­tor, 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 geo­neutrino analysis was calculated using the shift in the positions for the 241Am - 9Be source. The maximal result­ing 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. D. Geoneutrino signal at LNGS This section details the conversion of the number of geoneutrino events Ngeo, resulting from the spectral fits described in Sec. XI B, to the geoneutrino signal Sgeo expressed in TNU, the unit introduced in Sec. VB: Ngeo Ngeo Sgeo1TNU 11 E0 ; .42. Epp .geo · 1032 1032 where the detection efficiency .geo 1 0.8698 0.0150 (Table XIII) and the exposure Ep 1.1.29 0.05.× 1032 protons×yr (Sec. IIIA). We obtain Sbest 146.3 TNU, geo the median value Smed 1 47.0 TNU, and including the geo systematic uncertainties from Table XVI, the 68% coverage interval I68full 1.38.9–55.6.TNU. This results in a final Sgeo precision of our measurement of .18.3% with respect to -17.2 Smed geo . The comparison of the result, obtained assuming the FIG. 50. Comparison of the expected geoneutrino signal Sgeo.U .Th. at LNGS (calculated according to different BSE models, see Sec. VB) with the Borexino measurement. For each model, the LOC and FFL contributions are the same (Table VI), while the mantle signal is obtained considering an intermediate scenario [Fig. 16(b)]. The error bars represent the 1. uncertain­ties of the total signal S.U .Th.. The horizontal solid back line represents the geoneutrino signal Smed, while the grey band the geo I68full interval as measured by Borexino. Sgeo chondritic Th/U mass ratio of 3.9, with the expected geo­neutrino signal considering different geological models (Sec. VB)is showninFig. 50.Figure 51 shows the time evolution of the Borexino measurements of the geoneutrino signal Sgeo.U .Th.at LNGS from 2010 up to the current result. 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 Sec. XI B. It was shown in Sec. XC 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. E. 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. VB and given in Table XIV as 28.8.5.5 events. The corresponding -4.6 MC PDF was constructed from the PDFs of 232Th and 238U 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 232Th and 238U PDFs, but the applied Th/U signal ratio was 0.26, the value discussed in Sec. VB. The mantle signal, as well as the reactor antineutrino contribution were free in the fit. The best fit is shown in Fig. 52(a). It resulted in a mantle signal of Nbest mantle 123.1 events, with the median value Nmed mantle 123.7 events, and the 68% coverage interval I68stat Nmantle 1.13.7–34.4.events. The likelihood profile of the mantle signal is shown in Fig. 52(b). After considering the systematic uncertainties, the final mantle signal can be Sbest given as with the median value mantle 120.6 TNU, Smed and the 68% coverage interval mantle 121.2 TNU, I68full TNU, as shown also in Table XVII. Smantle 1112.2–30.8 The statistical significance of the mantle signal was studied using MC pseudoexperiments with and without a generated mantle signal as described in Sec XC. The qobs obtained from the spectral fit is 5.4479, and it is compared with the theoretical function f.qj0., described in Sec. XC, Eq. (39), as shown in Fig. 53. 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). The Borexino measurement constrains at 90(95)% C.L. a mantle compo­sition with amantle.U.> 13.9.ppb and amantle.Th.> 48.34.ppb assuming for the mantle homogeneous distri­bution of U and Th and a Th/U mass ratio of 3.7. F. 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 HLSp rad .U .Th .K.1 8.1.-11..49 TW (Table V), corresponding to ~17% of the total terrestrial heat power Htot 1.47 2.TW. The radiogenic heat from the unexplored mantle could embrace a wide range of Hmantle .U.Th.K.1.1.2–39.8.TW (Table VII), rad 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 M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) TABLE XVII. Summary of the number of geoneutrino and reactor antineutrino events and the corresponding signals in TNU as well as the fluxes obtained from this work. The systematic uncertainties (Table XVI) are given for the median values. Criteria Best fit Median Mean 68% C.L. stat. 99.7% C.L. stat. Sys. error 68% C.L. stat. & sys. Rs 1S.Th.=S.U.10.27 Ngeo [events] 51.9 52.6 53.0 44.0–62.0 28.8–82.6 .2.7 -2.1 43.6–62.2 Nrea [events] 92.5 93.4 93.6 82.6–104.7 63.6–129.7 .4.8 -5.1 81.6–105.8 NU [events] 40.5 41.1 41.4 34.4–48.4 22.5–64.5 .2.1 -1.6 34.0–48.6 NTh [events] 11.4 11.5 11.6 9.7–13.6 6.3–18.1 .0.6 -0.5 9.6–13.7 Sgeo [TNU] 46.3 47.0 47.3 39.3–55.4 25.7–73.8 .2.4 -1.9 38.9–55.6 Srea [TNU] 79.7 80.5 80.7 71.2–90.3 54.8–111.8 .4.1 -4.4 70.3–91.2 SU [TNU] 35.9 36.3 36.6 30.4–42.8 19.9–57.1 .1.9 -1.5 30.1–43.0 STh [TNU] 10.4 10.5 10.6 8.7–12.4 5.7–16.5 .0.6 -0.4 8.8–12.6 .U [106 cm-2 s-1] 2.8 2.8 2.9 2.4–3.4 1.6–4.5 .0.2 -0.1 2.4–3.4 .Th [106 cm-2 s-1] 2.5 2.6 2.6 2.2–3.0 1.4–4.1 .0.2 -0.1 2.1–3.1 Rs.lithosphere.1S.Th.=S.U.10.29; Rs.mantle.1S.Th.=S.U.10.26 Nmantle [events] 23.1 23.7 24.1 13.7–34.4 0.6–57.2 .1.2 -1.0 13.6–34.4 Smantle [TNU] 20.6 21.2 21.5 12.2–30.7 0.5–51.1 .1.1 -0.9 12.2–30.8 S.Th.and S.U.independent Ngeo [events] 48.9 50.4 51.3 28.4–73.9 1.1–124.1 .2.6 -2.0 28.2–74.0 Nrea [events] 95.8 96.7 97.1 85.2–109.0 65.1–136.1 .4.9 -5.3 84.2–110.1 NU [events] 27.8 29.0 29.7 16.1–43.1 0.6–73.7 .1.5 -1.2 16.0–43.2 NTh [events] 21.1 21.4 21.6 12.3–30.8 0.5–50.4 .1.1 -0.9 12.2–30.8 Sgeo [TNU] 43.7 45.0 45.8 25.4–66.0 1.0–110.8 .2.3 -1.8 25.2–66.1 Srea [TNU] 82.6 83.4 83.7 73.5–94.0 56.1–117.3 .4.2 -4.6 72.6–94.9 SU [TNU] 24.6 25.7 26.3 14.3–38.1 0.5–65.2 .1.3 -1.0 14.2–38.2 STh [TNU] 19.2 19.5 19.6 11.2–28.0 0.5–45.8 .1.0 -0.8 11.1–28.0 .U [106 cm-2 s-1] 1.9 2.0 2.1 1.1–3.0 0.04–5.1 .0.1 -0.1 1.1–3.0 .Th [106 cm-2 s-1] 4.7 4.8 4.9 2.8–6.9 0.1–11.3 .0.3 -0.2 2.7–6.9 the dependence of the expected mantle geoneutrino signal Smantle.U .Th. on the mantle radiogenic power Hmantle rad .U .Th.. The unequivocal relation between the radiogenic power and the HPEs’ masses can be expressed via the constant U and Th specific heats h.U.1 98.5 µW=kg and h.Th.126.3 µW=kg [26]: Hmantle .U .Th. rad 1h.U.· Mmantle.U..h.Th.· Mmantle.Th. .43. 11h.U..3.7 · h.Th. · Mmantle.U.; where Mmantle.U.is the U mass in the mantle (Table VII). In the last passage as well as in all calculations below, we assume the mantle Th/U mass ratio of 3.7. With this assumption, for a given detector site the ratio: ß 1Smantle.U .Th.=Hmantle .U .Th..44. rad depends only on U and Th distribution in the mantle. For Borexino, the calculated ß ranges between ßlow 1 0.75 TNU=TW and ßhigh 10.98 TNU=TW, obtained assuming the HPEs placed in an unique HPEs-rich layer just above the CMB [i.e., low scenario, Fig. 16(a)] and homogeneously distributed in the mantle [high scenario, Fig. 16(c)], respectively. Considering then Eq. (43), the linear relation between the mantle signal Smantle.U .Th. and radiogenic power Hmantle .U .Th.can be expressed: rad Smantle.U .Th. 1ß · 1h.U..3.7 · h.Th. · Mmantle.U. 1ß · Hmantle .U .Th..45. rad is reported in Fig. 54, where the slope of central line (i.e., ßcentr 10.86 TNU=TW, black line) is the average of ßlow (blue line) and ßhigh (red line). The area between the two extreme lines denotes the region allowed by all possible 0 2 4 6 8101214 q = -2(ln L{0} - ln L{µ}) Qp [p.e.] 0 (a) FIG. 53. The q0 distribution for 1.2 million pseudoexperiments without any generated mantle signal fitted with f.qj0.[Eq. (39)]. ×10-3 The vertical dashed line represents qobs obtained from the data. 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0.4 The indicated p-value, calculated following Eq. (41) by setting qmed 1qobs, represents the statistical significance of the Borexino observation of the mantle signal. 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), which central value agrees with our measurement at 2.4. level. Likelihood L Nmantle [events] (b) FIG. 52. (a) Spectral fit to extract the mantle signal after constraining the contribution of the bulk lithosphere. The grey shaded area shows the summed PDFs of all the signal and background components. (b) The likelihood profile for Nmantle, the number of mantle geoneutrino events. The vertical solid red line indicates the best fit, while the vertical solid black and green lines indicate the median and mean values of the distributions, respectively. The vertical dashed/dotted lines represent the 68%/ 99.7% confidence intervals of the distribution. 30 mantle(U+Th) [TW] Hrad U and Th distributions in the mantle, assuming that the abundances in this reservoir are radial, non-decreasing function of the depth and in a fixed ratio Mmantle.Th.= Mmantle.U.13.7. 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 among the different BSE models can be studied in the space Smantle.U .Th.vs Hmantle .U .Th.. In Figure 54, the solid rad black horizontal line represents the Borexino measurement, the median Smed mantle, which falls within prediction of the Geodynamical model (GD). The 68% coverage interval I68full Smantle, also represented in Fig. 54 by horizontal black FIG. 54. Mantle geoneutrino signal expected in Borexino as a function of U and Th mantle radiogenic heat: the area between the red and blue lines denotes the full range allowed between a homogeneous mantle [high scenario–Fig. 16(c)] and a unique rich layer just above the CMB [low scenario–Fig. 16(a)]. The slope of the central inclined black line (ßcentr 10.86 TNU=TW) is the average of the slopes of the blue and red lines. The blue, green, red, and yellow ellipses are calculated with the following U and Th mantle radiogenic power Hmantle .U .Th. (with 1. rad error) according to different BSE models: CC model (3.10.5) TW, GC model (9.51.9) TW, GD model (21.32.4) TW, and FR model (32.21.4) TW. For each model darker to lighter shades of respective colors 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). M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) The mantle signal measured by Borexino can be con­verted to the corresponding radiogenic heat by inverting the Eq. (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 radiogenic power from U and Th in the mantle Hmantle .U .Th. inferred from the Borexino signal rad Smantle.U .Th.can be obtained with: Hmantle rad .U .Th.1.1=ßcentr.· Smantle.U .Th. .46. 11.16 · Smantle.U .Th.: Adopting Smed mantle.U .Th.121.2 TNU together with the 68% C.L. interval including both statistical and systematic errors (Table XVII), we obtain: Hmantle-med .U .Th.124.6 TW rad I68full .U .Th.114.2–35.7 TW: .47. Hmantle rad Summing the radiogenic power of U and Th in the litho­sphere HLSp TW, the Earth’s radiogenic rad .U .Th.16.9.-11..26 power from U and Th is Hrad.U .Th.131.7.14.4 TW. -9.2 Assuming the contribution from 40K to be 18% of the total mantle radiogenic heat (Sec. II), the total radiogenic and Hmantle genic heat HLSp , respectively, and secular cooling rad rad HSC (blue). The labels on the x axis identify different BSE models (Table VII), while the last bar labeled BX represents the Borexino measurement. The lithospheric contribution HLSp 1 rad 8.1.TW (Table V) is the same for all bars. The amount of -1.4 HPEs predicted by BSE models determines the mantle radiogenic heat (Table VII), while for Borexino the value of 30.0.13.5 TW is -12.7 inferred from the extracted mantle signal. The difference between Htot and the respective total radiogenic heat is assigned to the heat from secular cooling of the Earth. mantle signal can be expressed as Hmantle .U .Th .K.1 rad 30.0.13.5 TW, where we have expressed the 1. errors with -12.7 respect to the median. If we further add the lithospheric HLSp contribution rad .U .Th .K.18.1-.11..49 TW, we get the 68% coverage interval for the Earth’s radiogenic heat Hrad.U .Th .K.138.2.13.6 TW, as shown in Fig. 55. -12.7 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. 55. 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 URCV 10.78.0.41 is -0.28 compared to the URCV predicted by different BSE models in Fig. 56. The Borexino geoneutrino measurement constrains at 90(95)% C.L. a mantle radiogenic heat power to be Hmantle .U .Th.> 10.7.TW and Hmantle .U .Th . rad rad K.> 12.2.8.6.TW and the convective Urey ratio URCV > 0.13.0.04.. G. Testing the georeactor hypothesis The georeactor hypothesis described in Sec. VE was tested by performing the spectral fit after constraining the expected number of reactor antineutrino events (Table XIV) to 97.61.7.stat. 5.2.syst.. The geoneutrino (Th/U fixed to chondritic mass ratio of 3.9) and georeactor contributions were left free in the fit. For each georeactor rad 16.8.-11..14 The blue, green, and red colors represent different BSE models (CC, GC, and GD; Table VII, respectively). Ngeorea [events] FIG. 57. Likelihood profiles for the number of georeactor events Ngeorea obtained in the spectral fits with the constrained number of reactor antineutrino events. The four profiles represent the fits that differ in the shape of the PDF used for the georeactor contribution (Fig. 34), corresponding to different source positions in the deep Earth [GR1, GR2, GR3, Fig. 21(a)] and to a no­oscillation hypothesis. The vertical dashed lines indicate the upper limit at 95% C.L. for Ngeorea as obtained in each of the four fits. In setting the upper limit on the power of the georeactor, that depends on the assumed location of the georeactor, we use conservatively the highest limit on Ngeorea equal to 21.7 events. location [Fig. 21(a)], we have used the respective georeac­tor 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. 57. The vertical lines represent the 95% C.L. limits for the number of georeactor events Ngeorea obtained in fits with different georeactor PDFs. In setting the upper limit on the power of the georeactor, that depends on the assumed location of the georeactor, we use conservatively the highest limit on Ngeorea equal to 21.7 events. The latter is transformed to the signal Sgeorea 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. Therefore, we exclude the existence of a georeactor with a power greater than 0.5=2.4=5.7 TW at 95% C.L., assuming its location at 2900=6371=9842 km distance from the detector. XII. CONCLUSION 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 7Be solar neutrinos. In addition to solar neutrino measurements, Borexino has proven to be able to detect also antineutrinos. 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 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 provides an in-depth motivation and description of geoneutrino measurement, as well as the geological inter­pretations 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. 232Th Borexino detects geoneutrinos from 238U and through inverse beta decay, in which electron flavour antineutrinos with energies above 1.8 MeV interact with free protons of the LS. The detection efficiency for optimized data selection cuts is (87.01.5)%. This inter-action 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)× 1032 protons × year. This represents an increase by a factor of two over the previous Borexino analysis reported in 2015. 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 discrimi­nation 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 .stat..2.7 .sys. geoneutrinos (68% -8.6 -2.1 interval) from 238U and 232Th, a geoneutrino signal of 47.0.8.4 .stat..2.4 .sys. TNU has been obtained. The total -7.7 -1.9 precision of .18.3% is found to be in agreement with the -17.2 M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) 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 238U and 232Th were both fit as free parameters. Importance of the knowledge of abundances and dis­tributions 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 mea­sured geoneutrino signal is found to be in agreement with the predictions of different geological models with a preference for those predicting the highest concentrations of heat producing elements. 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 thick, U and Th depleted sediments. We note that geophysical and geochemical observations constrain the Th=Umass 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.5 .stat..1.1 .sys. TNU. -9.0 -0.9 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 TW (68% interval). Assuming -10.4 the contribution of 18% from 40K in the mantle and adding the relatively-well known lithospheric radiogenic heat of 8.1.1.9 TW, Borexino has estimated the total radiogenic -1.4 heat of the Earth to be 38.2.13.6 TW. The latter is found to -12.7 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 . In conclusion, Borexino geoneutrino meas­ -0.28 urement has constrained at 90% C.L. the mantle compo­sition to amantle (U) > 13 ppb and amantle (Th) > 48 ppb, the mantle radiogenic heat power to Hmantle (U . Th) > rad 10 TW and Hmantle .U . Th . K. > 12.2 TW, as well as rad the convective Urey ratio to URCV > 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.4TW has been excluded at 95% C.L. In conclusion, Borexino confirms the feasibility of geoneutrino measurements as well as the validity of differ­ent 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 scin­tillator detectors has a strong potential to provide funda­mental information about our planet. ACKNOWLEDGMENTS The Borexino program is made possible by funding from Istituto Nazionale di Fisica Nucleare (INFN) (Italy), National Science Foundation (NSF) (USA), Deutsche Forschungsgemeinschaft (DFG) and Helmholtz-Gemeinschaft (HGF) (Germany), Russian Foundation for Basic Research (RFBR) (Grants No. 16-29-13014ofi-m, No. 17-02-00305A, and No. 19-02-00097A) and Russian Science Foundation (RSF) (Grant No. 17-12-01009) (Russia), and Narodowe Centrum Nauki (NCN) (Grant No. UMO 2017/26/M/ST2/00915) (Poland). We acknowl­edge the hospitality and support of the Laboratori Nazionali del Gran Sasso in Italy. The authors thank Matteo Alb`eri, Kassandra Raptis, and Andrea Serafini for their support. APPENDIX: LIST OF ACRONYMS . Alpha particle A BSE model Anderson, 2007 [54] ß Beta particle BDT Boosted decision tree BSE Bulk silicate Earth BTB Borexino trigger board BTB4 The same as MTB flag, see below CC Continental crust CC model Cosmochemical bulk silicate Earth model C.L. Confidence level CLM Continental lithospheric margin CMB Core-mantle boundary CT Central Tile DAQ Data acquisition DFV Dynamical fiducial volume DM Depleted mantle DMP Dimethylphthalate (DMP, C6H4.COOCH3.2) e- or ß- Electron e. or ß. Positron EM Enriched mantle Ep Energy of the prompt IBD candidate Ed Energy of the delayed IBD candidate FADC Flash analog-to-digital converter FEB Front end board FFL Far field lithosphere FR model Fully radiogenic bulk silicate Earth model FWFD Fast wave form digitizer . Gamma ray G Gatti parameter G4Bx2 GEANT4 based Borexino Monte Carlo code GC model Geochemical bulk silicate Earth model GD model Geodynamical bulk silicate Earth model GR1, GR2, 3 studied positions of georeactor inside the Earth GR3 Hrad Earth’s radiogenic heat (Table continued) (Continued) HCC rad Earth’s continental crust radiogenic heat Hmantle rad Earth’s mantle radiogenic heat HLSp rad Earth’s lithosphere radiogenic heat HSC Earth’s heat from the secular cooling Htot Integrated total surface heat flux of the Earth HPEs Heat producing elements HSc High scenario of the mantle signal prediction IBD Inverse beta decay ID Inner detector IDF Inner detector flag Inner vessel ISc Intermediate scenario of the mantle signal prediction J BSE model Javoy et al., 2010 [34] LF Load factor of nuclear power plants L & K BSE model Lyubetskaya & Korenaga, 2007 [52] LNGS Laboratori Nazionali del Gran Sasso LOC Local crust LS Liquid scintillator LSc Low scenario of the mantle signal prediction LSp Lithosphere µ Muon MC Monte Carlo MLP Multilayer perceptron M & S BSE model McDonough & Sun, 1995 [47] MTB Muon trigger board MTF Muon trigger flag m w.e. Meter water equivalent . Neutrino . —Antineutrino .—e Electron flavor antineutrino n Neutron (Table continued) (Continued) Nh NP Npe OC OD OV p Pee PC PDF p.e. PID PM PMNS PMTs PPO P&O Qp Qd RR SSS SVM T TMVA TNU T&S URCV W WE WT Number of detected hits Number of triggered PMTs Number of detected photoelectrons Oceanic crust Outer detector Outer vessel Proton Survival probability of electron flavor neutrino Pseudocumene liquid scintillator, C6H3.CH3.3, 1,2,4-trimethylbenzene Probability distribution function Photoelectron(s) Particle identification Primitive mantle Pontecorvo–Maki–Nakagawa–Sakata mixing matrix Photomultiplier tubes Fluorescent dye, C15H11NO, 2,5-diphenyloxazole BSE model Palme and O’Neil, 2003 [56] Charge of the prompt IBD candidate Charge of the delayed IBD candidate Rest of the region Stainless steel sphere Support vector machine BSE model Taylor, 1980 [53] Toolkit for multivariate data analysis Terrestrial neutrino unit BSE model Turcotte & Schubert, 2002 [57] Convective Urey ratio BSE model Wang et al., 2018 [55] Water extraction procedure of LS-purification Water tank [1] W. Winter, Nucl. Phys. B908, 250 (2016). [2] A. Donini, S. Palomares-Ruiz, and J. Salvado, Nat. Phys. 15, 37 (2019). [3] S. Bourret, J. A. B. Coelho, and V. Van Elewyck (KM3NeT Collaboration), J. Phys. Conf. Ser. 888, 012114 (2017). [4] IAEA nuclear data section, IAEA database, https://www­ nds.iaea.org/relnsd/NdsEnsdf/QueryForm.html. [5] G. Fiorentini, M. Lissia, and F. Mantovani, Phys. Rep. 453, 117 (2007). [6] G. Marx and N. Menyhárd, Mitt. Sternwärte Budapest 44, 1 (1960). [7] G. Eder, Nucl. Phys. 78, 657 (1966). [8] G. Marx, Czech. J. Phys. B 19, 1471 (1969). [9] L. M. Krauss, S. L. Glashow, and D. N. Schramm, Nature (London) 310, 191 (1984). [10] C. G. Rothschild, M. C. Chen, and F. P. Calaprice, Geophys. Res. Lett. 25, 1083 (1998). [11] R. S. Raghavan, S. Schnert, S. Enomoto, J. Shirai, F. Suekane, and A. Suzuki, Phys. Rev. Lett. 80, 635 (1998). [12] T. Araki et al., Nature (London) 436, 499 (2005). [13] S. Abe et al. (KamLAND Collaboration), Phys. Rev. Lett. 100, 221803 (2008). [14] A. Gando et al. (KamLAND Collaboration), Nat. Geosci. 4, 647 (2011). [15] A. Gando et al. (KamLAND Collaboration), Phys. Rev. D 88, 033001 (2013). [16] G. Bellini et al. (Borexino Collaboration), Phys. Lett. B 687, 299 (2010). [17] G. Bellini et al. (Borexino Collaboration), Phys. Lett. B 722, 295 (2013). [18] M. Agostini et al. (Borexino Collaboration), Phys. Rev. D 92, 031101 (2015). [19] N. Vinyoles, A. M. Serenelli, F. L. Villante, S. Basu, J. Bergstr, M. C. Gonzalez-Garcia, M. Maltoni, C. Pea Garay, and N. Song, Astrophys. J. 835, 202 (2017). [20] G. Bellini et al. (Borexino Collaboration), Phys. Rev. D 89, 112007 (2014). [21] G. Bellini et al. (Borexino Collaboration), Nature (London) 512, 383 (2014). M. AGOSTINI et al. PHYS. REV. D 101, 012009 (2020) [22] M. Agostini et al. (Borexino Collaboration), Nature (London) 562, 505 (2018). [23] M. C. Chen, Earth Moon Planets 99, 221 (2006). [24] F. An et al. (JUNO Collaboration), J. Phys. G 43, 030401 (2016). [25] L. Wan, G. Hussain, Z. Wang, and S. Chen, Phys. Rev. D 95, 053001 (2017). [26] S. Dye, Rev. Geophys. 50, RG3007 (2012). [27] M. Agostini et al. (Borexino Collaboration), Astropart. Phys. 97, 136 (2018). [28] R. M. Gaschnig, R. L. Rudnick, W. F. McDonough, A. J. Kaufman, J. Valley, Z. Hu, S. Gao, and M. L. Beck, Geochim. Cosmochim. Acta 186, 316 (2016). [29] National Research Council et al., Origin and Evolution of Earth: Research Questions for a Changing Planet (The National Academies Press, Washington, DC, 2008). [30] W. F. McDonough, in Treatise on Geochemistry, 2nd ed., edited by H. D. Holland and K. K. Turekian (Elsevier, New York, 2014), Vol. 3, pp. 559–577. [31] D. C. Rubie, D. J. Frost, U. Mann, Y. Asahara, F. Nimmo, K. Tsuno, P. Kegler, A. Holzheid, and H. Palme, Earth Planet. Sci. Lett. 301, 31 (2011). [32] G. Bellini, A. Ianni, L. Ludhova, F. Mantovani, and W. F. McDonough, Prog. Part. Nucl. Phys. 73, 1 (2013). [33] A. Deuss and J. Woodhouse, Science 294, 354 (2001). [34] M. Javoy et al., Earth Planet. Sci. Lett. 293, 259 (2010). [35] Y. Huang, V. Chubakov, F. Mantovani, R. L. Rudnick, and W. F. McDonough, Geochem. Geophys. Geosyst. 14, 2003 (2013). [36] N. A. Simmons, A. M. Forte, L. Boschi, and S. P. Grand, J. Geophys. Res. 115, B12310 (2010). [37] J. Korenaga, Rev. Geophys. 46 (2008). [38] C. Jaupart, S. Labrosse, F. Lucazeau, and J. C. Mareschal, in Treatise on Geophysics, Vol. 7, edited by G. Schubert (Elsevier B.V., New York, 2015), pp. 223–270. [39] A. M. Hofmeister and R. E. Criss, Tectonophysics 395, 159 (2005). [40] A. G. Crosby, D. McKenzie, and J. G. Sclater, Geophys. J. Int. 166, 553 (2006). [41] C. Jaupart and J. C. Mareschal, in Treatise on Geophysics, edited by G. Schubert (Elsevier B.V., New York, 2007), Vol. 6, pp. 217–251. [42] J. H. Davies and D. R. Davies, Solid Earth 1, 5 (2010). [43] D.L.WilliamsandR.P.VonHerzen, Geology 2, 327 (1974). [44] G. F. Davies, Rev. Geophys. 18, 718 (1980). [45] J. G. Sclater, C. Jaupart, and D. Galson, Rev. Geophys. 18, 269 (1980). [46] H. N. Pollack, S. J. Hurter, and J. R. Johnson, Rev. Geo­ phys. 31, 267 (1993). [47] W. F. McDonough and S. Sun, Chem. Geol. 120, 223 (1995). [48] A. Bouvier and M. Boyet, Nature (London) 537, 399 (2016). [49] A. Wohlers and B. J. Wood, Geochim. Cosmochim. Acta 205, 226 (2017). [50] I. Blanchard, J. Siebert, S. Borensztajn, and J. Badro, Geochem. Perspect. Lett. 5, 1 (2017). [51] B. A. Chidester, Z. Rahman, K. Righter, and A. J. Campbell, Geochim. Cosmochim. Acta 199, 1 (2017). [52] T. Lyubetskaya and J. Korenaga, J. Geophys. Res. 112, B03212 (2007). [53] S. Taylor, Proc. Lunar Planet. Sci. Conf. 11, 333 (1980). [54] D. L. Anderson, New Theory of the Earth, 2nd ed. (Cambridge University Press, Cambridge, England, 2007). [55] H. S. Wang, C. H. Lineweaver, and T. R. Ireland, Icarus 299, 460 (2018). [56] H. Palme and H. O’Neill, in Treatise on Geochemistry, edited by H. D. Holland and K. K. Turekian (Elsevier, New York, 2003), Vol. 2, pp. 1–38. [57] D. L. Turcotte and G. Schubert, Geodynamics, Applications of Continuum Physics to Geological Problems, 2nd ed. (Cambridge University Press, Cambridge, England, 2002). [58] O.Šrámek,W. F.McDonough,E. S.Kite,V.Lekić,S. T.Dye, and S. Zhong, Earth Planet. Sci. Lett. 361, 356 (2013). [59] F. Mantovani, L. Carmignani, G. Fiorentini, and M. Lissia, Phys. Rev. D 69, 013001 (2004). [60] R. Arevalo, W. F. McDonough, and M. Luong, Earth Planet. Sci. Lett. 278, 361 (2009). [61] J. Blichert-Toft, B. Zanda, D. S. Ebel, and F. Albar´ede, Earth Planet. Sci. Lett. 300, 152 (2010). [62] S. A. Wipperfurth, M. Guo, O. Šrámek, and W. F. McDonough, Earth Planet Sci. Lett. 498, 196 (2018). [63] V. Strati, S. A. Wipperfurth, M. Baldoncini, W. F. McDonough, and F. Mantovani, Geochem. Geophys. Geosyst. 18, 4326 (2017). [64] Y. Huang, V. Strati, F. Mantovani, S. B. Shirey, and W. F. McDonough, Geochem. Geophys. Geosys. 15, 3925 (2014). [65] M. Coltorti et al., Geochim. Cosmochim. Acta 75,2271 (2011). [66] J. M. Herndon, J. Geomagn. Geoelectr. 45, 423 (1993). [67] J. M. Herndon, Proc. Natl. Acad. Sci. U.S. A. 93, 646 (1996). [68] V. D. Rusov, V. N. Pavlovich, V. N. Vaschenko, V. A. Tarasov, T. N. Zelentsova, V. N. Bolshakov, D. A. Litvinov, S. I. Kosenko, and O. A. Byegunova, J. Geophys. Res. 112, 203 (2007). [69] R. Meijer and W. Westrenen, South Afr. J. Sci. 104, 111 (2008). [70] Y. Wang and L. Wen, J. Geophys. Res. 109, B10305 (2004). [71] B. Szczerbinska, D. Alyssa, and M. Dongming, Proc. Sci. 90, 13 (2011). [72] A. Cabrera et al., arXiv:1908.02859. [73] K. A. Goettel, Geophys. Surv. 2, 369 (1976). [74] G. Alimonti et al. (Borexino Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 600, 568 (2009). [75] M. Agostini et al. (Borexino Collaboration), J. Cosmol. Astropart. Phys. 02 (2019) 046. [76] H. Back et al. (Borexino Collaboration), J. Instrum. 7, P10018 (2012). [77] G. Bellini et al. (Borexino Collaboration), J. Instrum. 6, P05005 (2011). [78] H. Voss, A. Hocker, J. Stelzer, and F. Tegenfeldt, Proc. Sci., ACAT2007 (2007) 040. [79] C. Cortes and V. Vapnik, Mach. Learn. 20, 273 (1995). [80] Y. Freund and R. E. Schapire, J. Comput. Syst. Sci. 55, 119 (1997). [81] G. Lukyanchenko, Ph.D. thesis, National Research Center Kurchatov Institute, (in Russian) (2017), http://www.nrcki .ru/files/pdf/1490777684.pdf. [82] H. Back et al. (Borexino Collaboration), Nucl. Instrum. Methods Phys. Res., Sect. A 584, 98 (2008). [83] M. Agostini et al. (Borexino Collaboration), Astropart. Phys. 92, 21 (2017). [84] A. Strumia and F. Vissani, Phys. Lett. B 564, 42 (2003). [85] F. Capozzi, E. Lisi, and A. Marrone, Phys. Rev. D 89, 013001 (2014). [86] I. Esteban et al., Nufit 3.2, NuFIT website, http://www.nu­ fit.org/?q=node/166 (2018). [87] M. Baldoncini, I. Callegari, G. Fiorentini, F. Mantovani, B. Ricci, V. Strati, and G. Xhixha, Phys. Rev. D 91, 065002 (2015). [88] A. M. Dziewonski and D. L. Anderson, Phys. Earth Planet. Interiors 25, 297 (1981). [89] G. Fiorentini, A. Ianni, G. Korga, M. Lissia, F. Mantovani, L. Miramonti, L. Oberauer, M. Obolensky, O. Smirnov, and Y. Suvorov, Phys. Rev. C 81, 034602 (2010). [90] S. Enomoto, Geoneutrino spectrum and luminosity, online version https://www.awa.tohoku.ac.jp/~sanshiro/research/ geoneutrino/spectrum/. [91] CROP Project: Deep Seismic Exploration of the Central Mediterranean and Italy, edited by I. Finetti (Elsevier Science, New York, 2005), Vol. 1. [92] T. Plank, in Treatise of Geochemistry, 2nd ed., edited by H. D. Holland and K. K. Turekian (Elsevier, New York, 2014), Vol. 4, pp. 607–629. [93] C. Bassin, G. Laske, and G. Masters, EOS Trans AGU 81, F897 (2000). [94] G. Laske, G. Masters, and C. Reif, Crust 2.0. A new global crustal model at 2 × 2 degrees, CRUST 2.0 website, https://igppweb.ucsd.edu/~gabi/crust2.html (2001). [95] N. Shapiro and M. Ritzwoller, Geophys. J. Int. 151,88 (2002). [96] M. Reguzzoni and D. Sampietro, Int. J. Appl. Earth Observ. Geoinf. 35, 31 (2015). [97] G. Laske and T. Masters, EOS Trans. AGU 78, F483 (1997). [98] World Nuclear Association Database, “Reactor database,” Database website, https://world-nuclear.org/information­ library/facts-and-figures/reactor-database.aspx?id=27232 (2015). [99] X. B. Ma, W. L. Zhong, L. Z. Wang, Y. X. Chen, and J. Cao, Phys. Rev. C 88, 014605 (2013). [100] IAEA, Power reactor information system (PRIS), PRIS webpage, www.iaea.org/pris (2019). [101] F. P. An et al. (Daya Bay Collaboration), Phys. Rev. Lett. 116, 061801 (2016); 118, 099902(E) (2017). [102] Y. Abe et al. (Double Chooz Collaboration), J. High Energy Phys. 10 (2014) 086; 02 (2015) 074(E). [103] M. Y. Pac (RENO Collaboration), Proc. Sci., NU­FACT2017 (2018) 038 [arXiv:1801.04049]. [104] K. Siyeon (NEOS Collaboration), Proc. Sci., ICRC2017 (2018) 1024. [105] T. A. Mueller et al., Phys. Rev. C 83, 054615 (2011). [106] F. P. An et al. (Daya Bay Collaboration), Phys. Rev. Lett. 118, 251801 (2017). [107] G.Mention,M.Vivier,J.Gaffiot,T.Lasserre,A.Letourneau, and T. Materna, Phys. Lett. B 773, 307 (2017). [108] H. Ejiri, J. Suhonen, and K. Zuber, Phys. Rep. 797,1 (2019). [109] M. Honda, M. S. Athar, T. Kajita, K. Kasahara, and S. Midorikawa, Phys. Rev. D 92, 023004 (2015). [110] G. Battistoni, A. Ferrari, T. Montaruli, and P. R. Sala, Astropart. Phys. 23, 526 (2005). [111] R. Wendell, Prob3++ software for computing three flavor neutrino oscillation probabilities, Prob3++ website (2012). [112] C. Andreopoulos et al., Nucl. Instrum. Methods Phys. Res., Sect. A 614, 87 (2010). [113] K. Farley and E. Neroda, Annu. Rev. Earth Planet Sci. 26, 189 (1998). [114] S. T. Dye, Phys. Lett. B 679, 15 (2009). [115] V. D. Rusov et al., J. Mod. Phys. 4, 528 (2013). [116] J. M. Herndon and D. A. Edgerley, arXiv:hep-ph/0501216. [117] G. Bellini et al. (Borexino Collaboration), J. Cosmol. Astropart. Phys. 08 (2013) 049. [118] H. de Kerret et al. (Double Chooz Collaboration), J. High Energy Phys. 11 (2018) 053. [119] J. H. Kelley, J. E. Purcell, and C. G. Sheu, Nucl. Phys. A968, 71 (2017). [120] P. Mohr, Phys. Rev. C 97, 064613 (2018). [121] S. Harissopulos, H. W. Becker, J. W. Hammer, A. Lagoyannis, C. Rolfs, and F. Strieder, Phys. Rev. C 72, 062801 (2005). [122] A. J. Koning and D. Rochman, Nucl. Data Sheets 113, 2841 (2012). [123] G. Bellini et al. (Borexino Collaboration), Eur. Phys. J. A 49, 92 (2013). [124] P. Vogel and J. F. Beacom, Phys. Rev. D 60,053003 (1999). [125] Y. Prezado et al., Phys. Lett. B 618, 43 (2005). [126] D. M. Mei and A. Hime, Phys. Rev. D 73, 053004 (2006). [127] S. Westerdale, NeuCBOT (Neutron Calculator Based On TALYS), GitHub, https://github.com/shawest/neucbot. [128] S. Westerdale and P. D. Meyers, Nucl. Instrum. Methods Phys. Res., Sect. A 875, 57 (2017). [129] S. Westerdale (DEAP-3600 Collaboration), AIP Conf. Proc. 1921, 060002 (2018). [130] S. S. Westerdale, A study of nuclear recoil backgrounds in Dark Matter detectors, Ph. D. thesis, Princeton University, 2016, online version: https://doi.org/10.2172/1350520. [131] A. Koning, S. Hilaire, and M. Duijvestijn, Talys-1.9, TALYS website, http://www.talys.eu/download-talys/. [132] TALYS-1.0 Proceedings of the International Conference on Nuclear Data for Science and Technology–ND2007, edited by O. Bersillon, F. Gunsing, E. Bauge, R. Jacqmin, and S. Leray (EDP Sciences, Nice, France, 2008). [133] D. W. McKee, J. K. Busenitz, and I. Ostrovskiy, Nucl. Instrum. Methods Phys. Res., Sect. A 587, 272 (2008). [134] G. Alimonti et al., Nucl. Instrum. Methods Phys. Res., Sect. A 406, 411 (1998). [135] G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C 71, 1554 (2011); 73, 2501(E) (2013). [136] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, 030001 (2018). Correction: The order of Figures 14 and 15 was presented incorrectly and has been fixed. The previously published Figure 14(a) contained the wrong image and has been set right.