NEUTRINO MASS ORDERING STUDIES INCLUDING NEUTRINO/ANTINEUTRINO DISCRIMINATION WITH ICECUBE-DEEPCORE By Sebastian Enrique Sanchez Herrera A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics โ€“ Doctor of Philosophy 2022 ABSTRACT The observation of neutrino flavor oscillations by independent measurements of the Sudbury Neutrino Observatory (SNO) and the Super Kamiokande projects demonstrated the existence of neutrino masses, providing the first evidence of physics beyond the standard model (SM). To complete the picture of the neutrino oscillation formalism, one of the most important open questions is the ordering of the three neutrino mass eigenstates ๐‘š 1 , ๐‘š 2 and ๐‘š 3 . After incredible efforts by solar, atmospheric and accelerator neutrino experiments, this problem is now reduced to a simpler question: whether the third mass state ๐‘š 3 is above or below the other two. The former case corresponds to the so called normal ordering (NO) whereas the latter is called inverted ordering (IO). So far, there is not enough sensitivity to the neutrino mass ordering (NMO) to reach a definitive conclusion on this key oscillation parameter. Sensitivity to the NMO in atmospheric neutrino os- cillation experiments arises from distinct matter effects in the oscillation probabilities for neutrinos in NO or antineutrinos in IO, after propagation through the matter profile of the Earth. The other two combinations: antineutrinos in NO and neutrinos in IO, present negligible matter effects in their oscillation profiles. The IceCube neutrino observatory has performed a measurement on this parameter using the DeepCore detector, with results that are consistent with both orderings. A main challenge for NMO measurements with atmospheric neutrinos in IceCube is that sensitivity to the NMO is maximal when comparing detector expectations for NO/IO to experimental data for clean samples of neutrinos or antineutrinos, but this separation is typically not possible for very large-scale Cherenkov detectors. However, sensitivity to the NMO is achieved through statistical differences in the total neutrino + antineutrino count rates in the detector between orderings, arising from differences in the cross sections and atmospheric fluxes for neutrinos and antineutrinos. This work presents a study of the neutrino mass ordering with DeepCore that aims to directly separate ๐œˆ ๐œ‡ and ๐œˆยฏ ๐œ‡ , thus incorporating the oscillation matter effects to full potential sensitivity. The novel algorithm developed leverages the difference in signatures of muon decay for ๐œˆ ๐œ‡ and ๐œˆยฏ ๐œ‡ within the natural ice sheet. Within the sparsely instrumented DeepCore detector array, the ability to directly detect these low-energy signatures and identify the decay lifetimes was implemented. Dominant backgrounds in the DeepCore optical modules, including artifacts in the operation of the photomultiplier tubes (dark noise and afterpulsing) and the onboard electronics, were studied in detail and improved techniques, including improved models, to mitigate the effects were incorpo- rated in the study. A boosted decision tree was trained to optimize the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ separation power with the other key input parameter (neutrino energy, direction) of the mass ordering analysis, resulting in a modest improvement in the signature sensitivity (0.358+0.060 โˆ’0.057 ๐œŽ compared to 0.226+0.072 โˆ’0.068 ๐œŽ). A binary ensemble test utilizing 7.5 years of DeepCore data demonstrated a small preference for the โ€œnormalโ€ neutrino mass ordering over the inverted case with p-values ๐‘ NO = 0.72 and ๐‘ IO = 0.17. The ๐ถ ๐ฟ ๐‘  value of the result, including the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ separation parameter, is 0.61, consistent with both mass ordering scenarios. In the future, the fully developed analysis will be enhanced with data from the IceCube Upgrade which will incorporate an array design that will make it possible to enhance the detectability of the muon decay signature photons while significantly reducing the impact of internal optical module backgrounds. Copyright by SEBASTIAN ENRIQUE SANCHEZ HERRERA 2022 This thesis is dedicated to my mom and sister. v ACKNOWLEDGEMENTS I want to start by thanking my mother for all the opportunities she has allowed me to have and for guiding me throughout my life. It is because of her unbelievable efforts that I was able to pursue a career in physics. I also thank my sister for her constant advice and support, for motivating me unconditionally throughout the entirety of my PhD and for being my closest and dearest friend. Next I thank Darren Grant and Claudio Kopper for being outstanding advisors and mentors. Darren, I am extremely grateful for your guidance, invaluable insights and support. You always managed to identify the best course of action to complete my research tasks and I thank you for keeping me on track towards my academic and professional goals. Claudio, thank you for constantly helping me overcome the challenges that arised during the completion of my PhD. Many thanks for the countless conversations about the IceCube detector, software tools, particle physics and statistics. In addition, I thank my committee members Wade Fisher, Artemis Spyrou and Tyler Cocker for their excellent guidance and support and I also thank Kim Crosslan for all her help and patience in helping me identify and satisfy all requirements for the completion of my PhD. Finally, I want to thank the IceCube collaboration as a whole. I was able to learn an incredible amount of professional and personal skills while working for the collaboration and I truly enjoyed the time I spent working on my research project. Thanks to all IceCube researchers, professors, postdocs and students, in particular to those within the local research groups at Michigan State University and University of Alberta. Special thanks to Roger Moore, Tyce DeYoung, Kirsten Tollefson, Tom Stuttard, Ben Jones, Phillip Eller, Chris Weaver, Rob Halliday, Brian Clark, Shiqi Yu, Mehr Un Nisa, Joshua Hignight, Juan Pablo Yanez, Michael Larson, Jan Weldert, Sarah Nowicki, Tania Wood and Alexander Trettin. I thank you all for sharing your knowledge, wisdom and expertise, allowing me to successfully complete my research goals. vi PREFACE The work presented hereby is the product of original work by its Author except where noted below. Work from other authors is cited where relevant. The IceCube Neutrino Observatory is the result of the combined efforts of more than 300 scientists and collaborators around the world who designed, constructed, operate and perform the data analyses. It is typical for an author to have contributed to several projects and participated in multiple collaborative efforts. This work presents the first exploration of neutrino/antineutrino discrimination using muon decay signals within the IceCube collaboration, as well as the first effort within the collaboration of a neutrino oscillation analysis that included neutrino/antineutrino separation. Nevertheless, several of the applied methods and software tools were the product of the efforts and work of other collaborators. Specifically, the selection criteria for the atmospheric neutrino sample used in this work is the result of collaborative efforts by the oscillations working group within the IceCube collaboration, and information regarding the details of this work is currently only available within internal collaboration documents [1, 2]. The Monte Carlo sample used in this work was also produced by other collaborators within the oscillations working group, but the muon neutrino sample was reprocessed [by the author] for the purpose of this work, by adding a key physics process required for neutrino/antineutrino discrimination using muon decay signals, muon capture, to the simulation chain. The neutrino mass ordering study performed in this work was based on a previous neutrino mass ordering analysis performed by the collaboration using three years of DeepCore data [3], but updated to included neutrino/antineutrino discrimination. The methodology to measure the neutrino mass ordering required tools to calculate expected counts in the detector given a set of neutrino oscillation parameters and use these expectations to fit for the underlying oscillation parameters in simulated or experimental data. This software is called PISA and is outlined in [4]. Finally, the methodology used in this work to account for systematic uncertainties in the fits was also defined by the oscillation working group in [2], prior to this work. vii TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 The neutrino mass ordering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Towards an improved NMO sensitivity with IceCube-DeepCore (thesis overview) . 3 CHAPTER 2 PHYSICS BACKGROUND . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 The Standard Model of Particle Physics . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Neutrino oscillations in vacuum . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Global fits to the oscillation parameters . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Neutrino oscillations in matter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.5 Atmospheric neutrinos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.6 Neutrino interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.7 Physics background for neutrino/antineutrino discrimination . . . . . . . . . . . . 33 CHAPTER 3 THE ICECUBE AND DEEPCORE DETECTORS . . . . . . . . . . . . . . 41 3.1 Detection principle for Cherenkov detectors . . . . . . . . . . . . . . . . . . . . . 41 3.2 The IceCube Detector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 Neutrino event topologies in IceCube . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4 The Antarctic ice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.5 Background and noise signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 CHAPTER 4 NEUTRINO/ANTINEUTRINO CLASSIFICATION . . . . . . . . . . . . . 58 4.1 Muon decay signals in DeepCore . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.2 Neutrino/Antineutrino classifier . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 CHAPTER 5 THE DEEPCORE NEUTRINO MASS ORDERING ANALYSIS . . . . . . 74 5.1 Analysis method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.2 Asimov sensitivities to the Neutrino Mass Ordering . . . . . . . . . . . . . . . . . 92 5.3 Measurement of the Neutrino Mass Ordering using 7.5 years of DeepCore data . . 98 CHAPTER 6 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 APPENDIX A DATA ACQUISITION AND ONLINE PROCESSING IN ICECUBE . . 112 APPENDIX B THE ICECUBE UPGRADE . . . . . . . . . . . . . . . . . . . . . . . 120 APPENDIX C AFTERPULSES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 APPENDIX D DEEPCORE NMO ANALYSIS - EVENT SIMULATIONS AND WEIGHTING . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 APPENDIX E DEEPCORE NMO ANALYSIS - ANALYSIS SAMPLE . . . . . . . . . 139 APPENDIX F RECONSTRUCTION ALGORITHMS . . . . . . . . . . . . . . . . . . 153 viii APPENDIX G SANTAโ€™S TRACK AND CASCADE FITS . . . . . . . . . . . . . . . . 160 ix CHAPTER 1 INTRODUCTION Neutrinos were first hypothesised in 1930 by Pauli [5], in an attempt to explain the continuous energy spectrum of electrons produced in the ๐›ฝ decay process. Observed was a two body decay from a nuclei ๐‘1 to a lighter nuclei ๐‘2 of the form: A N โ†’ Z+1 AN + ๐‘’ โˆ’ , where charge was Z 1 2 evidently conserved. However, the expected final energy of the decay secondaries in a two body decay is fixed. To explain the continuous spectrum, Pauli proposed the existence of an unobserved, electrically neutral particle. This particle would receive its modern name, โ€œneutrinoโ€, when Fermi presented a formal theory of beta decay [5]. Evidence for the existence of neutrinos increased as weak decay processes were further studied. In particular, muon decay, where the observed process was ๐œ‡โˆ’ โ†’ ๐‘’ โˆ’ , once again showed a continuous electron energy spectrum, indicating the production of two invisible neutrinos. Neutrino interactions would not be directly observed until 1956, when (anti)neutrinos were detected via inverse beta decay processes (๐œˆยฏ + ๐‘ โ†’ ๐‘› + ๐‘’ + ) at the Hanford reactor experiment [6]. Soon after the direct observation of neutrinos, it was also hypothesised that there were two different neutrino types [5] following from the fact that the decay channel ๐œ‡ โ†’ ๐‘’ + ๐›พ was not observed. This implied a quantum number was not shared between the electron and the muon. The decay process ๐œ‡โˆ’ โ†’ ๐‘’ โˆ’ + ๐œˆ + ๐œˆ would imply one neutrino carried this quantum number (muon lepton number) and was therefore different from the neutrino observed in the beta decay process, that evidently didnโ€™t. In 1962, the first evidence of a two neutrino scenario was provided from measurements at Brookhaven [7]. In these measurements, antineutrinos produced in pion decay (๐œ‹ โˆ’ โ†’ ๐œ‡โˆ’ + ๐œˆ) ยฏ were detected, and it was shown that the reaction ๐œˆยฏ + ๐‘ โ†’ ๐‘› + ๐‘’ never occurred, whereas ๐œˆยฏ + ๐‘ โ†’ ๐‘› + ๐œ‡ was observed. Since then, the "muon neutrino" (๐œˆ ๐œ‡ ) and "electron neutrino" (๐œˆ๐‘’ ) were known to be different particles. Today three neutrino flavors have been detected: electron, muon and tau. The latter, however, would not be experimentally detected until the year 2000 [8]. The end of the 1960โ€™s witnessed a rapid evolution of the field with the rise of the so-called solar neutrino problem that described a deficit of the expected (๐œˆ๐‘’ ) neutrino flux produced in the Sun first 1 found by the Homestake experiment in 1968 [5]. This experiment was only sensitive to electron neutrinos through the interaction: 37๐ถ๐‘™ + ๐œˆ๐‘’ โ†’ 37 ๐ด๐‘Ÿ + ๐‘’ โˆ’ so a possible explanation proposed by Bruno Pontecorvo in the same year [5] was that as the neutrinos propagated to the Earth part of the electron neutrino flux was transformed into a different neutrino flavor. After more than 30 years, the Solar Neutrino Problem was solved by the Sudbury Neutrino Observatory (SNO), a heavy water (๐ท 2 ๐‘‚) Cerenkov detector located in Canada [9]. The experiment was sensitive to both the total neutrino flux and the individual electron neutrino flux and found that while the total flux was consistent with the expected flux of solar neutrinos, the electron neutrino flux only accounted for about a third. This demonstrated that the deficit in the Solar Neutrino Problem was a consequence of roughly two-thirds of the electron neutrinos produced at the Sun changing flavor after reaching the Earth [5]. The phenomena of neutrinos changing flavor as they propagate is now commonly referred as neutrino oscillations. The experimental confirmation of these oscillations is attributed both to the results obtained by SNO and by atmospheric neutrino measurements by Super Kamiokande [10], a water Cherenkov experiment located in Japan. In this experiment, muon and electron neutrinos produced within the Earthโ€™s atmosphere in air showers were detected. The experiment was capable of measuring the energy and directionality of the neutrino events, and to discriminate between ๐œˆ๐‘’ -like and ๐œˆ ๐œ‡ -like events. At multi-GeV energies, it was observed that neutrino fluxes arriving from above the detector (and thus traveling a relatively short distance) agreed with production model flux expectations, whereas measured fluxes corresponding to neutrinos traveling through the Earth (larger distances) presented a deficit of ๐œˆ ๐œ‡ events [11]. No deviation from expectations was observed for ๐œˆ๐‘’ events, and thus, the results were consistent with ๐œˆ ๐œ‡ โ†’ ๐œˆ๐œ oscillations. Neutrino oscillations are central to the work presented in this thesis which is to explore a novel method for measuring one of the key remaining unknown parameters of the process, the so-called neutrino mass ordering. 2 1.1 The neutrino mass ordering The discovery of neutrino oscillations by the SNO and Super Kamiokande demonstrated that these particles have a finite mass. An open question in particle physics is the ordering of the neutrino masses. Experimental measurements involving solar neutrino flavor oscillations have provided the ordering between the first and second states. Combined results from atmospheric and solar neutrino experiments have also provided measurements on the magnitude of the squared mass differences between the first and second, and second and third mass states. However, what is still unknown is whether the third neutrino mass state is heavier or lighter than the other two states. This is formally known as the neutrino mass ordering (NMO) problem. There are two possible mass orderings, normal ordering (NO), where ๐‘š 3 > ๐‘š 1,2 is satisfied, or inverted ordering (IO), where ๐‘š 3 < ๐‘š 1,2 . Figure 1.1 illustrates these two scenarios, and approximate ๐œˆ๐‘’ , ๐œˆ ๐œ‡ and ๐œˆ๐œ contributions of the neutrino mass states, based on current knowledge of neutrino oscillations. Current experimental sensitivities to the neutrino mass ordering based on global fits to the results of several neutrino oscillation experiments, show IO is currently only disfavoured at a 2 โˆ’ 3๐œŽ level [13]. 1.2 Towards an improved NMO sensitivity with IceCube-DeepCore (thesis overview) Neutrino oscillations for atmospheric neutrinos detected by IceCube-DeepCore are enhanced with respect to vacuum neutrino oscillation probabilities via matter effects, where the magnitude of this enhancement is determined by the density profile of the Earth. However, matter effects in the oscillation probabilities are visible only for neutrinos in the NO, or antineutrinos in the IO, for energies below โˆผ 20GeV. Thus, samples of atmospheric neutrinos or antineutrinos may be used to infer the NMO by looking for deviations from the expected vacuum oscillated event counts in the detector. The initial IceCube-DeepCore studies did not have the capability to directly separate neutrinos from antineutrinos, and instead utilized a combined ๐œˆ + ๐œˆยฏ sample. A larger compatibility 3 Figure 1.1: The neutrino normal ordering (left) and inverted ordering (right) [12]. Oscillation experiments are only sensitive to the mass square differences (๐‘š๐‘–2 โˆ’ ๐‘š 2๐‘— ) and therefore only these differences are illustrated. Also shown is the approximate flavor composition of the neutrino mass states. The absolute mass of the neutrino states is currently unknown. to matter oscillations was expected in the NO (where matter effects arise for neutrinos) than in the IO (where matter effects arise for antineutrinos). Completed an initial study of the NMO with 3 years of atmospheric neutrino events detected by DeepCore the measured sensitivity to the NMO was compatible with both orderings, and with a small observed preference for the NO below the 1๐œŽ level [3]. This proof of concept study demonstrated the potential to enhance the detectorโ€™s sensitivity. Presented in this body of work is the first steps towards realizing neutrino-antineutrino discrim- ination with IceCube-DeepCore, using muon decay signatures, to achieve a direct measurement of the NMO. An overview of the physics of neutrino oscillations (in both vacuum and matter), atmospheric neutrino oscillations, neutrino interactions, and a review on how neutrino and antineu- trinos can be separated via muon decay signals, is presented in chapter 2. Chapter 3 describes the IceCube and DeepCore detectors. Next Chapter 4 describes the classifier designed for sepa- ration of neutrino and antineutrino events developed for the analysis. Finally, the NMO analysis 4 with neutrino/antineutrino discrimination using DeepCore data is presented (Chapter 5) and final conclusions reported (Chapter 6). 5 CHAPTER 2 PHYSICS BACKGROUND The analysis presented in this work is a neutrino oscillation analysis, as discussed in chapter 1. Therefore, it is imperative to first introduce formally the main properties of neutrinos and the theoretical framework that describes their interactions. This framework, known as the Standard Model (SM), is presented in section 2.1 and it represents our current best description of fundamental particles and their interactions. Sections 2.2, 2.3 and 2.4, introduce the formalism of vacuum neutrino oscillations, global fits to the oscillation parameters and neutrino oscillations in matter, respectively. Since atmospheric neutrinos are the neutrino source used in the oscillation analyses with IceCube-DeepCore, including the NMO analysis presented in this thesis, their production mechanism will be covered in detail in section 2.5. Next, in section 2.6, the discussion will be focused towards the neutrino interactions that are relevant for the IceCube detector: nucleon- neutrino interactions. Finally, section 2.7 covers the main properties of muons that are important when studying their decay signals in matter and reviews why these signals can be used to separate ๐œˆ and ๐œˆยฏ populations. 2.1 The Standard Model of Particle Physics There are two types of particles within the SM. First there are bosons, integer spin particles that (with the exception of the Higgs boson, see below) act as force carriers in particle interactions. These include photons and gluons, that mediate electromagnetic (EM) and strong interactions respectively, and the W and Z bosons, responsible for the weak interaction. The second group of particles within the SM are called fermions, half odd integer spin particles that act as the building blocks of regular matter. These are further subdivided into quarks and leptons. Quarks come in three generation pairs: up (u) and down (d), charm (c) and strange (s), top (t) and bottom (b). The u, c, and t quark have positive charge (+2/3 e) while the d, s, and t have negative charges (-1/3 e). Particles made of quarks are referred to as hadrons, which include subgroups such as the mesons (quark-antiquark pairs) and baryons (three quark particles). In a similar fashion, leptons occur in three generations of particle pairs, where a charged lepton 6 Figure 2.1: A graphical representation of the Standard Model [14]. is always paired with a neutrino. A graphical representation of the SM, illustrating the different fermion generations as well as all bosons, can be found in figure 2.1. Within the SM, quarks undergo all types of interactions, whereas charged leptons canโ€™t par- ticipate in strong interactions and neutrinos are only able to interact via the weak force. When the force carrier interchanged is a (neutral) Z boson, the interaction is called a neutral current (NC) interaction. The Feynman diagram for these interactions involves an incoming and outgoing neutrino (see figure 2.2). On the other hand, weak interactions where a (charged) W boson is the mediator are called charged current (CC) interactions. In this type of interaction, a neutrino always couples to a charged lepton, which in turn defines its flavor. Therefore, neutrinos within the SM can have three flavors, electron (๐œˆ๐‘’ ), muon (๐œˆ ๐œ‡ ) or tau (๐œˆ๐œ ), defined by the charged lepton they couple to in CC interactions. Within the SM, there is also a scalar field called the Higgs field. Quarks, charged leptons and the W and Z bosons couple to this field, which gives them their masses (photons and gluons are massless). No mechanism within the SM satisfactorily accounts for neutrino masses, so within the SM, neutrinos are considered massless. However, given the definitive evidence for neutrino flavor oscillations (see chapter 1), neutrinos are now known to be massive particles, as discussed next. 7 ฮฝe,ยต,ฯ„ ฮฝe,ยต,ฯ„ ฮฝe,ยต,ฯ„ e, ยต, ฯ„ Z W Figure 2.2: Feynman diagrams for NC (left) and CC (right) weak interactions for neutrinos. Notice for CC interactions that neutrinos are always coupled to a charged lepton that defines their flavor. Diagrams produced with the LaTeX package TikZ [15]. 2.2 Neutrino oscillations in vacuum The discussion presented here is based on the descriptions provided in [5], [16], [11], [17] and [18]. The starting point when discussing neutrino oscillations is to note that since neutrino flavor states can change as they propagate, they cannot be eigenstates of the free (vacuum) propagation Hamiltonian, H. However, flavor states can be written in terms of the eigenstates of propagation (and vice-versa), commonly referred to as the mass eigenstates: โˆ‘๏ธ |๐œˆ๐›ผ โŸฉ = ๐‘ˆ๐›ผโˆ— ๐‘— ๐œˆ ๐‘— ๐‘— โˆ‘๏ธ (2.1) ๐œˆ๐‘— = ๐‘ˆ๐›ผ ๐‘— |๐œˆ๐›ผ โŸฉ , ๐›ผ where ๐‘ˆ๐›ผ ๐‘— is the change of basis matrix between the flavor and mass (propagation) eigenbases (and therefore unitary), ๐›ผ indicates the neutrino flavor (e, ๐œ‡ and ๐œ) and j denotes the mass states, named by convention 1, 2 and 3. A neutrino generated in a weak interaction will always be produced as a flavor eigenstate |๐œˆ๐›ผ โŸฉ. For free propagation, its time evolution will be defined by application of the (free) time evolution operator to this original state. Therefore, the neutrino state at a time t will be: โˆ’๐‘–๐ธ ๐‘— ๐‘ก โˆ‘๏ธ โˆ‘๏ธ |๐œˆ(๐‘ก)โŸฉ = ๐‘ˆ (๐‘ก) |๐œˆ๐›ผ โŸฉ = ๐‘’ โˆ’๐‘–๐ป๐‘ก ๐‘ˆ๐›ผโˆ— ๐‘— ๐œˆ ๐‘— = ๐‘’ ๐‘ˆ๐›ผโˆ— ๐‘— ๐œˆ ๐‘— , (2.2) ๐‘— ๐‘— 8 where ๐ธ ๐‘— is the energy of the j-th mass eigenstate. These energies can be approximated in terms of the neutrino masses by noticing that they are expected to be very small (neutrinos move at nearly the speed of light): โˆš๏ธƒ ๐‘š 2๐‘—  ๐ธ๐‘— = ยฎ2 | ๐‘| + ๐‘š 2๐‘— โ‰ˆ | ๐‘| ยฎ 1+ , (2.3) ยฎ2 2| ๐‘| where ๐‘ยฎ is the momentum of the neutrino state and ๐‘š ๐‘— is the mass of the j-th mass state. Now is possible to calculate the probability that the neutrino is found with a different flavor state ๐œˆ ๐›ฝ (๐›ผ โ‰  ๐›ฝ ) after traveling a distance L (โ‰ˆ ๐‘ก) in a vacuum: โˆ‘๏ธ โˆ’๐‘–๐‘š 2๐‘— ๐ฟ/2๐ธ โˆ’๐‘–| ๐‘|๐ฟ 2 ๐‘ƒ๐›ผโ†’๐›ฝ (๐ฟ) = | ๐œˆ ๐›ฝ ๐œˆ(๐ฟ) | 2 = ๐‘ˆ ๐›ฝ ๐‘— ๐‘ˆ๐›ผโˆ— ๐‘— ๐‘’ ๐‘’ ยฎ . (2.4) ๐‘— In the expression above, we have used the first order approximation | ๐‘| ยฎ โ‰ˆ ๐ธ, where E is the energy of the neutrino state. Notice that the term |๐‘’๐‘–| ๐‘|๐‘ก ยฎ | 2 will not contribute and, in the same fashion, no constant term in the Hamiltonian will contribute to the oscillation probabilities. 2.2.1 Two flavor oscillations Although in reality there are three known neutrino flavor states, is generally helpful to first introduce the main features of neutrino oscillations by simplifying the discussion to a two flavor case. In this case, there are only two flavor states ๐œˆ๐‘’ , ๐œˆ๐‘ฅ and two mass states ๐œˆ1 , ๐œˆ2 . Since in two dimensions, a unitary matrix is necessarily a 2x2 rotation matrix, parametrized by a single angle ๐œƒ, the change of basis matrix can be written as: ยฉ cos๐œƒ sin๐œƒ ยช ยฉ ๐‘ˆ๐‘’1 ๐‘ˆ๐‘’2 ยช ๐‘ˆ๐›ผ ๐‘— = ยญยญ ยฎ=ยญ ยฎ ยญ ยฎ. ยฎ (2.5) โˆ’sin๐œƒ cos๐œƒ ๐‘ˆ๐‘ฅ1 ๐‘ˆ๐‘ฅ2 ยซ ยฌ ยซ ยฌ Using this expression in equation 2.4, results in the following expression for the oscillation probabilities in the two flavor case: 2 โˆ’๐‘–๐‘š 2 ๐ฟ/2๐ธ โˆ’๐‘–๐‘š 2 ๐ฟ/2๐ธ ฮ”๐‘š 2 ๐ฟ  ๐‘ƒ๐œˆ๐‘’ โ†’๐œˆ๐‘ฅ (๐ฟ) = ๐‘ˆ๐‘ฅ1๐‘ˆ๐‘’1 ๐‘’ 1 + ๐‘ˆ๐‘ฅ2๐‘ˆ๐‘’2 ๐‘’ 2 = sin2 2๐œƒ sin2 , (2.6) 4๐ธ where ฮ”๐‘š 2 = ๐‘š 22 โˆ’ ๐‘š 12 . Several important properties of neutrino oscillations can be derived from equation 2.6: 9 โ€ข Oscillation probabilities in vacuum always depend on two fundamental properties of neutri- nos, the mixing angle between the flavor states (๐œƒ) and the squared difference between their masses (ฮ”๐‘š 2 ). The other quantities in equation 2.6 are the neutrino energy, that depends on their production mechanism, and their path length. These are the variables that can be measured by a neutrino detector when performing an oscillation analysis. โ€ข The amplitude of the oscillations is defined by the mixing angle. No oscillations occur at ๐œƒ = 0, and oscillations are enhanced as ๐œƒ โ†’ ๐œ‹/4, the so-called maximum mixing scenario. โ€ข The frequency of the oscillations depends on the squared mass difference. As a result, a requirement for oscillations is the existence of at least one neutrino state of finite mass. โ€ข Oscillations have no sensitivity to individual neutrino masses. In addition, vacuum oscilla- tions are not sensitive to the mass ordering (the sign of ฮ”๐‘š 2 ). In future chapters, neutrino oscillation studies with the IceCube experiment will be described. Those measurements are performed using neutrinos with typical path lengths of the order of thousands of kilometers and energies above 1 GeV. As a result ฮ”๐‘š 2 ๐ฟ ฮ”๐‘š 2 [eV2 ] ๐ฟ [km] [natural units] โ†’ 1.27 (2.7) 4๐ธ ๐ธ [GeV] is a useful conversion rule to keep for future discussion. 2.2.2 Three flavor oscillations Considering now the three flavor case, the change of basis (mixing) matrix has a special name, the PMNS (Pontecorvo-Maki-Nakawaga-Sakata) matrix [5]. For the purpose of neutrino oscillations, the PMNS matrix can be fully parametrized using three mixing angles ๐œƒ 13 , ๐œƒ 23 and ๐œƒ 12 and a complex phase ๐›ฟ๐ถ๐‘ƒ , which introduces CP violation into the lepton sector: ยฉ ๐‘ 12 ๐‘ 13 ๐‘ 12 ๐‘ 13 ๐‘ 13 ๐‘’ โˆ’๐‘–๐›ฟ๐ถ๐‘ƒ ยช ยญ ยฎ ๐‘ˆ๐›ผ ๐‘— = ยญ โˆ’๐‘ 12 ๐‘ 23 โˆ’ ๐‘ 12 ๐‘ 13 ๐‘ 23 ๐‘’๐‘–๐›ฟ๐ถ๐‘ƒ ๐‘ 12 ๐‘ 23 โˆ’ ๐‘ 12 ๐‘ 13 ๐‘ 23 ๐‘’๐‘–๐›ฟ๐ถ๐‘ƒ ยญ ยฎ ๐‘ 13 ๐‘ 23 ยฎ ยญ ยฎ ยญ ยฎ ๐‘  ๐‘  โˆ’ ๐‘ 12 ๐‘ 13 ๐‘ 23 ๐‘’๐‘–๐›ฟ๐ถ๐‘ƒ โˆ’๐‘ 12 ๐‘ 23 โˆ’ ๐‘ 12 ๐‘ 13 ๐‘ 23 ๐‘’๐‘–๐›ฟ๐ถ๐‘ƒ ๐‘ 13 ๐‘ 23 ยซ 12 23 ยฌ 10 ๏ฃฎ ๏ฃฏ1 0 0 ๏ฃน๏ฃฎ ๏ฃบ ๏ฃฏ ๐‘ 13 0 ๐‘  ๐‘’ โˆ’๐‘–๐›ฟ๐ถ๐‘ƒ ๏ฃน๏ฃบ ๏ฃฎ ๏ฃฏ ๐‘ 12 ๐‘ 12 0๏ฃบ ๏ฃน ๏ฃฏ ๏ฃบ๏ฃฏ 13 ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ = ๏ฃฏ0 ๐‘ 23 ๐‘ 23 ๏ฃบ ๏ฃฏ ๏ฃฏ ๏ฃบ ๏ฃฏ 0 1 0 ๏ฃบ ๏ฃบ ๏ฃฏโˆ’๐‘  ๏ฃฏ 12 ๐‘ 12 0๏ฃบ , ๏ฃบ (2.8) ๏ฃฏ ๏ฃบ๏ฃฏ ๏ฃบ ๏ฃฏ ๏ฃบ ๏ฃฏ0 โˆ’๐‘ 23 ๐‘ 23 ๏ฃบ ๏ฃฏโˆ’๐‘ 13 ๐‘’ ๏ฃฏ ๏ฃบ๏ฃฏ ๐‘–๐›ฟ ๐ถ๐‘ƒ 0 ๐‘ 13 ๏ฃบ ๏ฃบ ๏ฃฏ 0 ๏ฃฏ 0 1๏ฃบ ๏ฃบ ๏ฃฐ ๏ฃป๏ฃฐ ๏ฃป ๏ฃฐ ๏ฃป ๐ด๐‘ก๐‘š๐‘œ๐‘ ๐‘โ„Ž๐‘’๐‘Ÿ๐‘–๐‘ ๐‘…๐‘’๐‘Ž๐‘๐‘ก๐‘œ๐‘Ÿ ๐‘†๐‘œ๐‘™๐‘Ž๐‘Ÿ where ๐‘ ๐‘– ๐‘— = sin(๐œƒ๐‘– ๐‘— ) and ๐‘๐‘– ๐‘— = cos(๐œƒ๐‘– ๐‘— ). Notice the PMNS matrix can be written as the product of three unitary matrices, each described by a different mixing angle. Different neutrino experiments have larger sensitivity to different mixing angles, and equation 2.8 indicates the type of neutrino experiments that are customarily used to probe each. Oscillation probabilities will ultimately only depend on the squared mass differences, instead of the individual neutrino masses. As a result, we will have two additional free parameters in three flavor oscillations, since the third squared mass difference can always be found using the other two: 2 + ฮ”๐‘š 2 + ฮ”๐‘š 2 = 0 ฮ”๐‘š 32 (2.9) 21 13 With equations 2.4 and 2.8, the following expression is found for the vacuum oscillation probabilities in the three flavor case: โˆ‘๏ธ  ฮ”๐‘š๐‘–2๐‘— ๐ฟ  ๐‘ƒ๐›ผโ†’๐›ฝ (๐ฟ) = ๐›ฟ๐›ผ๐›ฝ โˆ’ 4 Re(๐‘ˆ๐›ผ๐‘–โˆ— ๐‘ˆ ๐‘ˆ ๐‘ˆ โˆ— ) sin2 ๐›ฝ๐‘– ๐›ผ ๐‘— ๐›ฝ ๐‘— 4๐ธ ๐‘–> ๐‘— (2.10) โˆ‘๏ธ  ฮ”๐‘š๐‘–2๐‘— ๐ฟ  +2 Im(๐‘ˆ๐›ผ๐‘–โˆ— ๐‘ˆ ๐‘ˆ ๐‘ˆ โˆ— ) sin , ๐›ฝ๐‘– ๐›ผ ๐‘— ๐›ฝ ๐‘— 2๐ธ ๐‘–> ๐‘— Again, the vacuum oscillation probabilities depend directly on the ratio L/E, that can be tuned by neutrino experiments performing oscillation analyses. However, there is now a small sensitivity to the sign of the squared mass splittings as long as there is non-zero CP violation in the lepton sector (๐›ฟ๐ถ๐‘ƒ โ‰  0). In other words, in the case of no CP violation, vacuum oscillations are insensitive to the NMO. Again, the amplitude of the oscillations is given by the mixing matrix, whereas the frequency of the oscillations depends on the squared mass splitting. The oscillation probabilities above apply to neutrino states. For anti-neutrinos, one needs to replace ๐‘ˆ by ๐‘ˆ โˆ— , equivalent to the following conversion: ๐‘ƒ( ๐›ผยฏ โ†’ ๐›ฝ) ยฏ = ๐‘ƒ(๐›ผ โ†’ ๐›ฝ; ๐›ฟ๐ถ๐‘ƒ โ†’ โˆ’๐›ฟ๐ถ๐‘ƒ ). (2.11) 11 As a result, vacuum oscillation probabilities are equal for neutrinos and antineutrinos in the absence of CP violation, as would be expected. 2.3 Global fits to the oscillation parameters 2 , ฮ”๐‘š 2 and To experimentally measure the six free oscillation parameters (๐œƒ 13 , ๐œƒ 23 , ๐œƒ 12 , ฮ”๐‘š 21 32 ๐›ฟ๐ถ๐‘ƒ ), two types of experiments can be designed. The first is a disappearance experiment, where the probabilities of survival โˆ‘๏ธ ๐‘ƒ(๐œˆ๐›ผ โ†’ ๐œˆ๐›ผ ) = 1 โˆ’ ๐‘ƒ๐›ผโ†’๐›ฝ , (2.12) ๐›ฝโ‰ ๐›ผ for a particular neutrino flavor are probed by measuring a deficit in the expected flux (of that flavor). The second is an appearance analysis, where a detector measures an excess of neutrinos of flavor ๐œˆ ๐›ฝ โ‰  ๐œˆ๐›ผ , where ๐œˆ๐›ผ is the flavor generated at the neutrino production site. An enormous interest followed after the discovery of oscillations by SNO [9] and Super Kamiokande [10], where several experiments set out to measure flavor oscillations and contribute to the experimental constraints of the oscillation parameters. A complete picture of these efforts is beyond the scope of this work. Instead, discussion will be restricted to the oscillation parameters that were ultimately used in the NMO analysis with ๐œˆ/๐œˆยฏ separation presented in this work, specif- ically the results of a global fit of neutrino oscillations data performed by the NuFit collaboration in 2018 [19, 20] (see table 2.1). Results on the mass splitting ฮ”๐‘š 21 2 and the mixing angle ๐œƒ were derived mainly from solar 12 neutrino experiments, including the measurements from the Homestake experiment [21], Super Kamiokande [22, 23, 24], Sage [25], Gallex [26], SNO [27] and Borexino [28, 29, 30], and the reactor neutrino experiment KamLAND [31]. For these experiments, the ratio ๐ฟ/๐ธ is quite large, so that, in general, 1 โ‰ช ๐ฟ|ฮ”๐‘š 3๐‘™ 2 |/๐ธ is   satisfied. This means that the oscillation terms sin2 ฮ”๐‘š 3๐‘™ 2 ๐ฟ/4๐ธ average out to 1/2 [11] and, as a result, the oscillation probabilities in equation 2.10 become independent of the ฮ”๐‘š 3๐‘™ 2 mass splittings. In fact, under those conditions, once the oscillation probabilities are evaluated (and assuming no CP violation for simplicity), the following expression is obtained for the ๐œˆ๐‘’ survival 12 Parameter NO best fit (ยฑ 1 ๐œŽ) IO best fit (ยฑ 1 ๐œŽ) 105 ยท ฮ”๐‘š 21 2 /eV2 7.39+0.21 โˆ’0.20 7.39+0.21 โˆ’0.20 103 ยท ฮ”๐‘š 3๐‘™ 2 /eV2 2.525+0.033 โˆ’0.031 โˆ’2.512+0.034 โˆ’0.031 ๐œƒ 12 /โ—ฆ 33.82+0.78 โˆ’0.76 33.82+0.78 โˆ’0.75 ๐œƒ 13 /โ—ฆ 8.61+0.12 โˆ’0.13 8.65+0.12 โˆ’0.13 ๐œƒ 23 /โ—ฆ 49.7+0.9 โˆ’1.1 +0.9 49.7โˆ’1.0 ๐›ฟ๐ถ๐‘ƒ /โ—ฆ 217+40 โˆ’28 +25 280โˆ’28 Table 2.1: Results from the global fit to the neutrino oscillation parameters performed by the NuFit collaboration in 2018, evaluating both mass orderings. The mass splitting difference ฮ”๐‘š 3๐‘™ 2 2 > 0 and ฮ”๐‘š 2 < 0 and for NO and IO, respectively [19, 20]. is defined as ฮ”๐‘š 31 32 probabilities [32]: h  ฮ”๐‘š 2 ๐ฟ i 21 ๐‘ƒ๐‘’โ†’๐‘’ (๐ฟ) = cos4 (๐œƒ 13 ) 1 โˆ’ sin2 (2๐œƒ 12 ) sin2 + sin4 (๐œƒ 13 ), (2.13) 4๐ธ which can be approximated as  ฮ”๐‘š 2 ๐ฟ  21 ๐‘ƒ๐‘’โ†’๐‘’ (๐ฟ) โ‰ˆ 1 โˆ’ sin2 (2๐œƒ 12 ) sin2 , (2.14) 4๐ธ since ๐œƒ 13 is small. It is clear, once we have performed these approximations, why these experiments would be sensitive to the 1-2 sector. In addition, notice we have recovered the results that would be expected from a two flavor mixing scenario (equation 2.6). The global fit values for ๐œƒ 13 arise primarily from data provided by reactor neutrino experiments, including Daya Bay [33], Reno [34], and Double Chooz [35], where ๐œˆยฏ ๐‘’ , produced in nuclear reactors via the beta-decay (๐‘› โ†’ ๐‘ + ๐œˆยฏ ๐‘’ + ๐‘’ โˆ’ ) of fission products [36], are detected using inverse beta decay (๐‘ + ๐œˆยฏ ๐‘’ โ†’ ๐‘’ โˆ’ + ๐‘›) to probe the oscillation probabilities associated to ๐œˆยฏ ๐‘’ disappearance. Lastly, constraints to ๐œƒ 23 and ฮ”๐‘š 3๐‘™ 2 were derived mainly from two different types of neutrino experiments. Accelerator neutrino experiments, including Minos [37, 38], Nova [39] and T2K [40], where muon neutrino beams are produced, and used to probe oscillation probabilities for ๐œˆ ๐œ‡ disappearance and ๐œˆ๐‘’ appearance for both neutrinos and antineutrinos. Atmospheric neutrino experiments, including Super Kamiokande [41] and the IceCube-DeepCore detector [42], are able to probe the oscillation probabilities via ๐œˆ ๐œ‡ disappearance. 13 Both atmospheric and accelerator neutrino experiments can probe L/E values for which ๐ฟ|ฮ”๐‘š 3๐‘™2 |/๐ธ โ‰ˆ ๐‘‚ (1), and therefore are sensitive to single oscillation periods in the 3-1 and 3-2 sectors. The direct consequence of this is they probe L/E values where ๐ฟ|ฮ”๐‘š 21 2 |/๐ธ โ‰ช 1, and thus are insensitive to ฮ”๐‘š 212 . Using the approximation ๐ฟ|ฮ”๐‘š 2 |/๐ธ โ‰ช 1 in equation 2.10, is possible to 21 obtain the following approximations to the ๐œˆ ๐œ‡ oscillation probabilities in vacuum (again assuming no CP violation for simplicity): 2 ๐ฟ ฮ”๐‘š 31 sin2 (๐œƒ ) sin2 (2๐œƒ ) sin2  ๐‘ƒ(๐œˆ ๐œ‡ โ†’ ๐œˆ๐‘’ ) โ‰ˆ 23 13 , (2.15) 4๐ธ 2 ๐ฟ ฮ”๐‘š 31 ๐‘ƒ(๐œˆ ๐œ‡ โ†’ ๐œˆ๐œ ) โ‰ˆ cos4 (๐œƒ 13 ) sin2 (2๐œƒ 23 ) sin2  . (2.16) 4๐ธ Notice how the main contribution comes from oscillations into tau neutrinos, since ๐œƒ 13 is small. The muon survival probabilities can then be approximated as: 2 ๐ฟ ฮ”๐‘š 31 ๐‘ƒ(๐œˆ ๐œ‡ โ†’ ๐œˆ ๐œ‡ ) โ‰ˆ 1 โˆ’ sin2 (2๐œƒ 23 ) sin2  , (2.17) 4๐ธ and once again, the approximation emulates the expectations for two flavor oscillations. The approximation in equation 2.17 to the vacuum oscillation probabilities holds quite well for the path lengths and energies studied in the analysis performed in this thesis, as shown in figure 2.3. According to equation 2.17 there is a degeneracy in the oscillation probabilities, so that: ๐‘ƒ(๐œƒ 23 ) โ‰ˆ ๐‘ƒ(90โ—ฆ โˆ’ ๐œƒ 23 ) (2.18) for ๐œƒ 23 โˆˆ [0, 90] โ—ฆ . Equation 2.18 only holds approximately, and corrections from the ๐œˆ๐‘’ appearance term (equation 2.15), CP violation and matter effects, will break the symmetry. However, the symmetry is present in the leading term in the oscillation probabilities, and disentanglement of the correct octant is still an open issue, commonly referred to as the octant problem. For example, although the results from the global fits in 2018 (see table 2.1) show a preference for the right octant (๐œƒ 23 > 45โ—ฆ ), the 3๐œŽ range for these results, as reported in [19, 20], covers both octants and also the maximum mixing case ๐œƒ 23 = 45โ—ฆ . 14 Path length = 12760.0 km 1.4 3 flavor (vacuum) 2 flavor (vacuum) 1.2 Analysis range 1.0 0.8 P 0.6 0.4 0.2 0.01 2 4 6 8 10 20 E (GeV) Figure 2.3: Showing the two flavor approximation in equation 2.17 agreement with the full three flavor oscillation probabilities in vacuum (for the oscillation parameters listed in table 2.1). The approximation depends on the condition ๐ฟ|ฮ”๐‘š 21 2 |/๐ธ โ‰ช 1 and the plot considers the largest possible path travelled by an atmospheric neutrino before detection at an experiment (i.e. IceCube). The vertical line shows the lowest neutrino energy considered in the NMO analysis performed in this work. The approximation holds quite well in the analysis energy range. Similar to the ๐œˆ ๐œ‡ disappearance case, it is possible to show that the oscillation probabilities for 2 |/๐ธ โ‰ช 1 can be approximated to: ๐œˆ๐‘’ when ๐ฟ|ฮ”๐‘š 21 2 ๐ฟ ฮ”๐‘š 31 1 โˆ’ sin2 (2๐œƒ ) sin2  ๐‘ƒ(๐œˆ๐‘’ โ†’ ๐œˆ๐‘’ ) โ‰ˆ 13 . (2.19) 4๐ธ These probabilities are therefore lower than those found for ๐œˆ ๐œ‡ disappearance (since their amplitudes are given by ๐œƒ 13 instead of ๐œƒ 23 ). As will be shown in subsequent sections, atmospheric neutrinos contain fluxes of both ๐œˆ๐‘’ and ๐œˆ ๐œ‡ and analyses will therefore include both the leading effects from ๐œˆ ๐œ‡ disappearance and the sub-leading effects from ๐œˆ๐‘’ disappearance. 15 ฮฝe,ยต,ฯ„ ฮฝe,ยต,ฯ„ ฮฝe e ฮฝฬ„e ฮฝฬ„e W Z W e, p, n e, p, n e ฮฝe e e Figure 2.4: Feynman diagrams for neutrino forward scattering in matter. All neutrino (anti- neutrino) flavors will experience a weak NC interaction potential ๐‘‰๐‘๐ถ in matter, with contributions from electrons, neutrons and protons (left). Notice this potential is independent of flavor. On the other hand, only ๐œˆ๐‘’ and ๐œˆยฏ ๐‘’ will have a weak CC interaction potential with the electrons in matter (center and right). 2.4 Neutrino oscillations in matter The oscillation probabilities in vacuum are independent on the sign of the ฮ”๐‘š 31 2 mass splitting. Sensitivity to the NMO arises, however, from matter effects in the oscillation probabilities. The discussion presented here is based on descriptions found in [43], [11], [17] and [18]. 2.4.1 The Mikheyevโ€“Smirnovโ€“Wolfenstein effect Propagation through matter involves weak interactions of neutrinos with the electrons and nucleons in the medium. For the purpose of neutrino oscillations, only coherent forward scattering (where the direction of propagation remains unchanged after the interaction) will be relevant. The Feynman diagrams for these interactions are shown in figure 2.4. Note there will be an asymmetry in the CC interaction potential for neutrinos, since matter is made of electrons (no muons or taus) and as result, the interaction potential for ๐œˆ๐‘’ will be different than that for ๐œˆ ๐œ‡ and ๐œˆ๐œ . The perturbed Hamiltonian in matter can then be written, in its flavor representation, as: ๐ป = ๐‘ˆ๐ป๐‘š 0 ๐‘ˆโ€  + ๐‘‰ , (2.20) ๐‘“ where ๐ป๐‘š 0 = ๐ธ + diag(๐‘š 2 , ๐‘š 2 , ๐‘š 2 )/2๐ธ is the mass representation of the unperturbed (free) Hamil- 1 2 3 tonian, U is the PMNS (mixing) matrix and ๐‘‰ ๐‘“ = diag(๐‘‰๐‘๐‘ , 0, 0) + ๐‘‰๐‘๐ถ , (2.21) 16 the interaction potential in the flavor representation. The first term is the interaction potential contribution for CC interactions (which is only relevant for ๐œˆ๐‘’ ), and the second term corresponds to NC interactions (and will not contribute to the oscillation probabilities since its a constant term). The CC interaction potential ๐‘‰๐ถ๐ถ is different for neutrinos and antineutrinos, and is given by: โˆš ๐‘‰๐ถ๐ถ = ยฑ 2๐บ ๐‘“ ๐‘ ๐‘’ , (2.22) where ๐‘ ๐‘’ is the electron density in the medium, ๐บ ๐‘“ (โ‰ˆ 1.16 ร— 10โˆ’5 GeVโˆ’2 ) is Fermiโ€™s constant and the sign is positive (+) for neutrinos and negative (-) for anti-neutrinos [17]. The eigenstates of this new Hamiltonian will no longer be the mass eigenstates |๐œˆ1 โŸฉ , |๐œˆ2 โŸฉ , |๐œˆ3 โŸฉ. They can, however, be found in the flavor representation by diagonalizing the Hamiltonian in equation 2.20. For simplicity, we can restrict this discussion to the two flavor case, which, as shown in section 2.3, gives a good approximation to the leading contributions of the three flavor oscillation probabilities. Using equations 2.20 and 2.21, as well as equation 2.5 for the mixing matrix U in the two flavor case, one obtains the following expression for the flavor representation of the perturbed Hamiltonian: 2๐ธ๐‘‰๐‘๐‘ ฮ”๐‘š 2 ยฉยญ โˆ’cos2๐œƒ + ฮ”๐‘š 2 sin2๐œƒ ยฎ + ๐ผ2๐‘ฅ2 (๐‘‰๐‘๐‘ /2 + ๐ธ 1 + ฮ”๐‘š 2 /4๐ธ) , ยช ๐ป= (2.23) 4๐ธ ยญ sin2๐œƒ cos2๐œƒ + 2๐ธ๐‘‰๐‘๐‘ ยฎ ยซ ฮ”๐‘š 2 ยฌ where the NC potential has already been dropped. The constant terms do not contribute to the oscillation probabilities. Diagonalization of this Hamiltonian results in two new energy eigenstates |๐œˆ1๐‘š โŸฉ and |๐œˆ2๐‘š โŸฉ, with associated (energy) eigenvalues ๐ธ 1๐‘š and ๐ธ 2๐‘š . These new eigenstates will be related to the neutrino flavor states through an effective mixing matrix: ยฉ cos๐œƒ ๐‘š sin๐œƒ ๐‘š ยช ๐‘ˆ = ยญยญ ยฎ, ยฎ (2.24) โˆ’sin๐œƒ ๐‘š cos๐œƒ ๐‘š ยซ ยฌ where the mixing angle ๐œƒ ๐‘š is defined by the relation: sin2๐œƒ sin2๐œƒ ๐‘š = โˆš๏ธ‚ . (2.25)  2 2๐ธ๐‘‰ sin2 2๐œƒ + cos2๐œƒ โˆ’ ๐‘๐‘ ฮ”๐‘š 2 17 The oscillation probabilities will therefore be the same as in equation 2.6, but the oscillation frequency (originally defined by the difference between the vacuum neutrino energies ฮ”๐‘š 2 /2๐ธ, obtained from equation 2.3 and using the approximation | ๐‘| ยฎ โ‰ˆ ๐ธ) is now defined by the difference between the new energy eigenvalues: โˆš๏ธ„ ฮ”๐‘š 2  2๐ธ๐‘‰๐‘๐‘  2 ๐ธ 2๐‘š โˆ’ ๐ธ 1๐‘š = sin2 2๐œƒ + cos2๐œƒ โˆ’ . (2.26) 2๐ธ ฮ”๐‘š 2 In addition, the vacuum mixing angle is replaced by the effective mixing angle ๐œƒ ๐‘š . Therefore: โˆš๏ธ„ ! ฮ”๐‘š 2 ๐ฟ  2๐ธ๐‘‰ 2 ๐‘ƒ๐œˆ๐‘’ โ†’๐œˆ๐‘ฅ (๐ฟ) = sin2 2๐œƒ ๐‘š sin2 ๐‘๐‘ sin2 2๐œƒ + cos2๐œƒ โˆ’ . (2.27) 4๐ธ ฮ”๐‘š 2 In matter, the amplitude of the oscillations is now given by taking the square of equation 2.25. This new oscillation amplitude, however, has a dependence on the sign of the squared mass difference ฮ”๐‘š 2 and, as a result, oscillation experiments are sensitive to the mass ordering. The oscillation probabilities shown in equation 2.27 have the following important properties: 2๐ธ๐‘‰๐ถ๐ถ โ€ข They recover the vacuum oscillation probabilities if ๐‘‰๐ถ๐ถ = 0, or for โ‰ช cos 2๐œƒ. ฮ”๐‘š 2 โ€ข As pointed out in [11], in a extremely dense medium where ๐‘‰๐ถ๐ถ โ†’ โˆž the oscillation amplitudes vanish (the amplitude in equation 2.25 tends to zero). โ€ข The matter oscillation amplitude is maximal (=1) when cos2๐œƒ โˆ’ 2๐ธ๐‘‰๐‘๐‘ /ฮ”๐‘š 2 = 0. This is called the Mikheyevโ€“Smirnovโ€“Wolfenstein (MSW) resonance. โ€ข Since the sign of ๐‘‰๐‘๐‘ changes for neutrinos and antineutrinos, differences in the oscillation probabilities for the two arise in matter, even in the absence of CP violation. โ€ข The NMO will define the sign of ฮ”๐‘š 2 . It turns out the MSW resonance condition will be satisfied only for neutrinos or antineutrinos depending on the NMO and the sign of cos2๐œƒ. As pointed in [43], this ultimately allows the determination of the sign of ฮ”๐‘š 21 2 and might eventually enable a determination of the NMO. 18 As discussed in [17], once the three flavor matter propagation Hamiltonian is diagonalized numerically, it is found that MSW resonances appear in the 1-2 sector and the 1-3 sector, but not in the 3-2 sector. As a result, the MSW resonance condition relevant to ๐œˆ ๐œ‡ and ๐œˆ๐‘’ disappearance is: 2๐ธ๐‘‰๐‘๐‘ cos2๐œƒ 13 = . (2.28) ฮ”๐‘š 31 2 This condition is satisfied either for neutrinos in the normal ordering or by antineutrinos in the inverted ordering. The importance of an NMO analysis with neutrino/antineutrino discrimination is then clearly apparent, since deviations from vacuum oscillations for either particle type are a clear indicator of the mass ordering. This will be discussed in more depth in section 2.5. 2.4.2 Parametric resonance The MSW effect is introduced in the context of a constant matter profile. However, an additional resonance effect, called parametric resonance, arises once a varying matter profile is introduced. In that case the Hamiltonian is modified via a non constant CC interaction potential that changes with time (path length), taking the form (in the two flavor case, only considering non constant terms): 2๐ธ๐‘‰๐‘๐‘ (๐‘ก) ฮ”๐‘š 2 ยฉยญ โˆ’cos2๐œƒ + ฮ”๐‘š 2 sin2๐œƒ ยช ๐ป (๐‘ก) = ยฎ. (2.29) 4๐ธ ยญ sin2๐œƒ cos2๐œƒ + 2๐ธ๐‘‰๐‘๐‘ (๐‘ก) ยฎ ยซ ฮ”๐‘š 2 ยฌ Although the evolution of the neutrino states is now found by solving a time dependent Schrรถdinger equation, for which there are no general analytic solutions, approximate solutions for systems with periodic matter profiles have been studied (e.g.[44]), and is specifically in those systems that the resonance effects have been found. Contrary to the MSW effect, the resonance arises from modifications to the oscillation phases, rather than from enhanced oscillation ampli- tudes, but a more detailed review of the mechanism is beyond the scope of this work. Parametric resonance effects are found in the oscillation probabilities of neutrinos crossing the Earthโ€™s core [44], since propagation through the Earthโ€™s matter profile can be modelled as the passage through a (single period) periodic step function potential that increases at the core and then decreases in the mantle. 19 2.5 Atmospheric neutrinos As mentioned in previous sections, atmospheric neutrinos constitute the neutrino source used in oscillation analyses with IceCube-DeepCore. As a result, an overall picture on the produc- tion mechanism of these neutrinos, their oscillation probabilities and their expected fluxes at an experiment site (i.e. IceCube) is necessary. 2.5.1 Production mechanisms Atmospheric neutrinos are produced in hadronic showers resulting from cosmic ray (highly energetic atomic nuclei of extraterrestrial origin) interactions with particles in the atmosphere. Particles generated in the showers, including charged pions and kaons, further interact with air particles, or proceed to decay, producing atmospheric neutrinos [45]. For charged pions, the dominant decay mode is to muons: ๐œ‹ ยฑ โ†’ ๐œ‡ยฑ + ๐œˆ ๐œ‡ ( ๐œˆยฏ ๐œ‡ ), (2.30) with a branching ratio of โˆผ 99.99% [16], and thus, virtually all charged pions will decay through this channel. For kaons, on the other hand, there are several relevant decay channels into neutrinos. The main channel is: Kยฑ โ†’ ๐œ‡ยฑ + ๐œˆ ๐œ‡ ( ๐œˆยฏ ๐œ‡ ), (2.31) with a branching ratio of 63.56% [16], while the secondary channels are Kยฑ โ†’ ๐œ‹ 0 + ๐‘’ ยฑ + ๐œˆ๐‘’ ( ๐œˆยฏ ๐‘’ ), (2.32) with a branching ratio of 5.07% [16] and Kยฑ โ†’ ๐œ‹ 0 + ๐œ‡ยฑ + ๐œˆ ๐œ‡ ( ๐œˆยฏ ๐œ‡ ), (2.33) with a branching ratio of 3.352% [16]. In addition, the atmospheric muons produced via pion and Kaon decay will further decay into neutrinos through the channel: ๐œ‡ยฑ โ†’ ๐‘’ ยฑ + ๐œˆยฏ ๐œ‡ (๐œˆ ๐œ‡ ) + ๐œˆ๐‘’ ( ๐œˆยฏ๐‘’ ). (2.34) 20 The atmospheric neutrino flux is dominated by its ๐œˆ ๐œ‡ component, with a estimated ratio ๐œˆ ๐œ‡ : ๐œˆ๐‘’ โ‰ˆ 2 : 1. However, the ๐œˆ๐‘’ contribution decreases further with energy since as energy increases muon lifetimes also increase via time dilation, resulting on a larger fraction of muons reaching the Earth surface and losing substantial energy before decaying [45]. 2.5.2 Atmospheric neutrino flux Modeling of the atmospheric neutrino flux requires several inputs, including: the energy spectrum of the cosmic ray primaries, the properties of the Earthโ€™s atmosphere, the geomagnetic field and modeling of the interactions and decays of the shower secondaries as they propagate through the atmosphere [45]. The cosmic ray energy spectrum has several identified features, as shown in figure 2.5. The flux follows a power law ๐‘‘ฮฆ โˆ’๐›พ with index ๐›พ โ‰ˆ 2.6 up to โˆผ 3 PeV energies. Up to these energies, ๐‘‘๐ธ โˆ ๐ธ cosmic rays are expected to originate primarily from galactic sources [46]. The region between PeV and EeV energies corresponds to the transition from galactic to extra-galactic sources [46]. At this point the spectrum becomes steeper and closely follows a second power law until energies โˆผ 7 EeV. At higher energies the spectrum again reflects a ๐›พ โ‰ˆ 2.6 spectrum due to the galactic component falling off until the spectrum is dominated by the extragalactic component [16]. A GZK cutoff at โˆผ 5 ร— 1010 GeV arises from extremely high energy cosmic rays interacting with the cosmic microwave background (CMB) [16]. The power law nature of the primary spectrum will result in an atmospheric neutrino flux that also follows a power law. As explained in [45], only the cosmic rays that interact before being deflected back into space by the geomagnetic field contribute to the atmospheric neutrino flux. This field effectively filters lower energy cosmic rays and makes the primary spectrum anisotropic [45]. The yield of atmospheric neutrinos at a particular energy can be found solving a set of coupled differential equations, know as cascade equations (see reference [48]). The system of equations is typically solved numerically using dedicated software and includes a large set of interaction coefficients that are calculated using a specific hadronic model. The total yield of neutrinos will be related to the production rates of pions and kaons and how often these particles will decay instead 21 Figure 2.5: The cosmic ray spectrum as measured by several experiments [47]. The cosmic ray flux approximately follows a power-law spectrum, with steepening at the regions known as the knee and second-knee (believed to be a consequence of the transition between cosmic rays of galactic and extragalactic origin), as well as the ankle, where the extragalactic component is believed to dominate [16]. of interacting in the atmosphere, which in turn will depend on the hadronic model used when modeling the shower evolution. The strategy used when modeling non-oscillated neutrino fluxes in the NMO analysis presented in this work is discussed in section 5.1.1. It is however helpful to present here the atmospheric neutrino flux that was ultimately used as the nominal flux for the NMO analysis, the HKKM model [49]. The expected unoscillated neutrino flux at the South Pole under this model is shown in figure 2.6. Similar to the cosmic ray spectrum, the atmospheric neutrino energy spectrum falls steeply with energy, although the spectral index is closer to ๐›พ โ‰ˆ 3.7. As discussed in [50], the power law spectrum arises from the cosmic ray spectrum with an additional ๐ธ โˆ’1 factor arising from the fact that decay lengths have a direct dependence on energy (๐œ† ๐‘‘๐‘’๐‘ โ‰ˆ ๐‘๐œ๐ธ/๐‘š; where m and ๐œ are the mass and lifetime of the decaying particle respectively). The ๐œˆ ๐œ‡ component becomes more and more dominant as energy increases and the neutrino flux is higher than the antineutrino flux since 22 Figure 2.6: The expected atmospheric neutrino flux at the South Pole, according to the HKKM model [49]. The left plot shows the energy spectrum for atmospheric neutrinos averaged over all directions and over one year (given large seasonal variations in the atmospheric density profile in the polar region). The right image shows the dependence of the flux on the zenith angle, where cos(๐œƒ ๐‘ง ) = 1 referrers to neutrinos produced in the southern sky, directly above the detection point at the South Pole. the cosmic ray primaries (atomic nuclei) are positively charged and the nuclei in the atmosphere are predominantly positively charged as well, resulting in positive mesons (pions and kaons) being produced more often [50]. Finally, the zenith dependence of the fluxes presented in figure 2.6 arises from the atmospheric density experienced by the pions and kaons as they propagate. For zenith angles close to the horizon cos(๐œƒ ๐‘ง ) โ‰ˆ 0, the particles travel through thin layers of the atmosphere for a longer time, and thus, have a higher chance to decay instead of interacting [43]. The effect is amplified for ๐œˆ๐‘’ since they require two subsequent decays [43]. 2.5.3 Atmospheric neutrino oscillations Atmospheric neutrinos provide an abundant source for oscillation studies since neutrinos pro- duced over a broad range of energies in the atmosphere will travel several different path lengths across the Earth to reach a detector. Several combinations of energy and path length can there- 23 fore be leveraged to probe oscillation probabilities. In addition, atmospheric neutrinos crossing the Earth in their path towards a detector will experience matter effects in their oscillation probabilities, and as result, oscillation experiments using atmospheric neutrinos provide sensitivity to the mass ordering. To discuss the oscillation probabilities expected for atmospheric neutrinos, the density profile of the Earth must first be modeled. In the NMO analysis performed in this work, a 12-layer simplification of the Preliminary Reference Earth Model (PREM) [51] is used (see figure 2.7). Neutrinos produced in the atmosphere are typically generated at heights between 10 and 20 km [45] and the matter density profile "seen" by a neutrino in its path from production to e.g. IceCube detector, depends on its production point within the atmosphere. The zenith angle (๐œƒ ๐‘ง ) defining the arrival direction of an atmospheric neutrino (shown in figure 2.8) can be used to infer the path length traversed as: โˆš๏ธƒ ๐ฟ(๐œƒ ๐‘ง ) = โˆ’(๐‘… โˆ’ ๐‘‘) cos(๐œƒ ๐‘ง ) + (๐‘… โˆ’ ๐‘‘) 2 cos2 (๐œƒ ๐‘ง ) + (๐‘… + โ„Ž) 2 โˆ’ (๐‘… โˆ’ ๐‘‘) 2 , (2.35) where d is the detector depth (see section 3.2.1), R is the radius of the Earth and h is the altitude of the production point in the atmosphere. The three flavor survival probabilities for muon neutrinos and antineutrinos produced in the atmosphere (see section 2.5) after traveling through the Earthโ€™s density profile, are shown in figure 2.9, for both NO and IO. The first effect that can be noticed is the oscillation probabilities for neutrinos in the IO and antineutrinos in the NO are largely consistent with the oscillation probabilities expected for vacuum oscillations, where the probabilities are constant for a given value of the ratio L/E. This is largely due to the characteristics of the MSW resonance condition in equation 2.28, that can only be fulfilled by neutrinos in the NO or antineutrinos in the IO. In those cases, the matter effects can be observed clearly in the oscillation probabilities. To understand this better, is appropriate to explore equation 2.15. Notice that beyond the multiplicative factor sin2 (๐œƒ 23 ), this expression is just the 2 flavor oscillation probabilities in equation 2.6 when considering mass states ๐œˆ1 and ๐œˆ3 . Once matter effects are introduced, in constant matter density, the 2-flavor oscillation probabilities are modified as described in section 2.4.1. Thus, 24 Figure 2.7: The Preliminary Reference Earth model (PREM) compared to the internal structure of the Earth, broadly defined by 4 geological layers [52]. The 4-layer and 12-layer approximations to the PREM model are also depicted. The density profile resembles a step profile, and therefore allows for both the MSW effect and parametric resonance effects in the oscillation probabilities. approximately: 2 ๐ฟ โˆš๏ธ„ ! ฮ”๐‘š 31  2๐ธ๐‘‰๐‘๐‘  2 ๐‘ƒ๐œˆ๐‘š๐œ‡ โ†’๐œˆ๐‘’ (๐ฟ) โ‰ˆ sin2 (๐œƒ 23 ) sin2 2๐œƒ 13 ๐‘š sin2 sin2 2๐œƒ 13 + cos2๐œƒ 13 โˆ’ (2.36) ๐ธ 2 ฮ”๐‘š 31 where: ๐‘š = sin2๐œƒ 13 sin2๐œƒ 13 โˆš๏ธ‚  2 . (2.37) sin2 2๐œƒ 13 + cos2๐œƒ 13 โˆ’ 2๐ธ๐‘‰๐‘๐‘ /ฮ”๐‘š 312 Although the oscillation probabilities to ๐œˆ๐‘’ are clearly suppressed by the smallness of ๐œƒ 13 in vacuum, once the resonance condition in equation 2.28 is satisfied, these oscillations become maximal. In turn, the MSW resonances will result in a drop in the survival probabilities for ๐œˆ ๐œ‡ . Since the Earth profile has several different density layers this effect will be seen for several different MSW resonances and the clear discontinuity observed in the oscillation probabilities at 25 Figure 2.8: A schematic illustration of the different path lengths traversed by atmospheric neutrinos from their generation in the atmosphere towards the IceCube detector. The zenith direction shown in the diagram can be converted to a path length across the Earth using equation 2.35. ๐ฟ โ‰ˆ 11000๐‘˜๐‘š is due to the large change in density between the Earthโ€™s outer core and the mantle. Parametric resonance effects are responsible for the distortions in the survival probabilities for neutrinos crossing the core (L > 11000 km) [43]. The matter effects provide a clear signature for the NMO, including that deviations from vacuum oscillations for neutrinos would indicate the NO is correct. In general, experiments that can separate neutrinos and antineutrinos may most directly leverage this effect. Sensitivity to the NMO can also be achieved using a sample of both neutrinos and antineutrinos for experiments where the two particle types can not be separated, as long as one type dominates the signal. The overall sensitivity, however, is suppressed. Finally, it is helpful to evaluate the effect of ๐œƒ 23 and ๐›ฟ๐ถ๐‘ƒ on oscillation analyses. For the mixing angle, the octant degeneracy from equation 2.18 will play a role when fitting oscillation parameters 26 Figure 2.9: Oscillation probabilities for atmospheric muon neutrinos after propagating a path length L through the earth for neutrinos and antineutrinos. Oscillation probabilities are only computed for neutrinos with zenith angles satisfying cos(๐œƒ ๐‘ง ) < 0.1. The plots show the oscillations probabilities under both the NO (top) and the IO (bottom) cases. Matter effects are visible for neutrinos in the NO and antineutrinos in the IO. since the probabilities under each octant are comparable, resulting in two local minima. Differences between the probabilities for the octants arise from sub-leading effects (given the smallness of ๐œƒ 13 ). However, for energies close to the MSW resonances, oscillations channels usually suppressed by the small value of ๐œƒ 13 are now apparent, enhancing the differences between the expected oscillation patterns obtained in each octant. This can be seen in figure 2.10, where the difference in the oscillation probabilities in each octant are shown both for neutrinos and antineutrinos in the NO. On the other hand, the effect of ๐›ฟ๐ถ๐‘ƒ in the oscillation probabilities is not pronounced at energies above 6 GeV, as shown in figure 2.11. Since the NMO analysis performed in this work studies oscillations at energies above this range, the effect of this parameter on the analysis is minimal. For 27 (a) (b) Figure 2.10: Differences in the oscillation probabilities between the left and right ๐œƒ 23 octants, for neutrinos (a) and antineutrinos (b) in the NO. The degeneracy is more pronounced for antineutrinos, where oscillation probabilities are largely consistent with vacuum oscillations. (a) (b) Figure 2.11: The differences in the oscillation probabilities between no CP violation and maximum CP violation, for neutrinos (a) and antineutrinos (b) in the NO. The effect is not visible above โˆผ 6 GeV. simplicity, only the ๐›ฟ๐ถ๐‘ƒ = 0 (no CP violation) scenario was ultimately studied. 2.6 Neutrino interactions As described in section 2.1, neutrinos can only interact via weak interactions. These interactions have the lowest cross sections across all fundamental interactions within the SM. As a result, 28 neutrino interactions are rare and tracking neutrinos directly to measure their properties is not feasible. Neutrino detectors usually detect secondary interaction particles, and use these to infer the properties of the primary neutrino. At the energy range relevant for neutrino oscillation analyses in IceCube-DeepCore (a few GeV to approximately 100 GeV), a number of interactions apply: โ€ข Elastic and quasi-elastic scattering: In this case the neutrino scatters off the entire nucleon since it doesnโ€™t have enough energy to release individual quarks by overcoming the binding force. This, in turn, usually frees nuclei from within the target atom. Charged current inter- actions of this type are referred to as quasi-elastic (QE), whereas neutral current interactions are referred to as elastic [53]. There are two main QE interactions: ๐œˆ๐‘™ + ๐‘› โ†’ ๐‘ + ๐‘™ โˆ’ and ๐œˆยฏ๐‘™ + ๐‘ โ†’ ๐‘› + ๐‘™ + where l is a lepton flavor. For elastic scattering, the relevant processes are: ๐œˆ๐‘™ ( ๐œˆยฏ๐‘™ ) + ๐‘›, ๐‘ โ†’ ๐œˆ๐‘™ ( ๐œˆยฏ๐‘™ ) + ๐‘›, ๐‘. โ€ข Resonance production: Neutrino interactions where the target nucleon is excited to a baryonic resonance [53]. The resonance state is short lived, and decays into mesons and nucleons [53]. โ€ข Deep inelastic scattering (DIS): This is the dominant interaction channel above a few GeV. In this case, neutrinos have enough energy to overcome the binding energy of the nucleon and interact with individual quarks [53]. The recoiled quark then undergoes hadronization producing a hadronic shower. The Feynman diagrams for these interactions are shown in figure 2.12. It is noted that interactions with atomic electrons are highly suppressed when compared to nucleon-neutrino interactions [53]. The only exception for this rule is the so called Glashow resonance: the resonant production of a W boson via the following process ๐œˆยฏ ๐‘’ + ๐‘’ โ†’ ๐‘Š at 6.3 PeV [53]. However, these energies are beyond the scale for oscillation analyses, so effectively, only neutrino-nucleon interactions will contribute. 29 ฮฝe,ยต,ฯ„ e, ยต, ฯ„ ฮฝe,ยต,ฯ„ ฮฝe,ยต,ฯ„ W Z Hadronic Shower Hadronic Shower N N Figure 2.12: Feynman diagrams for CC (left) and NC (right) neutrino-nucleon DIS. A diagram showing the cross section (averaged over results for neutrons and protons) for CC neutrino-nucleon interactions is shown in figure 2.13 and DIS dominates above โˆผ 10 GeV and is the only relevant channel above 100 GeV. In addition, in this region the interaction cross section for DIS scales linearly with energy. In fact, the theoretical differential cross sections for CC neutrino (and antineutrino) DIS off a nucleon are [54]: ๐‘‘๐œŽ ๐œˆ๐‘ = ๐œŽ0 [๐‘„ + (1 โˆ’ ๐‘ฆ) 2 ๐‘„], ยฏ (2.38) ๐‘‘๐‘ฆ ยฏ ๐‘‘๐œŽ ๐œˆ๐‘ = ๐œŽ0 [๐‘„ยฏ + (1 โˆ’ ๐‘ฆ) 2 ๐‘„], (2.39) ๐‘‘๐‘ฆ where ๐œŽ0 โ‰ˆ 1.5 ร— 10โˆ’42 (๐ธ/๐บ๐‘’๐‘‰), (2.40) ยฏ is the fraction of the nucleon momentum carried by quarks (anti-quarks), and ๐‘„ (๐‘„) ๐ธ ๐œˆ๐›ผ โˆ’ ๐ธ ๐‘™ ๐›ผ ๐ธ ๐‘๐‘Ž๐‘ ๐‘๐‘Ž๐‘‘๐‘’ ๐‘ฆ= = (2.41) ๐ธ ๐œˆ๐›ผ ๐ธ ๐œˆ๐›ผ is known as the inelasticity, or the fraction of the incoming neutrino energy that is transferred to the hadronic cascade (the rest of the energy is transferred to the scattered lepton, ๐‘™ ๐›ผ ). Equation 2.38 has a suppression term (๐‘ฆ โˆ’ 1) accompanying the ๐‘„ยฏ factor. This arises from ๐œˆ ๐‘žยฏ interactions being forbidden for ๐‘ฆ = 1, which can be understood from helicity considerations [54], described below. For the purpose of cross section calculations for DIS, the leptons, neutrinos and quarks 30 Figure 2.13: Cross sections for CC neutrino-nucleon interactions in the energy range relevant for neutrino oscillation analyses with IceCube-DeepCore [53]. The top diagram is for neutrinos and the bottom one for antineutrinos. At high-energies DIS dominates and cross sections increase linearly with energy. The neutrino cross section is higher, since quark-antineutrino (and neutrino-antiquark) interactions are partially suppressed. 31 ฮฝฮฑ qฬ„ lฮฑ qฬ„ Figure 2.14: Helicities for the initial (upper diagram) and final (lower diagram) states of the ๐œˆ๐›ผ + ๐‘žยฏ โˆ’1/3 โ†’ ๐‘™ ๐›ผ + ๐‘žยฏ +2/3 scattering in the center of mass frame, for the case ๐‘ฆ = 1. For each particle, the upper arrow shows the direction of motion, while the lower arrow shows the helicity. Notice that the total helicity is reversed and therefore not conserved. As a result, the interaction is suppressed. are generally considered massless, since DIS occurs at energy ranges where ๐ธ โ‰ซ ๐‘š 2 is satisfied. Massless fermions (anti-fermions) have negative (positive) helicity. In the center of mass frame, a total transfer of momentum from a neutrino to a quark requires a violation of helicity [54], as shown in figure 2.14. As a result, the case ๐‘ฆ = 1 is forbidden. The same happens in ๐œˆ๐‘ž ยฏ interactions. However, the fractions ๐‘„ and ๐‘„ยฏ are quite different. Both neutrons and protons are made of three valence quarks. Anti-quarks, on the other hand, are only found within the nucleons as virtual particles. As a result, the cross sections for antineutrinos have a larger suppression factor than those for neutrinos. This explains why the cross sections for neutrinos (obtained by integrating equations 2.38 and 2.39 over all available inelasticity values) are larger than those for antineutrinos, as seen in figure 2.13. For the analysis presented in this work, it was found that reconstructed event inelasticities provide a means to - at a modest level - separate neutrinos and antineutrino populations. This can be understood from equations 2.38 and 2.39. Evidently, the inelasticity distributions for antineutrinos will show a faster decrease in number of events as ๐‘ฆ increases that those for neutrinos ยฏ This is clearly seen in figure 2.15. Finally, cross sections for NC interactions have (since ๐‘„ > ๐‘„). similar shape than those for CC interactions, although with lower magnitude [55]. 32 Figure 2.15: Showing the differential CC cross sections for deep inelastic scattering for both neutrinos and antineutrinos (equations 2.38 and 2.39) [50]. Event counts in a detector like IceCube will decrease faster as a function of inelasticity for antineutrinos. 2.7 Physics background for neutrino/antineutrino discrimination One of the main components of the research presented in this work was to explore the possibility to separate neutrino and antineutrino events in the IceCube detector, using muon decay signals. 2.7.1 Muon propagation in matter Muons propagating through a medium (say the Antarctic ice sheet) will lose energy primarily via electromagnetic interactions. The mean muon energy lost per unit distance traveled, also known as stopping power, for muons in ice is provided in figure 2.16. Muon energy losses at intermediate energies, between 10 MeV and 100 GeV, are due mostly to interactions between the muons and the atomic electrons, leading to ionization and atomic excitation [16]. The stopping power in this regime is well described by the Bethe-Bloch formula [16]: D 1 ๐‘‘๐ธ E ๐‘˜ ๐‘ h  2 ๐‘š ๐‘’ ๐›ฝ2 ๐›พ 2๐‘Š๐‘š๐‘Ž๐‘ฅ  i โˆ’ = ln โˆ’ 2๐›ฝ2 โˆ’ ๐›ฟ(๐›ฝ๐›พ) , (2.42) ๐œŒ ๐‘‘๐‘ฅ 2๐›ฝ2 ๐ด ๐ผ๐‘’2 ๐‘“ ๐‘“ where ๐›ฝ is the velocity of the traveling muon, ๐›พ = (1 โˆ’ ๐›ฝ2 ) 1/2 , Z/A is the (effective) atomic number to atomic mass ratio for ice, ๐œŒ is the density of ice, ๐‘˜ โ‰ˆ 0.307 MeVcm2 /mol, ๐ผeff is the 33 Mass Stoping Power for muons in ice 102 Total Ionization Bremsstrahlung Pair Production Critical (MeV cm2/g) Energy Photonuclear 101 dE 1 dx 100 Minimum ionization 10 2 10 1 100 101 102 103 104 KE(GeV) Figure 2.16: The different contributions to the total stopping power of muons in ice generated using tables provided by the pdg group [56]. Above the minimum ionization energy, the losses can be approximated using equation 2.45. effective ionization potential for ice (averaged over all atomic electrons), ๐›ฟ is the so called density correction, accounting for corrections in the extension of the electromagnetic fields in the media given by polarization, relativistic effects, among others [16], and ๐‘Š๐‘š๐‘Ž๐‘ฅ is the maximum possible energy transfer in a single interaction, given by: 2๐‘š ๐‘’ ๐›ฝ2 ๐›พ 2 ๐‘Š๐‘š๐‘Ž๐‘ฅ = . (2.43) 1 + 2๐›พ๐‘š ๐‘’ /๐‘š ๐œ‡ + (๐‘š ๐‘’ /๐‘š ๐œ‡ ) 2 The Bethe-Bloch equation has a minimum value at โˆผ 0.314 GeV in ice, referred to as the minimum ionization energy. Muons above these energies (and below 100 GeV) are usually referred to as minimum ionizing muons. Notice that for minimum ionizing muons, ๐›ฝ โ‰ˆ 1 for the muon mass ๐‘š ๐œ‡ โ‰ˆ 105.6 MeV. In this range, the only dependence of equation 2.42 on energy is through ๐›พ = ๐ธ/๐‘š ๐œ‡ . Since the density corrections are small (see [16]), the Bethe-Bloch equation in the minimum ionizing regime 34 approximately satisfies: D 1 ๐‘‘๐ธ E โˆ’ โˆ ๐‘™๐‘›(๐ธ), (2.44) ๐œŒ ๐‘‘๐‘ฅ ๐‘–๐‘œ๐‘›๐‘–๐‘ง๐‘Ž๐‘ก๐‘–๐‘œ๐‘› that is, it exhibits a logarithmic increase with energy. This behavior can be clearly seen in figure 2.16 and muon energy losses due to ionization are fairly continuous. In addition, given relation 2.44, it is customary to assume a constant energy deposition rate for these particles. At higher energies, highly stochastic processes start to dominate the energy losses for muons. These processes are: โ€ข Photonuclear interactions: high-energy electromagnetic interactions where the muon trans- fers enough energy through a virtual photon to interact directly with the nucleons in the atom [16]. โ€ข Pair production: electromagnetic interactions with nuclei in the medium where particle pairs (๐‘’ โˆ’ ๐‘’ + ) are produced. โ€ข Bremsstrahlung: electromagnetic interaction of a muon with the atomic field, where the muon is decelerated and the lost energy is emitted as a photon. The average effect of each of these processes at the energy ranges relevant for this work is shown in figure 2.16. These contributions tend to increase - to very good approximation - linearly with energy. As a result, the average energy losses of muons at high energies can be approximated by [16]: D 1 ๐‘‘๐ธ E โˆ’ โ‰ˆ ๐‘Ž + ๐‘๐ธ, (2.45) ๐œŒ ๐‘‘๐‘ฅ where a is the contribution from ionization (almost constant above the minimum ionization energy) and bE is the combined contribution of the radiative processes. In reality ๐‘ = ๐‘(๐ธ), although similar to the ionization case this factor increases fairly slowly with energy. At an energy of 1.03 TeV, average losses from radiative effects and ionization contribute equally. This energy is referred to as the muon critical energy. In practice, however, when simulating muon propagation, ionization processes are usually considered continuous whereas radiative processes 35 are simulated as discrete processes, given the stochastic nature of the energy depositions for these interactions. The average range of muons in ice can be approximated using equation 2.45. Upon integration, its found that: 1  ๐‘ฅ(๐ธ 0 ) โ‰ˆ ln 1 + ๐ธ 0 /๐ธ ๐‘ , (2.46) ๐œŒ๐‘ where ๐ธ ๐‘ is the muon critical energy and ๐ธ 0 is the initial energy of the muon. A typical value of ๐‘ โ‰ˆ 3.3 ร— 106 cm2 gโˆ’1 is utilized by IceCube (see [43, 50]) providing fairly good approximation to the range of muons in ice. The importance of relation 2.46 for a decay signal analysis is the ability to infer the energy ranges at which muons (produced in ๐œˆ ๐œ‡ CC interactions) are expected to stop within the detector, since only those will provide a โ€œdecay signalโ€. 2.7.2 Muon atomic and nuclear capture Although energy depositions in ice are similar for muons and anti muons, ๐œ‡โˆ’ may capture in an atom and form a bound system. The process starts with the muons losing energy via electromagnetic interactions as they travel through a medium. Once their energies reach the keV-scale and their velocities fall below those of the valence electrons in the medium, they interact with these electrons and very rapidly (10โˆ’13 ๐‘  in insulators), come to rest [57, 58]. After the muons have come to rest, they are captured by the atom on the outermost shell [57, 58]. Following their capture, they rapidly cascade down to the 1s state (โˆผ 10โˆ’13 ๐‘ ), releasing Auger electrons and photons as they transition between orbital states [57, 58]. A consequence of the process described above is that a ๐œ‡โˆ’ "stopping" in a medium will always be captured by an atom as defined in the software used in this analysis to simulate atomic capture, Geant4 [59]. As a result, the relevant quantities for modeling atomic capture of muons are not the absolute atomic capture rates, but instead the relative capture rates of the atoms in a medium. As explained in [57], relative capture rates between two atoms with atomic numbers ๐‘1 and ๐‘2 and probabilities of muon capture ๐‘Š (๐‘1 ) and ๐‘Š (๐‘2 ), are usually reported experimentally using the notation ๐ด(๐‘1 /๐‘2 ): ๐‘  1 ๐‘Š (๐‘1 ) ๐ด = (2.47) ๐‘2 ๐‘Š (๐‘2 ) 36 Modeling of this quantity started with the Z-law approximation by Fermi and Teller [60]. Since muons ultimately stop by interacting with the atomic electrons, the approximation inferred a capture rate proportional to the total number of electrons in an atom ( the stopping power of the atoms [57]). That is: ๐‘  1 ๐ด โˆ ๐‘1 /๐‘2 (2.48) ๐‘2 This approximation had limited success in fitting experimental results and several corrections have been applied to it over the years [61]. An empirical relation developed by Stanislaus, and discussed in [61, 62] has provided a better agreement with experimental data for oxides (๐‘2 = 8) and chlorides (๐‘2 = 17), specifically: ๐‘  1 ๐ด = 0.6๐œŒ(1 + ๐›ผ๐œŒ)(๐‘2 /๐‘1 ) 1/8 (1 + 5.53๐‘‰ 5.45 ยท 10โˆ’5 ) (2.49) ๐‘2 where V is the valence of the element with value ๐‘1 , ๐œŒ is the density of the compound (in g/cm3 ) and ๐›ผ is a model parameter equal to -0.164 for oxydes with ๐‘1 < 18, -0.222 for chlorides and zero for ๐‘1 โ‰ฅ 18. It is noted that equation 2.49 is not implemented in Geant4 and the code instead uses the Z-law approximation, where the proportionality constant in equation 2.48 is derived from the review presented in [63] and equal to 0.56 for oxides, 0.66 for halogens and 1 otherwise. For the purpose of decay signals in IceCube, however, the vast majority of the medium is Hydrogen and Oxygen. As mentioned in both [57] and [61], capture in Hydrogen compounds is substantially different to the picture above, since once bound to Hydrogen, muons immediately transfer to the higher Z atom, and thus, water (ice) is effectively an Oxygen target. This effect is not accounted for in Geant4, but introducing transfers between Hydrogen and Oxygen only corrects the amount of discrimination between ๐œ‡+ and ๐œ‡โˆ’ events, not the strength of decay light signals. Re-weighting of simulated events to correct the atomic capture ratios was not performed, but it can explored in a future iteration of the NMO analysis. Once bound on the 1s state, the captured muons undergo either bound decay or nuclear capture, and the competition between the two results in a shortened "effective" lifetime. Indeed, the 37 Z Element Lifetime Capture rate ร—10โˆ’3 (๐‘ โˆ’1 ) Huff factor (Q) ๐œ‡+ 2194.90 0.450 1. 1 1H 2194.90 0.450 1. 1 2H 2194.53 0.470 1. 8 12 C 2028 37.9 1 8 13 C 2037 35.0 1 8 14 N 1919 66 1 8 16 O 1796 102.5 0.998 8 18 O 1844 88.0 0.998 Table 2.2: Nuclear capture rates and effective lifetimes for bound muons on different isotopes (based on information from [61]). Also shown is the lifetime of ๐œ‡+ in any medium (the "free" muon lifetime). interaction rate for bound muons can be written as [61]: ฮ›๐‘ก = ฮ›๐‘ + ๐‘„๐œ† ๐‘‘ , (2.50) where ฮ›๐‘ก is the total rate of interaction rate for bound ๐œ‡โˆ’ , ฮ›๐‘ is the nuclear capture rate, Q is the so-called Huff factor that accounts for differences in the intrinsic decay times of bound muons (i.e. time dilation corrections, since they are not stationary), and ฮ›๐‘ is the free decay rate in the rest frame of the muon. Nuclear capture rates are well documented for several different nuclei and table 2.2 provides the capture rates for several elements as presented in reference [61]. Since stopping ๐œ‡+ will not capture in atoms and will simply decay as a free muon at rest, the effective lifetimes for muons and anti-muons in ice are defined roughly by ๐œ๐œ‡+ = 2.197 ๐œ‡s and ๐œ๐œ‡โˆ’ โ‰ˆ ฮ›๐‘กโˆ’1 ( 16 O) = 1.796 ๐œ‡s. There is therefore an โˆผ 400 ns difference in the lifetimes of ๐œ‡+ and ๐œ‡โˆ’ in ice. 2.7.3 Decay and nuclear capture signals Decay times can be estimated in a detector using signals generated by the decay secondaries. The Feynman diagram for muon decay can be found in figure 2.17a. Notice that the only visible secondary is the outgoing electron (positron) and these are commonly referred to as "Michel electrons" ("Michel positrons"). Their energy spectrum is presented in figure 2.18 for bound decays on several nuclei, and for free decays. The maximum allowed Michel electron energy in free decays is about half the mass of the muon (๐‘š ๐œ‡ โ‰ˆ 105.66 MeV). 38 ยตโˆ’ ฮฝยต ฮฝ ยต (ฮฝยต ) ยต+(โˆ’) e+ (ฮฝ e ) W W +(โˆ’) u d p u u n ฮฝe (eโˆ’ ) d d (a) (b) Figure 2.17: Feynman diagrams for muon decay (a) and muon nuclear capture (b), the two competing processes for bound muons. Figure 2.18: Spectrum of Michel electrons for decays in several materials and for free decays [61]. 39 In the case of ๐œ‡โˆ’ nuclear capture, energy is transferred to both the outgoing neutrino and the nucleus. Usually the capture results in an excited nucleus that decays emitting several neutrons [61]. However, since only a couple of neutrons are produced and they can only produce a signal in a Cherenkov detector like IceCube (since they are neutrally charged, see section 3.1) via neutron capture, muon capture signals are observed in only a relatively small number of stopping muon events. 40 CHAPTER 3 THE ICECUBE AND DEEPCORE DETECTORS The analysis presented in this work is a neutrino mass ordering measurement performed with atmospheric neutrino events in DeepCore, a sub-array within the IceCube neutrino detector. It is therefore imperative to review the IceCube and DeepCore detectors before describing the analysis and that is accomplished in this chapter. To start, the detection principle for these detectors is introduced in section 3.1. Afterwards, the detector geometries and main components are reviewed in section 3.2. The discussion continues by introducing the event topologies produced in the detector by different neutrino interactions (see section 3.3), whereas reconstruction of the neutrino physical quantities is reviewed independently in Appendix F. The optical properties of the detector medium are then introduced in section 3.4, to aid in the understanding of detector systematics used in the NMO analysis and finally, background events for atmospheric neutrinos are reviewed in section 3.5 (the removal of these backgrounds from experimental data to obtain the atmospheric neutrino sample used in this work is reviewed independently in appendix E). 3.1 Detection principle for Cherenkov detectors As a relativistic particle travels in a dielectic medium, light is radiated coherently in a conical wavefront [64]. This type of radiation is known as Cherenkov radiation and the conical wavefront can be interpreted as the coherent interference of multiple spherical wavefronts originating at the instantaneous position of the particle as it travels (see figure 3.1). This type of radiation is known as Cherenkov radiation and the conical wavefront can be inter- preted as the coherent interference of multiple spherical wavefronts originated at the instantaneous position of the particle as it travels, as shown in figure 3.1. The angle at which radiation is emitted relative to the travel direction of the relativistic particle can be derived from the components of the electric field parallel and perpendicular to the particleโ€™s velocity, and is found to be [16, 64]: 1 cos ๐œƒ ๐‘ = (3.1) ๐›ฝ๐‘› where ๐›ฝ is the velocity of the particle (in natural units) and ๐‘› is the phase index of refraction of the 41 Figure 3.1: A schematic showing the Cherenkov waveform for a particle traveling in a dielectric medium with a speed larger than the speed of light in that medium [65]. medium. For example, the phase index of refraction in the Antartic ice surrounding the IceCube detector (see section 3.2 ) is โˆผ 1.32 at 400 nm, resulting in a Cherenkov angle ๐œƒ ๐‘ โ‰ˆ 41โ—ฆ . There is a slight variation of the phase index of refraction with wavelength, as shown in [66]. Cherenkov photons propagate in the direction of this angle with a velocity that is given by the group index of refraction of the medium:  ๐œ† ๐‘‘๐‘›  ๐‘›๐‘” = ๐‘› ยท 1 + , (3.2) ๐‘› ๐‘‘๐œ† It is noted that the conical waveform and the light emission direction are not generally perpen- dicular, unless ๐‘› is independent of the radiation wavelength ๐œ† (๐‘› = ๐‘›๐‘” ) [67]. The number of Cherenkov photons generated per unit distance travelled, in a wavelength range [๐œ†, ๐œ† + ๐‘‘๐œ†], is described by [16]: ๐‘‘2 ๐‘ 2๐œ‹๐›ผ  1  = 1โˆ’ , (3.3) ๐‘‘๐‘ฅ๐‘‘๐œ† ๐œ†2 ๐›ฝ2 ๐‘›2 (๐œ†) 42 where ๐›ผ is the fine structure constant. The Cherenkov radiation spectrum, as described in equation 3.3, is generally peaked at the ultraviolet but has a significant contribution in the optical range. High energy particle experiments can exploit this feature using photon sensors to detect the radiation as a charged particle crosses a transparent medium and relying on the known Cherenkov angle ๐œƒ ๐‘ to infer the direction of propagation of the particles. Since the amount of light emitted is proportional to the path length of the particle in the medium, the number of detected photons also serve as a proxy of the particleโ€™s energy. 3.2 The IceCube Detector The IceCube neutrino observatory is a multi-purpose Cherenkov detector comprised of three main components: the baseline IceCube detector, optimized for high-energy astrophysical neutrino searches, the DeepCore low-energy sub-array, designed for particle physics studies, and the Ice-Top surface array designed for detection of cosmic ray showers and also used for vetoing atmospheric neutrinos and muons. A schematic illustrating the size and spacial configuration of the main components can be found in figure 3.2. 3.2.1 The IceCube Optical Array The baseline IceCube detector consists of a cubic kilometer array of digital optical modules (DOMs) embedded within the Antarctic ice, at depths between 1450 m to 2450 m. The deep ice provides an optically transparent medium that is optimal for Cherenkov detection Charged secondary particles produced in neutrino interactions within and close to the instrumented volume produce Cherenkov radiation that is detected by the optical modules. The large detection volume is necessary given the low neutrino cross section, the expectation that the astrophyical flux would follow a steep power law, and to capture sufficient information to reconstruct events of TeV-scale and above [68]. The IceCube DOMs consist of a 10-inch photomultiplier tube (PMT) encapsulated in a glass pressure sphere, together with an on-board autonomous data acquisition system (DAQ). DOMs are instrumented in the ice using so called "strings", a arrangement of 60 DOMs coupled to a single vertical copper power and communications cable. The cable provides means for neighboring 43 Figure 3.2: An overall schematic of the IceCube Neutrino Observatory [69]. It includes its three main components: the baseline detector hexagonal array within the Antarctic ice, the DeepCore sub-array, and the IceTop surface array. Also shown is the location of the IceCube laboratory at the surface. DOMs to exchange information with each other and for overall communication between the DOMs and the main IceCube laboratory (counting house), located at the surface. There are 78 strings in the baseline IceCube detector, organized in an hexagonal array with horizontal distances between strings of approximately 125 m and vertical distances between DOMs in a string of 17 m [68]. This optimized geometrical configuration allows detection and reconstruction of neutrino events with energies above 100 GeV. 3.2.2 DeepCore The DeepCore detector (see figure 3.3) is a sub-array of DOMs deployed near the center of IceCube in the clearest region of the deep ice sheet. The array lowers the threshold of the observatory 44 Figure 3.3: Schematic of the IceCube detector, focused on DeepCore [70]. The top view image showcases the 8 "DeepCore" strings that are deployed at the center of the IceCube array. These, together with the 7 innermost IceCube strings, define the DeepCore detector. An extended version of DeepCore, also shown, is occasionally used when applying event selection cuts. In addition, the image showcases a so called corridor, a region in the baseline IceCube detector through which atmospheric muons can travel largely undetected until they arrive to DeepCore. The side-view image shows the optical properties of the ice as a function of depth. DeepCore DOMs were deployed in the region with best optical properties for detection of Cherenkov light. 45 to a few GeV, opening the possibility for a wide spectra of new physics analyses including neutrino oscillation measurements and dark matter searches. Each of the 8 dedicated DeepCore strings has 60 DOMs in total, where 50 DOMs are located below depths of โˆผ2100 m and tightly packed with a vertical inter DOM distance of 7 m. The remaining 10 DOMs are located above 2000 m depth, and serve as a veto for down-going atmospheric muons. The region separating the two groups of DeepCore DOMs is commonly referred as the dust layer and has the largest measured scattering and absorption in the detector. In addition, six out of the eight strings also contain high quantum efficiency (HQE) DOMs, utilizing PMTs with a โˆผ35% higher quantum efficiency than regular DOMs [68, 71]. For the purpose of this work the DeepCore detector is defined as the complete 15 string sub- array. On average, the horizontal distance between the strings is approximately 72 m [68, 71] (about 1.5 times the effective scattering length of the deep ice within the sub-array [71]). The higher quantum efficiency, cleaner ice and compact geometry all contribute towards DeepCoreโ€™s reduced energy threshold. 3.2.3 IceTop The IceTop surface array is comprised of pairs of cylindrical tanks filled with clear ice, each containing two standard DOMs at the same vertical position and with a horizontal distance of 0.58 m between the DOM centers (each at a distance of 0.29 m to the center of the tank) [72]. The DOMs are embedded within the ice and detect cosmic ray shower particles via Cherenkov radiation. There are in total 81 IceTop stations, each station consisting of a pair of tanks separated by a 10 m distance [72]. The stations are located on a similar grid to IceCube, as shown in figure 3.4. Information from multiple stations can be used to detect an incoming shower, infer its overall size, and provide a veto to atmospheric muons. 3.2.4 The IceCube Lab (ICL) The ICL is the main operations center for the IceCube detector, located on the surface of the ice approximately above the center of the detector (see figure 3.2). The building contains the so-called South Pole System (SPS), a set of computer hubs used for data acquisition and online processing 46 Figure 3.4: Schematic of the IceTop surface array [72], showing the positions of the 81 IceTop stations relative to the baseline IceCube strings (marked as holes in the figure). The In-fill refers to a more densely instrumented area in IceTop, approximately located in the region above DeepCore. The figure also shows the location of the IceCube Lab. at the South Pole. String cables enter two towers at each side of the structure to reach the server room containing the SPS [68]. 3.2.5 Detector Coordinates This document refers to "detector coordinates" when discussing neutrino event selection cuts and reconstructions. For reference, the z-coordinate in IceCube is the vertical direction (pointing upwards towards the Southern sky), where z=0 corresponds roughly to the detector center at a depth of 1948 m (see figure 3.3). The x and y detector coordinates are defined as in figure 3.4, where the positive y direction follows the prime meridian [73] and the origin is roughly at the detector center. Radial directions in IceCube are defined as within the xy-plane. The zenith angle is defined as the angle relative to the z-axis, as shown in figure 2.8. Azimuthal directions are defined such that a 47 particle traveling from positive to negative x is 0โ—ฆ azimuth, and a particle traveling from positive to negative y is 90โ—ฆ azimuth. Notice that the zenith and azimuth definitions follow typical spherical coordinates, but used to describe the direction from which particles are arriving. 3.3 Neutrino event topologies in IceCube Cherenkov radiation produced by charged secondaries from neutrino interactions close or within the detector volume is detected by the DOMs and can be used to infer the type of interaction that occurred, as well as several physical quantities of the neutrinos. A detailed discussion of the reconstruction algorithms required for the NMO analysis is provided in Appendix F, including a detailed review of energy, vertex and directional reconstruction of neutrino events in DeepCore. This section presents a more general review of the main light deposition topologies expected from neutrino interactions within or close to the detector volume. Event by event identification of the topologies using detector data is addressed on detail in Appendix E. 3.3.1 Cascades This event topology is characterized by a quasi-spherical light emission pattern, similar to the signal that would be expected from a point source of light at the neutrino interaction vertex. They are typically the result of hadronic or EM showers produced in neutrino interactions, since for all relevant energies in IceCube and DeepCore (few GeV to several PeV), these showers have lengths of the order of 10 m or less (see figure 3.5) and are therefore effectively point-like when compared to the overall dimensions of the detector. Cascades are the expected signature for typical (DIS) NC interactions within IceCube, where the only visible secondaries are charged particles within the hadronic shower generated from the recoiled nucleus. For typical (DIS) CC interactions, on the other hand, the event topology will also include the light deposited by the outgoing leptons: โ€ข For ๐œˆ๐‘’ CC events: the outgoing electron develops an EM shower and the cascade topology remains, although as a superposition of both hadronic and EM showers. In the special case of a Glashow resonance, the W will decay into hadrons with a branching ratio of 67.41% [16], and cascades are therefore the main topology expected for these interactions. 48 Figure 3.5: The mean ranges of heavy charged leptons and cascades (particle showers) in water. Notice range only increases significantly as a function of energy for the charged leptons (the exception being electrons, which tend to produce EM showers). In general, cascade dimensions reach up to about 10 m and are therefore effectively point-like when compared to the dimensions of the IceCube detector (from [74]). โ€ข For ๐œˆ ๐œ‡ CC events: only low energy outgoing muons will have a short enough path length to generate a cascade-like signature. DeepCore analyses in the past have used an upper limit of 50 m (roughly muon energies below 10 GeV) for cascade-like events [70]. โ€ข For ๐œˆ๐œ CC interactions: the tau decay time and channel play a mayor role in the event topology. Higher mass leptons have a lower energy loss rate in matter [65]. However, taus will only be more penetrating than muons at high energies, where time dilation allows taus to travel for a long enough time before decaying (the tau lifetime is โˆผ 0.29 ps). Tau light deposition will therefore remain more cascade-like than muons at higher energies (up to โˆผ 1 PeV) as seen in figure 3.5. There are several different decay channels for taus, and two that result in a cascade topology; the preferred channel into hadrons (๐œ โ†’ ๐œˆ๐œ + hadrons; branching ratio 49 Figure 3.6: A schematic showing a cascade topology [43]. Light is emitted outwards from a nearly point-like source within the detector, producing a quasi-spherical light detection pattern. Earlier PMT pulses containing a large fraction of the total deposited charge are detected close to the interaction point. Cascade events within IceCube will generally exhibit subtle asymmetries arising from the finite dimensions and directionality of the light source. of โˆผ 64.8 %), where the final result will be a hadronic shower, and the decay channel into electrons ( ๐œ โ†’ ๐œˆ๐œ + ๐‘’ + ๐œˆยฏ ๐‘’ ; branching ratio of โˆผ 17.8 %), where the final result is an EM shower. An example of a high energy cascade event in IceCube can be found in figure 3.7. Although quasi-spherical, these events still contain directional information on the incoming neutrino in the form of asymmetries in the light deposition pattern, although they typically have worst directional resolution compared to track-like events. An advantage of cascade events is their energy resolution, since is common for these events to have significant fractions of their charged interaction secondaries deposit their energy in the detector volume. Neutrino energy estimation is best for CC events, 50 Figure 3.7: A visual representation of a track event (left), a cascade event (center) and a double- bang event (right) [46]. High energy interactions are depicted to facilitate identification of the main features of each topology. Colored spheres represent a DOM that has detected light in the event. The size of the spheres indicates the amount of charge (light) observed by the DOMs. Time is represented using a color scale, with early detection depicted in red and late detection depicted in blue. since in the case of NC events, a large fraction of the original neutrino energy is carried by an outgoing neutrino and the probability that this neutrino subsequently interacts within the detector is negligible. 3.3.2 Tracks Track signatures are elongated, linear light depositions resulting from the passage of a charged particle through the detector. Most commonly, ๐œˆ ๐œ‡ CC interactions above โˆผ10 GeV may produce an outgoing muon of sufficient energy to produce this signature as do ๐œˆ๐œ CC interactions when the tau decays into muons ( ๐œ โ†’ ๐œˆ๐œ + ๐œ‡ + ๐œˆยฏ ๐œ‡ ; branching ratio of โˆผ 17.4 %). In both cases, there will be a cascade at the neutrino interaction vertex and a track signature that starts at the interaction vertex, from the light produced by an outgoing muon. An example of a high energy track can be found in figure 3.7. At high energies the outgoing lepton can be energetic enough to leave the detector before stopping. As a result, these events have poorer energy resolution than cascades, since they are not fully contained within the detector, but the directional resolution is significantly better because their path through the detector can be traced over hundreds of meters. At low energies (below approximately 100 GeV), the differences between tracks and cascades are substantially more challenging to discern, both because of the limited number of detected photons per event at these energies and as a result of shorter muon path lengths (see figure 3.8). Since this energy range is key for performing an atmospheric neutrino oscillation analysis, a detailed 51 Figure 3.8: Showing ๐œˆ๐‘’ and ๐œˆ ๐œ‡ CC events at low energies (โˆผ 50 GeV), where identification of the event shape is substantially more challenging given the low amount of light and the short muon path lengths [43]. The analysis presented in this work applied a boosted decision tree trained to discriminate between true track (๐œˆ ๐œ‡ CC) and cascade (NC and ๐œˆ๐‘’ , ๐œˆ๐œ CC) events. description on how to perform topological identification in an event by event basis at low energies is provided in Appendix E.6. Directional and energy resolutions are also comparable at very low energies for tracks and cascades, since muon paths are contained and short. The reconstruction algorithms used for the NMO analysis are described in Appendix F, where resolutions at low energies are also presented. 3.3.3 Tau Neutrino โ€œDouble-bangsโ€ The โ€œdouble-bangโ€ signature is a particular mixture of the track and cascade topologies, found in very high energy ๐œˆ๐œ events. It consist of the original hadronic shower at the neutrino interaction vertex, an outgoing track, produced by a high energy tau, and a second cascade at the point where the tau decays. The final topology is therefore that of two cascades joined by a single track (see figure 3.7). In general, double-bang signatures will not be relevant for the (low energy) analysis described. 52 3.4 The Antarctic ice Modeling of the South Pole ice is extremely important to properly simulate the propagation of Cherenkov light within the detector as photon propagation will be affected by both scattering and absorption within the ice. The properties of the ice itself are related to how it was deposited throughout several geological epochs. Snow is deposited over thousands of years, and the snow layers at larger depths were compressed over time, until they became dense ice [50]. In this process, different concentrations of dust were also deposited, with the amount being correlated to the composition of the atmosphere throughout the epochs. As a result of this process, the optical properties and dust concentrations of the ice vary with depth. Early ice models in IceCube defined the ice as a series of stacked quasi-horizontal layers, with constant optical properties within each layer [75]. The reason why these layers were not perfectly horizontal (but instead had a slight tilt), is because the ice is deposited over a bedrock (shown in figure 3.2), which isnโ€™t perfectly horizontal [50]. However, there are also differences in the optical properties at constant depth, arising from the glacial flow of the ice [50] which has a depth dependent velocity [76]. Current ice models, including that applied in this work, include this anisotropy within each ice layer. The ice model properties can be summarized in terms of two quantities. The first one is the absorption coefficient: 1 ๐‘Ž= , (3.4) ๐œ†๐‘Ž where ๐œ† ๐‘Ž is the absorption length. The second one is the scattering coefficient: 1 ๐‘๐‘’ = , (3.5) ๐œ†๐‘’ where ๐œ† ๐‘’ is the effective scattering length or the length scale over which the direction of propagation of the scattered photon is randomized [75]. The depth dependence of the scattering and absorption for the model used in this work is shown in figure 3.9. These scattering and absorption parameters affect the expected detection rate of neutrino events in the ice, and as a result, become nuisance parameters in the physics analysis. The modeling of the ice, however, is not complete in this picture. In order to deploy IceCube strings into the ice, a hot water jet was used to melt the Antarctic ice. Strings were then lowered 53 Figure 3.9: The South Pole scattering and absorption coefficients as a function of depth for the ice model used in this work [77]. Best optical properties are obtained at the depth range defined by DeepCore. into these water-filled holes with the water refreezing to fix the strings in place. The properties of this refrozen ice are significantly different from those of the natural ice. This โ€œhole iceโ€ is usually modeled and included in IceCube simulation as a corrected angular acceptance curve for the optical modules. The angular acceptance curve of the DOMs is a function of ๐œ‚, the arrival angle of photons (where ๐œ‚ = 0 correspond to up-going photons), the shape of which depends on two independent parameters: ๐‘ 0 and ๐‘ 1 . These are ultimately systematic parameters in physics analyses. An example of the angular acceptance curves obtained for different sets of these parameters can be found in figure 3.10. 3.5 Background and noise signals IceCube oscillation analyses are performed using low energy atmospheric neutrino events detected in the DeepCore volume. The main backgrounds to these events are atmospheric muons and noise triggered events. In the following section we present a brief review of each background and how it can be mitigated. The specifics of the selection criteria for the neutrino sample used in 54 Figure 3.10: The effect of changing the two shape parameters ๐‘ 0 (left) and ๐‘ 1 (right) used to model the angular acceptance curve of DOMs when including the effects from the refrozen hole ice surrounding them [78]. the NMO analysis presented in this work is outlined in Appendix E. In addition, as discussed below afterpulse signals originating in the DOM PMTs also impact the identification of decay light signatures for discriminating neutrinos from antineutrinos. 3.5.1 Atmospheric muons Atmospheric muons represent the primary background to detecting neutrino events in the IceCube detector. Muons produced in the northern hemisphere will lose their energy before they can reach the detector and are therefore not a concern for neutrino oscillation analyses. However, high energy atmospheric muons produced in the southern hemisphere may reach the detector, losing only a fraction of their energy while traveling through the ice immediately above IceCube. Although several muons are generally produced in a single shower, most range out by the time they arrive at the detector volume. A muon bundle event occurs when several muons from a single shower arrive simultaneously to the detector. It is also possible that muons produced in two independent showers arrive to the detector within the same trigger window, referred to as a coincident event. The detection rate of atmospheric muons is several orders of magnitude higher than that of atmospheric neutrinos at the simple trigger level. Several techniques can be employed to remove this background, including vetoing them using DOMs in the IceCube array and looking for well defined starting events where the reconstructed neutrino interaction vertex is well within the detector fiducial volume. The specific selection criteria used in this analysis to remove these background 55 events is discussed in Appendix E. 3.5.2 Dark Noise Dark Noise refers to signals detected by the optical modules when no external light source is present. Four main components of dark noise are currently modeled in IceCube simulations: electronics noise (arising within the readout logic of each DOM), thermionic emission in the photo- cathode, radioactive decays within the DOM glass sphere and subsequent scintillation/fluorescence light in the DOM glass [68]. The number of photons from thermal noise and radioactive decays in an event are usually simulated using a Poisson model, with average rates for both processes being 20 Hz (with depth variations, since the rate depends on temperature) and โˆผ 250 Hz respectively, whereas their times are modelled using a uniform distribution. These first noise components are usually referred to as uncorrelated noise. On the other hand, scintillation is referred to as correlated noise, since the time of these pulses is related to the time of their precursor pulse. Separation in time between subsequent scintillation pulses is simulated using a log-normal distribution, based on an empirical model. The description above on the modeling of the noise components was based on the work presented in [79]. 3.5.3 PMT effects When a photon arrives at the PMT photocathode within a DOM, it may free an electron via the photoelectric effect. The photoelectron will be accelerated though a series of dynodes where at each step it collides with an electrode producing more electrons. The process continues until a shower of electrons is produced. The charge distribution from a single photoelectron will typically peak at the amplification factor, called the gain, which is defined by the PMT voltage. The transit time of the photoelectrons through the dynode chain will also be defined by this voltage. Uncertainty in the arrival time of the photons as a result of the spread of the arrival time distribution at a given voltage is referred to as time jitter. There are PMT effects that produce additional types of pulses, different from the typical ones described before. Early pulses, referred to as prepulses, can be produced when a photon travels 56 through the photocathode and instead releases a photoelectron on the first dynode. This electron will then have a shorter transit time. When pulses are reconstructed from digitized data produced by the DOM DAQ [68], transit times are subtracted to infer the original signal times. For prepulses, this reconstructed time will be earlier than the actual photon arrival time. Several prepulses can therefore have an effect on event reconstruction. The time scale of this effect is of the order of โˆผ 30 ns [80]. Late pulses, on the other hand, can be produced as a result of photoelectrons scattering from the first dynode back into the photocathode. The photon arrival times inferred for late pulses can be delayed up to about 200 ns, with respect to the real ones. Finally, PMTs can produce exceedingly delayed secondary pulses, called afterpulses, that will require a more detailed description, given their importance in the neutrino/antineutrino discrimi- nation method presented in this work. Afterpulses are reviewed in detail in Appendix C 57 CHAPTER 4 NEUTRINO/ANTINEUTRINO CLASSIFICATION The NMO analysis presented in this thesis was performed on a highly pure sample of track-like atmospheric neutrino events in DeepCore, with a livetime of 7.5 years [1]. The selection criteria used to realize this neutrino sample (see Appendix E for details) was developed by the IceCube collaboration specifically for oscillation analyses and was tested on a dedicated Monte Carlo (MC) sample (see Appendix D) that included: a total number of neutrino events equivalent to 70 years of effective detector livetime for each neutrino flavor; simulated background events (atmospheric muons and noise triggered events) [1]. A modified version of this MC sample was used in the NMO analysis presented in this work, where the ๐œˆ ๐œ‡ events were reprocessed to include muon capture physics since this process allows differences in the decay times for muons and antimuons, as explained in section 2.7.2. The reprocessed ๐œˆ ๐œ‡ sample contained 22.6 years of effective detector livetime. The neutrino sample selection includes reconstructions (see Appendix F) to recover neutrino energy, direction and particle identification (separation of tracks and cascades), from the light signals deposited in the detector. It is possible to probe the oscillation probabilities of track and cascade events using these reconstructed variables, since direction translate directly to atmospheric neutrino path length. However, as explained in section 2.5.3, sensitivity to the NMO arises from matter effects in the oscillation probabilities of neutrinos in the NO and antineutrinos in the IO. Pure neutrinos (antineutrino) samples allow for a simple criteria for discrimination of the mass ordering: if oscillations are consistent with matter effects then the NO (IO) is correct, otherwise, the result is the IO (NO). Since discrimination between neutrinos and antineutrinos is generally not possible in the Ice- Cube detector, previous measurements of the NMO relied on the fact that DeepCore detects more neutrinos than antineutrinos given the increased cross section of neutrinos and the differences in their atmospheric fluxes. Therefore, in the NO (where matter effects appear for neutrinos), a sample of ๐œˆ + ๐œˆยฏ has a larger contribution from matter effects than in the IO, and this could be used to obtain 58 a small but definite sensitivity to the NMO. The work presented in this chapter introduces the development of a neutrino/antineutrino classifier for ๐œˆ ๐œ‡ CC (track-like) events in DeepCore using muon decay signals. As explained in section 2.7, differences in the decay times of ๐œ‡+ and ๐œ‡โˆ’ produced in atmospheric neutrino interactions would provide a method for ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ discrimination. However, for this to be possible, decay times related to the neutrino events must be extracted, in turn requiring the identification and extraction of decay signals in the detector. While it was found in the course of this work that decay signals measured in DeepCore cannot be extracted at a sufficient level to significantly influence the obtained NMO sensitivity (see section 5.3), the study sets a promising method for future use of the technique in the IceCube Upgrade (see Chapter 6). 4.1 Muon decay signals in DeepCore Using the ๐œˆ ๐œ‡ events in the analysis sample presented in Appendix E, a comprehensive study on decay signal extraction in DeepCore was performed. Events prior to the final cuts at L7 (see Appendix E.6) were studied to increase MC statistics available for the studies, while retaining the distributions found in a pure atmospheric neutrino sample. The first and simplest step in this study was a check on the number of ๐œˆ ๐œ‡ events where a muon decay secondary (Michel electron or positron) produced light that generated a digitized pulse within a DeepCore DOM. The study indicated this occurs in about 6.1% of all events in the neutrino sample. The relatively low percentage of events with detected decay light is explained by the low typical energy of the Michel electrons (or positrons), โˆผ 52 MeV (equivalent path lengths of about 17 cm), which is several orders of magnitude below the energy threshold for DeepCore. This low energy is translated to an extremely low detection probability per DOM. For example, a simple calculation can be used to verify the likelihood to detect a pulse at a DeepCore DOM located at a distance r from the decay point, demonstrating the modest detection probability per DOM. Integrating the differential Cherenkov emission rate, given by equation 3.3, in the wavelength acceptance range for IceCube modules ([300,650] nm) returns a total emission rate of 350 photons/cm, or an average of 5950 photons per decay. The wavelength dependent quantum efficiency curve for DeepCore DOMs 59 can be incorporated in this integration to find the percentage of these photons that would be detected in an optical module, upon reaching its surface. The result is approximately 560 (416) photons using DeepCore (IceCube baseline) DOMs. Once including the DOM angular acceptance curve, the number is reduced by approximately half. Finally including the geometric effect from photon propagation. for example, assuming the photons are produced isotropically, then the fraction of photons detected by a DOM with radius ๐‘… ๐ท๐‘‚ ๐‘€ (โˆผ 16.5 cm) at a distance r from the light source (decay point) can be approximated as: ๐‘… 2๐ท๐‘‚ ๐‘€ ๐‘“ = , (4.1) 4๐‘Ÿ 2 and thus, an average number of โˆผ 1.9๐‘Ÿ โˆ’2 (r in meters) photons detected per decay is expected. This simple calculation returns about one photon detected every 50 decays, for a DOM at 10 m. Notice that although all DOMs within the detector will contribute to the detection of decay pulses, only a small number on average will be close enough to have a sizeable chance of detection. An initial conclusion of the above is: regardless of how efficient the extraction of decay signals becomes, only up to 6% of all events may be classified using this information. The next step consisted on testing whether the muon capture simulation was properly generated, and there were indeed differences in the true decay times of muons in ๐œˆ ๐œ‡ and ๐œˆยฏ ๐œ‡ events. As seen in Figure 4.1, the difference in the decay times is evident, although the separation between the two populations is potentially challenging. Identifying variables that could be used for extraction of the decay signals in ๐œˆ ๐œ‡ CC events makes uses of the following properties: โ€ข Most decay pulses will be detected by the DOMs closest to the decay point, as evident in equation 4.1. โ€ข A considerable fraction of the decay pulses can be expected at largely delayed times with respect to the main signal in ๐œˆ ๐œ‡ CC events - Cherenkov light produced by the outgoing muon and the hadronic cascade -, especially since at rest decays always happen after all muon light has been generated. 60 Figure 4.1: The truth decay time distributions (relative to the truth muon stopping time in each event) for ๐œˆ ๐œ‡ and ๐œˆยฏ ๐œ‡ events. Notice ๐œ‡โˆ’ have a shorter effective lifetime than ๐œ‡+ , and as result, earlier decays are more likely to happen in ๐œˆ ๐œ‡ events. The separation is complicated, however, since the distributions are very similar, given the small difference in the effective lifetimes for ๐œ‡ยฑ (about 400 ns). The underlying idea is to separate pulses with these characteristics from pulses generated by other mechanisms. The other types of pulses that could be expected in ๐œˆ ๐œ‡ CC events included: โ€ข Muon track pulses (outgoing muon) and cascade pulses (from the hadronic cascade at the neutrino interaction point): The main light deposition from these sources is attributed to unscattered Cherenkov photons and happens mostly within the first 430 ns after the neutrino interaction; however, a subdominant contribution from scattered light can be found at later times, where it becomes a measurable background to the decay signals. โ€ข Noise pulses: An approximate uniform distribution in time and thus can be found at the 61 delayed time window where decay pulses are expected. โ€ข Afterpulses: Produced by precursor pulses (see appendix C) and detected in a time window up to 11.3 ๐œ‡๐‘  after the original pulse. A main challenge of defining variables for the decay signal extraction was identifying the physical processes responsible for each observed pulse since once pulses are digitized within a DOM, and subsequently reconstructed (see Appendix A.2), there is no one-to-one mapping from reconstructed pulses to their light sources in the IceCube simulation chain. Individual pulses that are close in time are usually merged in the digitization and are no longer recognized as individual pulses by the reconstruction algorithms. To solve this problem, a simple tagging algorithm was designed for the purpose of decay signal extraction. It was based on time matching between reconstructed pulses and PMT pulses prior to digitization (called MCPulses in this context), and consisted of the following steps: โ€ข A brief study was performed with 10000 simulated MCPulses. Each pulse was passed through the digitization and pulse reconstruction algorithms, assuming digitization at the DOM ATWD (FADC) [68]. The detection time of the resulting reconstructed pulses was compared to that of the original pulses, and a standard deviation on the time resolution of 2.1 ns (4.2 ns) was identified. โ€ข For each event, each MCPulse was directly matched to a light source: decay, muon, cascade, noise or afterpulse. โ€ข For each DOM, the detection time of each reconstructed pulse was compared to the detection time of all MCPulses, by checking: ๏ฃฑ ๏ฃด ๏ฃฒ [โˆ’4.2, 4.2]ns for ATWD ๏ฃด ๏ฃด ๏ฃด ๐‘ก ๐‘€๐ถ๐‘ƒ๐‘ข๐‘™๐‘ ๐‘’ โˆ’ ๐‘ก๐‘Ÿ๐‘’๐‘๐‘œ โˆˆ (4.2) ๏ฃด ๏ฃด [โˆ’8.4, 8.4]ns for FADC ๏ฃด ๏ฃด ๏ฃณ 62 where the time windows were selected based on the resolutions identified in the first step above, covering the vast majority of reconstructed pulses while avoiding excessive overlap between subsequent pulses within a single DOM launch (data taking instance) [68]. โ€ข For each reconstructed pulse, a list of all MCPulses satisfying the criteria in the previous step was identified. All light sources associated to these MCPulses were then assigned to the reconstructed pulse. As a result, reconstructed pulses were separated into the following tagged categories: decay, muon, cascade, noise, afterpulse, mixed with decay light, mixed without decay light, undefined. The last category represents all reconstructed pulses that couldnโ€™t be matched to any MCPulses. Figure 4.2 illustrates the number of reconstructed pulses identified in each of the tagging categories. Although there is a large number of pulses in the unidentified category, extending the matching time window to try to include the missing pulses only marginally decreased the fraction of unidentified pulses. The missing pulses are therefore associated with the finite efficiency of the pulse reconstruction algorithm and digitizer noise. The figure also shows the effect of performing a first cut that improves the extraction of decay signals, namely, searching for signals on DOMs within 100 m from the reconstructed muon endpoint. Is evident that this removes only a marginal number of decay pulses, while substantially removing noise and other background pulses. To test whether the remaining pulses were indeed properly tagged, time distributions comparing the MCPulses and tagged reconstructed pulses under each category were generated, as shown in figure 4.3. It can be seen that there is a good match between the two, with discrepancies arising from the 6.4๐œ‡s readout window of the FADC. Once the reconstructed pulses were tagged, it was finally possible to identify per pulse variables that could be used for a binary classifier between decay pulses and other pulses: โ€ข Pulse charge: MPE afterpulses (see Appendix C) can be excluded using this variable. โ€ข The time difference between pulse detection and the neutrino interaction time: muon light should be found at earlier times than the average decay time window. 63 108 Pulses within 100 m All pulses 107 106 105 104 103 102 Afterpulse MixedDecay MixedNoDecay Michel Muon Noise Cascade NuclearCapture NotDefined Figure 4.2: Showing the number of pulses identified on each category of the reconstructed pulse tagging algorithm. To accomplish the goal of extracting decay (Michel) pulses, a first cut was introduced to study only pulses detected within a 100 m radius from the reconstructed decay point (muon endpoint). This removes a considerable number of noise pulses, while only marginally removing Michel pulses. โ€ข The delay between pulse detection time and the detection time of the earliest pulse in the receiving DOM: introduced to aid in the removal of afterpulses, since they are generated by precursor pulses. โ€ข The distance between the detected pulse DOM and the reconstructed muon endpoint: this helps enhance the decay signal, that becomes stronger in DOMs closest to the decay point. โ€ข The distance between the detected pulse DOM and the reconstructed interaction vertex: this variable helps in the removal of the Cherenkov light produced by the hadronic cascade at the interaction vertex. 64 Muon Cascade 107 Truth 106 Truth Reco Reco 105 105 Charge (PE) Charge (PE) 104 103 103 102 101 101 100 1.0 1 Truth Reco 0.5 Truth Reco 0 Truth Truth 1 0.0 10000 5000 0 5000 10000 10000 5000 0 5000 10000 Reco_vertex_rel_time (ns) Reco_vertex_rel_time (ns) Noise Afterpulse 106 Truth 105 Truth 105 Reco Reco 104 Charge (PE) Charge (PE) 104 103 103 102 102 101 101 100 1.00 1.0 Truth Reco Truth Reco Truth Truth 0.6 0.50 0.4 0.25 10000 5000 0 5000 10000 10000 5000 0 5000 10000 Reco_vertex_rel_time (ns) Reco_vertex_rel_time (ns) Decay Truth Reco 103 Charge (PE) 102 101 1.00 Truth Reco 0.50 Truth 0.25 10000 5000 0 5000 10000 Reco_vertex_rel_time (ns) Figure 4.3: Showing cumulative event distributions of all pulses over all neutrino events in the analysis sample for detection times of pulses relative to the reconstructed neutrino interaction times. Distributions for both true pulses (at the PMT anode) and reconstructed pulses (after digitization within a DOM) are compared. The light source that generated the true pulses can be retrieved directly, whereas matching reco pulses to their light sources requires a dedicated tagging algorithm, developed in this work. The distributions for pulses produced by the muon, interaction cascade, noise, afterpulses, and decays show good agreement between reconstructed and true pulses up to differences introduced by the readout window of the FADC digitizers (6.4 ๐œ‡s after a DOM launch) and the trigger window (see Appendix A.2). 65 โ€ข Pulse based muon field angle: This variable is defined in a coordinate system centered at the muon endpoint in each event. For each pulse we determine the position vector rDOM of the receiving DOM and the unitary direction vector of the muon track rฬ‚. If the angle between these vectors is smaller than the Cherenkov angle, direct light from the muon track should not reach the receiver DOM, and thus, light is more likely to have been produced by a decay secondary. The lack of track light is particularly important because it also prevents the generation of afterpulses in these DOMs. The distributions related to each variable were found to be different for all classification cate- gories, as shown in figure 4.4, and thus, promising in the extraction of decay signals. Since the best linear combination of these variables was not evident, it was estimated by implementing and training a BDT. After several attempts at implementing an optimal decay signal classifier, however, it was clear that no combination of the variables allowed for an efficient identification of decay pulses. The reason can be seen in the score distributions for the best BDT classifier found by using these training variables, as shown in figure 4.5. This was a binary classifier between decay pulses and background pulses. It was trained with 70% of the pulses in the full baseline analysis sample, and tested in the remaining 30% of the pulses. Given a set of input variables, the classifier calculates a probability score that the input features correspond to a decay pulse. Normalizing the distribution for the probability score for several background pulses and the decay pulses, it was clear that the classifier was working as intended (i.e, assigning high probabilities to true decay pulses, and low probabilities to background pulses). However, there were underlying issues identified in the classification for DeepCore events: โ€ข The extremely low detection probability for decay pulses: the detection probability of decay signals is strongly dependent on the density of optical modules in the close proximity of the decay point. For DeepCore, it is found that decay signals have a significantly lower detection rate than other pulse types, reducing the effective discrimination power provided by the BDT classifier. 66 Decay 10 1 Decay Decay 100 Muon Muon Muon 10 2 Noise Noise Noise Cascade 10 2 Cascade Cascade Afterpulse Afterpulse 10 1 Afterpulse 10 3 10 3 10 4 10 4 10 2 10 5 10 5 10 3 10 6 10 6 5000 2500 0 2500 5000 7500 10000 10 4 10000 5000 0 5000 10000 0 20 40 60 80 100 tpulse tinteractionvextex (ns) tpulse tfirst pulse inDOM (ns) DOM to endpoint distance (m) 100 Decay Decay Muon 102 Muon Noise Noise 10 1 Cascade Cascade Afterpulse Afterpulse 101 10 2 100 10 3 10 1 10 4 0 50 100 150 1.0 0.5 0.0 0.5 1.0 DOM to muon vertex dist (m) cos ( field) Figure 4.4: Normalized distributions of the per pulse variables used to train the decay pulse classifier used in this work. Pulses farther than 100 m of the reconstructed muon endpoint were not used for training. Differences in the distributions for decay pulses (blue-filled histogram) and other pulse types can be found for all variables. โ€ข The difficulty in separating decay pulses from afterpulses: the distribution of BDT scores for decay pulses peaks at maximum probability, while most of the backgrounds light sources (noise included) generate pulses with BDT scores that peak at lower values with the exception of afterpulses. These represent the main background to decay signal extraction. To understand this better it is important to look back into the main features of the decay signals and the reason afterpulses could not be rejected successfully. The main feature of decay pulses is their delayed detection time, relative to the main light depositions in an event. However, afterpulses also appear at these late time windows, at a significantly higher rate. To reject afterpulses, the strategy envisioned was to look for single pulses in DOMs, given the fact that afterpulses are generated from precursor pulses. However, from figure 4.4 looking at the ๐‘ก ๐‘๐‘ข๐‘™๐‘ ๐‘’ โˆ’๐‘ก ๐‘“ ๐‘–๐‘Ÿ ๐‘ ๐‘ก ๐‘๐‘ข๐‘™๐‘ ๐‘’ ๐‘–๐‘› ๐ท๐‘‚ ๐‘€ distribution afterpulses are detected as the original pulse in a DOM in many instances, shown by the peak at zero in the afterpulse curve. This anomaly was investigated and it was found that 67 Unweighted_Train 104 Decay Muon Pulse counts (per year) 103 Noise Cascade Afterpulse 102 MixedOnlyBkg MixedDecayBkg 101 Unlabeled 100 10 1 10 2 0.0 0.2 0.4 0.6 0.8 1.0 Score Unweighted_Train 10 Decay Muon Noise Normalized counts 8 Cascade Afterpulse MixedOnlyBkg 6 MixedDecayBkg Unlabeled 4 2 0 0.0 0.2 0.4 0.6 0.8 1.0 Score Figure 4.5: Score distributions for the decay pulse classifier. The top image presents the total pulses per year expected as a function of score for each pulse type, whereas the bottom image presents normalized distributions for each pulse type. The blue distribution corresponds to decay pulses. Notice this distribution increases with score, as expected. Backgrounds tend to eventually decrease at high scores, with the exception of afterpulses, the main background for decay pulse extraction. In addition, the top figure shows that the extraction of decay pulses is negatively impacted by their total rate, significantly lower than that of other contributions. Future analyses will require both higher decay signal rates as well as better discrimination variables against afterpulses to be able to exploit decay signals. 68 although afterpulses are generated from precursor hits, is possible for their precursor to not pass the discriminator threshold to start a data readout instance (launch) in a DOM (0.25 PE). Thus, the original pulse fails to start DOM readout, but afterpulses can be generated with higher charges than their precursor. As a result the precursor is lost in some cases, diminishing the effectiveness of discriminating afterpulse events. 4.2 Neutrino/Antineutrino classifier Although it was determined that extraction of decay pulses in DeepCore events is currently marginal, a neutrino/antineutrino classifier based on this information was implemented to demon- strate the capability of the technique for future analyses in the IceCube Upgrade [81], where decay signal extraction will be enhanced by: an increased optical module density, increasing the probabil- ity for decays to occur in the proximity of one or more DOMs, and; improved variables for rejection of afterpulses in Upgrade modules, i.e. using the fact that a decay close to a mDOM generates several coincident hits within the same time window, whereas afterpulses in different modules are generated at random times. The ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ classifier leverages the limited decay signal extraction power from the classifier described in the previous section. To do this, time information from the most decay-like pulses in each event was used. However, given the low level of discrimination provided by this method, an exploration was performed to identify other features, beyond decay times, that could help separate neutrino and antineutrino signals. Differences in the reconstructed lengths for the two event types were ultimately identified and attributed to differences in the expected inelasticity distributions for neutrinos and antineutrinos that arise at low energies (see section 2.6). Figure 4.6 shows the true (figure 4.6a) and reconstructed (figure 4.6b) inelasticity distributions in the full baseline analysis sample prior to the final cuts at level 7. An advantage of the inelasticity parameter is the fact that it can be reconstructed in every single neutrino event, while decay signals are only recorded in a subset of events. However, the difference in the expected number of neutrino and antineutrino events in IceCube (about 2:1) results in an impossibility to utilize inelasticity to obtain a sub-sample where antineutrinos are dominant. It is instead only possible to obtain sub-samples with a higher 69 (a) (b) Figure 4.6: Inelasticity distributions for ๐œˆ ๐œ‡ and ๐œˆยฏ ๐œ‡ in the analysis sample, both from retrieving the true Monte Carlo information (a) and from the reconstructed cascade and track energies using LEERA. or lower ๐œˆ ๐œ‡ dominant population than in the total sample. This is also true for discrimination using muon decay signals, and is a feature that is ultimately transferred to the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ classifier. The full set of per event training features for the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT classifier are: โ€ข The time delay between the pulse detection and the neutrino interaction for the first, second and third most decay-like pulses in each event (3 features). โ€ข The charge of the first, second and third most decay-like pulses in each event (3 features). โ€ข The decay probability score for the first, second and third most decay-like pulses in each event (3 features). โ€ข Reconstructed inelasticity โ€ข Reconstructed length of the outgoing lepton (muon) โ€ข Quality metric for LEERA reconstruction 70 Most decay-like pulse per event Second most decay-like pulse per event Third most decay-like pulse per event 102 102 102 Probability density Probability density Probability density 101 101 101 100 100 100 10 1 10 1 10 1 10 2 10 2 10 2 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 pulse charge (PE) pulse charge (PE) pulse charge (PE) Most decay-like pulse per event Second most decay-like pulse per event Third most decay-like pulse per event 10 2 10 2 10 2 10 3 Probability density Probability density Probability density 10 4 10 3 10 3 10 5 10 4 10 4 10 6 10 5 10 5 10 7 10 8 10 6 10 6 10 9 10 7 10 7 2500 0 2500 5000 7500 10000 2500 0 2500 5000 7500 10000 2500 0 2500 5000 7500 10000 tpulse tinteractionvertex (ns) tpulse tinteractionvertex (ns) tpulse tinteractionvertex (ns) Most decay-like pulse per event Second most decay-like pulse per event Third most decay-like pulse per event 102 102 102 Probability density Probability density Probability density 101 101 101 100 100 10 1 100 10 2 10 1 10 3 10 1 10 2 10 4 0.0 0.2 0.4 0.6 0.8 1.0 10 3 0.0 0.2 0.4 0.6 0.8 1.0 10 2 0.0 0.2 0.4 0.6 0.8 1.0 decay BDT score decay BDT score decay BDT score 10 1 Probability density 10 2 10 3 10 4 10 5 0 200 400 600 reco muon length (m) Figure 4.7: Distributions of pulse time, decay BDT score and charge for the three most decay-like pulses per event. Small discrimination power from decay signals is only visible for the score and charge variables related to the most decay-like pulse per event and for these pulses the time distribution is the closest to an exponential decay curve. Distributions for reconstructed length are also shown and is noted that this variable provides a modest discrimination power. 71 (a) (b) Figure 4.8: Schematic showing (a) BDT scores for the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ classifier for both true ๐œˆ ๐œ‡ and true ๐œˆยฏ ๐œ‡ populations and (b) fraction of ๐œˆยฏ ๐œ‡ events found when splitting the sample in two bins using a single score threshold. Scores represent a probability that the particle type that generated the event was a ๐œˆยฏ ๐œ‡ . The threshold selected for the analysis was at 0.605. This value provided different ๐œˆ ๐œ‡ : ๐œˆยฏ ๐œ‡ ratios in the two identified bins (above and below the threshold), while maintaining significant MC statistics in each bin. The last feature, the quality metric for LEERA (the reconstruction algorithm used in the analysis to estimate interaction vertex, outgoing lepton path length and event energy), is simply the hit/no hit likelihood for the LEERA reconstructed particle (see Appendix F.4), included since events where the reconstruction has a high likelihood have generally better inelasticity estimators. Distributions of these variables are shown in figures 4.6b and 4.7, and reflect that inelasticity and length currently provide the largest separation power for ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ discrimination. The classifier was a BDT, trained in 70% of all ๐œˆ ๐œ‡ + ๐œˆยฏ ๐œ‡ events, and tested in the remaining 30%. The resulting classifier calculates a probability score that the particle type that generated the event was a ๐œˆยฏ ๐œ‡ . The distribution of these scores can be found in figure 4.8. No significant improvement in performance was found in the testing set when compared to the training set, as shown in figure 4.9. As a result, the full MC sample was used to perform the NMO analysis. A simple test performed with this classifier was to test whether ๐œˆ๐‘’ and ๐œˆ๐œ events showed any classification power. It was found that this was not really the case (see table 4.1), and thus, the classifiers only separated ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ events, as expected. The binning based on this classifier that was ultimately used in the NMO analysis was based on a single cut in the probability score distribution. The cut was at 0.605, and provided two samples with different ๐œˆ ๐œ‡ : ๐œˆยฏ ๐œ‡ ratios (events with scores 72 Figure 4.9: ROC curves for the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ classifier. Similar performances were obtained with the training and testing samples. As a result, the full sample was used in the NMO analysis with ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ discrimination. Event type ๐œˆ : ๐œˆยฏ (below threshold) ๐œˆ : ๐œˆยฏ (above threshold) ๐œˆ๐‘’ CC 0.687 : 0.313 0.681 : 0.319 ๐œˆ ๐œ‡ CC 0.710 : 0.290 0.668 : 0.332 ๐œˆ๐œ CC 0.721: 0.279 0.705 : 0.295 ๐œˆ NC 0.751 : 0.249 0.756 : 0.244 Table 4.1: Table describing the fractions ๐œˆ : ๐œˆยฏ above and below the score threshold at 0.605. Maximum change in the ratios is found for ๐œˆ ๐œ‡ CC events, as expected. above and below this threshold). A summary of the ratios in the two bins for each neutrino type can be found in table 4.1. As can be noted, the difference is only marginal and as will be shown in chapter 5 it only provides an small improvement to the NMO sensitivities. 73 CHAPTER 5 THE DEEPCORE NEUTRINO MASS ORDERING ANALYSIS This chapter presents an overall description of the NMO analysis. The methodology is based on a previous NMO analysis performed with three years of DeepCore data [3]. The NMO is measured by selecting a highly pure sample of atmospheric neutrino events in Monte Carlo simulations and experimental data (see appendix E), applying the same selection criteria, and comparing event distributions in simulated data, predicted under each mass ordering to those found with experimental data. Comparisons are performed using binned distributions representing event counts as a function of reconstructed energy and path length (zenith angle, given equation 2.35). In addition, since the oscillation patters are different for ๐œˆ๐‘’ , ๐œˆ ๐œ‡ and ๐œˆ๐œ , the analysis presented in [3] also used a third binning dimension, PID, that roughly separated cascade-like and track-like events. The difference in the oscillation patterns for the NO and the IO arises primarily from the fact that matter effects impact neutrinos in the NO and antineutrinos in the IO. Differences in the oscillation patterns in previous analyses relied on differences in the cross section and flux of neutrinos and antineutrinos, resulting in an event ratio ๐œˆ : ๐œˆ โ‰ˆ 2 : 1 [3]. As a result, matter effects were more prominent in the combined distributions for ๐œˆ + ๐œˆ in the NO, allowing for a sensitivity to the NMO. It is expect that the sensitivity would be maximized when the distributions for ๐œˆ/๐œˆ are separable. In the analysis presented here, a measurement using the additional binning for ๐œˆ ๐œ‡ events, based on the ๐œˆ/๐œˆ classifier introduced in chapter 4, is explored. 5.1 Analysis method The analysis presented in this work was based on a physics analysis software package in IceCube called PISA [4]. PISA was used to retrieve expected atmospheric neutrino event distributions (plus background contributions) in the IceCube-DeepCore detectors given a set of neutrino oscillation parameters. The software calculates these distributions by re-weighting MC events. Several additional systematic parameters that affected event expectations were also added via weighting, including those related to: atmospheric neutrino fluxes; modeling of neutrino-nucleon interactions; optical properties of the Antarctic ice and the refrozen hole ice; the overall detection efficiency of 74 the DOMs, and; the overall ratio between atmospheric muon and atmospheric neutrino fluxes. Event expectations found with PISA may be utilized to realize desired physical observables. For the purpose of the NMO analysis the observables of interest are the reconstructed energy and reconstructed zenith of the neutrino events since these variables allowed probing into the oscillation probabilities, as shown in figure 2.9. Another observable of interest is PID since generally the oscillation probabilities for tracks and cascades are different (each neutrino flavor has different oscillation probabilities). A final observable employed was neutrino/antineutrino scoring since oscillation probabilities for these particle types are different under each neutrino mass ordering for atmospheric neutrinos (given matter effects as they crossed the Earth). In the NMO analysis performed in this work, oscillation parameters in experimental (or simu- lated) data were fit with PISA, using the following procedure: โ€ข Select a binning over each of the observables of interest (energy, zenith, PID, neutrino/antineutrino scoring). โ€ข Obtain the multi-dimensional data distribution from experimental or simulated data using the identified binning. โ€ข Define a set of nominal values ( ๐‘ยฎnominal ) that reflects current best knowledge of the oscillation and nuisance parameters in the fit. โ€ข Obtain test distribution from PISA under the nominal set of parameters. โ€ข Compare test and data distributions using an appropriate metric, e.g. a Poisson log-likelihood or a Pearsonโ€™s ๐œ’2 . โ€ข Use a minimization algorithm (SLSQP [82]) to determine the physics and nuisance parameters for the defined metric. Use PISAโ€™s weighting scheme to define distributions at the new values explored by the minimizer. โ€ข Once the convergence criteria for the minimizer is satisfied, best fit values for the oscillation parameters are determined. 75 In addition, prior knowledge on the oscillation probability space is used to further aid the fitting process. For example, it is known that there is a degeneracy in the oscillation probabilities given the octant of the ๐œƒ 23 mixing angle (see figure 2.10). As a result, the likelihood space contains two local minima, one in each octant. To account for this, fitting is performed independently within each octant and, subsequently, the global minima is identified. In a similar fashion, the normal 2 > 0) and inverted (ฮ”๐‘š 2 < 0) ordering hypotheses were fit individually. There are a total (ฮ”๐‘š 32 32 of 4 fits performed to each dataset in order to determine the best represented set of parameters. A measurement of the NMO was performed in this work by fitting experimental data to both the NO and the IO using PISA and determining the best fit ordering. Both the experimental data and the MC simulation data passed to PISA were previously processed to describe the analysis sample of neutrino events. Several details related to the fitting scheme using PISA will be covered in this chapter, including: the treatment of systematic parameters in the fits and the nominal values used to weight the MC events with PISA (subsection 5.1.1); the binning used in the analysis (subsection 5.1.2); and the modified Poisson likelihood function that was used to perform the fits (subsection 5.1.3). In addition to the fits to the NMO, a hypothesis test was performed to determine the statistical significance of the result. In the context of the NMO, this is a non-nested binary hypothesis test, to determine, given a best fit hypothesis (denoted true ordering, or TO), the rejection power against the alternative hypothesis (denoted wrong ordering, or WO). The hypothesis test followed the methodology used in a previous NMO analysis [3] and is introduced and explained in detail in subsection 5.1.4. Preliminary tests to the analysis using only MC distributions are reviewed in section 5.2, and the final analysis results with experimental data are presented in section 5.3. 5.1.1 Nominal set and treatment of systematics The PISA framework provides the means to extract oscillation parameters from experimental (or MC) data by fitting event expectations in the data to those obtained by re-weighting using PISA. The physics and systematic parameters considered in this work included: โ€ข Atmospheric flux systematic parameters (introduced in this section). 76 โ€ข Oscillation parameters: ฮ”๐‘š 21 2 , ฮ”๐‘š 2 , ๐œƒ , ๐œƒ , ๐œƒ , ๐›ฟ . 32 12 12 23 ๐ถ๐‘ƒ โ€ข Detector systematic parameters: ๐‘ 0 , ๐‘ 1 , DOM efficiency factor (see Appendix D.4). โ€ข Neutrino-nucleon interaction systematics: ๐‘€ ๐ด ๐ถ๐ถ ๐‘…๐ธ ๐‘† and ๐‘€ ๐ถ๐ถ ๐‘„๐ธ ๐ด (see Appendix D.1). โ€ข Normalization of the distributions calculated for NC events (referred to as NC norm) โ€ข Overall normalization of the neutrino distribution obtained by adding up all NC and CC components (๐ด๐‘’ ๐‘“ ๐‘“ ) โ€ข Overall normalization of the muon distribution (referred to as the weight scale) PISA calculates the distributions for neutrinos and muons separately, and multiplies a normal- ization factor to each before adding them. Muon distributions required the generation of dedicated Kernel Density Functions from the raw MC counts, given their low statistics [2]. As mentioned in section 2.5.2, the nominal atmospheric neutrino flux model used in the analysis was the HKKM model [49]. The atmospheric neutrino flux is ultimately dependent on the hadronic model used to describe cosmic ray air showers since this choice changes the production rates of pions and kaons (and thus neutrinos). Uncertainties in the hadronic model were introduced in the analysis using a modified version of the Barr scheme [83]. In this scheme, uncertainties in the production rates for kaons and pions were estimated for different phase space regions in energy of the incident parent particle (๐ธ๐‘– ) and fraction of energy transferred to the meson secondary (๐‘ฅ ๐‘™๐‘Ž๐‘ ). A reduced version of this scheme was defined by the collaboration for oscillation studies using the analysis sample chosen for the NMO analysis presented in this work, where phase space regions with similar effects were merged [2]. The kaon uncertainties were defined independently for ๐พ ยฑ , whereas ๐œ‹ ยฑ both used the same uncertainties and differences in their production rates were estimated by introducing a single parameter describing the uncertainty in the ๐œ‹ + /๐œ‹ โˆ’ ratio. In addition, uncertainties related to the primary cosmic ray spectrum are added by multiplying an energy dependent correction weight to the nominal neutrino fluxes (and added to the event 77 weights, as in equation D.6), given by:  ๐ธ  ฮ”๐›พ ฮ”ฮฆ๐‘›๐‘œ๐‘š = (5.1) ๐ธ ๐‘๐‘–๐‘ฃ๐‘œ๐‘ก where ๐ธ is the neutrino energy and ๐ธ ๐‘๐‘–๐‘ฃ๐‘œ๐‘ก is fixed (โˆผ24.1 GeV). The corrections to the atmospheric flux introduced by the uncertainties in meson production ratios and the spectral index of the cosmic ray spectrum modified the nominal flux weights according to the following relation: โˆ‘๏ธ  ๐‘‘ฮฆ  ฮฆ๐‘š๐‘œ๐‘‘ = (ฮฆ๐‘›๐‘œ๐‘š ยท ฮ”ฮฆ๐‘›๐‘œ๐‘š ) + ๐‘๐‘– , (5.2) ๐‘‘๐‘๐‘– ๐‘– where ฮฆ๐‘›๐‘œ๐‘š is the nominal flux, the summation is over each phase space region (and meson type) in the reduced Barr scheme, and b is the value of the uncertainty in the i-th Barr region. The derivatives of the fluxes with respect to each systematic ๐‘๐‘– were previously calculated, interpolated and stored in order to allow fast calculation of flux weights in PISA, given a particular set of flux systematic parameters (ฮ”๐›พ plus the Barr systematics). Unoscillated fluxes are corrected in PISA to include the effects of neutrino oscillations. The nominal set of oscillation parameters used in the NMO analysis is drawn from the global fits performed by the NuFit collaboration in 2018 [19, 20] (see table 2.1). However, the MC expectations were fixed to assume no CP violation (๐›ฟ๐ถ๐‘ƒ = 0) since the corrections to the oscillation parameters introduced by this parameter are negligible in the analysis energy range, as discussed in section 2.5.3. Given the low sensitivity of atmospheric neutrino experiments to the solar oscillation parameters 2 were also kept fixed when calculating event count expectations with (the 1-2 sector), ๐œƒ 12 and ฮ”๐‘š 21 PISA. Treatment of detector systematics require dedicated systematic MC sets. For all other systemat- ics, weights were applied to each event to modify the event distributions calculated with PISA. For detector systematics the binning for the analysis distributions were determined in advance. PISA was then used to calculate the expected counts in the detector assuming the nominal values for all non-detector systematics using each of the detector systematic MC sets described in Appendix D.4. Each MC set returned different counts in each of the analysis bins and the variations in the bin 78 counts were then fit to a linear relation: 2 โˆ‘๏ธ ๐‘“ ( ๐‘0, ๐‘1, ๐‘2) = ๐‘ + ๐‘š๐‘– ฮ”๐‘๐‘– (5.3) ๐‘–=0 where the DOM efficiency correction factor is described by the ๐‘ 2 parameter, ฮ”๐‘๐‘– represents the systematic variation of each parameter with respect to its nominal value, ๐‘š๐‘– are the slopes of the function over each systematic parameter, and c is an offset. The factor ๐‘“ ( ๐‘ 0 , ๐‘ 1 , ๐‘ 2 ) is ultimately a re-weighting factor per bin, that modifies counts calculated by PISA using the baseline MC set to expected counts given a new set of values ๐‘ 0 , ๐‘ 1 and ๐‘ 2 for the detector systematics. The set of all re-weighting factors over all bins is referred to as the systematic "hypersurfaces". These hypersurfaces were further calculated for different values of the physics parameter ฮ”๐‘š 32 2 since it was found in previous studies [2] that hypersurfaces calculated with the analysis sample chosen for the NMO analysis varied as a function of this parameter. As a result, the hypersurfaces are interpolated over different ฮ”๐‘š 32 2 values. In addition, this process is performed for the NO and the IO independently since the fit range ultimately chosen for ฮ”๐‘š 32 2 was not continuous. The hypersurface fits performed using the analysis binning for the NMO analysis presented in this work are discussed in section 5.1.2. Finally, conservative assumptions on the validity range for the neutrino and muon normalization factors were used to define prior ranges for these systematics. Table 5.1 showcases the nominal values used for each of the systematic parameters considered in the NMO analysis, as well as their priors and the ranges to which the fitted values were restricted. Previous studies on the impact of the systematic parameters on a muon neutrino disappearance analysis performed with the analysis sample [2] were used as a basis when deciding what parameters to keep fixed and which to vary when measuring the NMO. 5.1.2 Analysis binning When choosing a histogram binning with PISA, notice that the finer the binning the better the observable space (energy, direction, PID and nunubar scores) can be probed to identify features of the oscillation probabilities. This comes with an inherent limitation, however, since the reliability of the expectations provided by PISA will be small once the number of MC events within each bin 79 Parameter Nominal value Fit range Fixed? Description (ยฑ prior) ฮ”๐›พ 0 ยฑ 0.1 nominal ยฑ 5๐œŽ NO Shifting of atm. ๐œˆ spectral index ๐œ‹ + /๐œ‹ โˆ’ 0 ยฑ 0.05 nominal ยฑ 5๐œŽ YES Pion production ratio uncer- tainty Barr af Pi 0.0 ยฑ 0.63 nominal ยฑ 5๐œŽ NO ๐œ‹ + production uncertainty in phase space regions A-F Barr g Pi 0.0 ยฑ 0.03 nominal ยฑ 5๐œŽ NO Like Barr af Pi, but in region G Barr h Pi 0.0 ยฑ 0.15 nominal ยฑ 5๐œŽ NO Like Barr af Pi, but in region H Barr i Pi 0.0 ยฑ 0.112 nominal ยฑ 5๐œŽ YES Like Barr af Pi, but in region I Barr w K 0.0 ยฑ 0.4 nominal ยฑ 5๐œŽ NO ๐พ + production uncertainty in phase space region W Barr x ๐พ + 0.0 ยฑ 0.1 nominal ยฑ 5๐œŽ YES Like Barr w ๐พ + , but in region X Barr y ๐พ + 0.0 ยฑ 0.3 nominal ยฑ 5๐œŽ NO Like Barr w ๐พ + , but in region Y Barr z ๐พ + 0.0 ยฑ 0.122 nominal ยฑ 5๐œŽ YES Like Barr w ๐พ + , but in region Z Barr w ๐พ โˆ’ 0.0 ยฑ 0.4 nominal ยฑ 5๐œŽ NO ๐พ โˆ’ production uncertainty in phase space region W Barr x ๐พ โˆ’ 0.0 ยฑ 0.1 nominal ยฑ 5๐œŽ YES Like Barr w ๐พ โˆ’ , but in region X Barr y ๐พ โˆ’ 0.0 ยฑ 0.3 nominal ยฑ 5๐œŽ YES Like Barr w ๐พ โˆ’ , but in region Y Barr z ๐พ โˆ’ 0.0 ยฑ 0.122 nominal ยฑ 5๐œŽ YES Like Barr w ๐พ โˆ’ , but in region Z ๐ด๐‘’ ๐‘“ ๐‘“ 1 [0,3] NO ๐œˆ (all) distribution normalization Weight scale 1 [0,3] NO ๐œ‡ distribution normalization NC norm 1 ยฑ 0.2 nominal ยฑ 5๐œŽ NO ๐œˆ NC distribution normalization ๐‘€๐ด ๐ถ๐ถ ๐‘…๐ธ ๐‘† 0.0 ยฑ 1.0 nominal ยฑ 2๐œŽ NO Variation to nominal axial mass ๐ถ๐ถ ๐‘„๐ธ ๐‘€๐ด 0.0 ยฑ 1.0 nominal ยฑ 2๐œŽ NO Variation to nominal axial mass ๐‘0 0.101569 [-2,1] NO Refrozen hole ice systematic ๐‘1 -0.049344 [-0.2,0.2] NO Refrozen hole ice systematic DOM eff 1.0 ยฑ 0.1 [0.8, 1.2] NO DOM efficiency correction fac- tor ฮ”๐‘š 212 /eV2 7.39 ร— 10โˆ’5 - YES Same for both orderings ๐œƒ 12 /โ—ฆ 33.82 - YES Same for both orderings ๐œƒ 13 /โ—ฆ (NO) 8.61 ยฑ 0.13 [8.22, 8.98] YES Only for fits over NO ๐œƒ 13 /โ—ฆ (IO) 8.65 ยฑ 0.13 [8.27, 9.03] YES Only for fits over IO ๐œƒ 23 /โ—ฆ 49.7 [0, 90] NO Same for both orderings ฮ”๐‘š 312 /eV2 (NO) 2.525 ร— 10โˆ’3 [1, 7] ร—10โˆ’3 NO Only for fits over NO ฮ”๐‘š 312 /eV2 (IO) โˆ’2.438 ร— 10โˆ’3 [-7, -1] ร—10โˆ’3 NO Only for fits over IO Table 5.1: The nominal values of the physics and systematic parameters accounted for when fitting with PISA. Appropriate Gaussian priors and fit ranges were introduced for several parameters. Definitions of the different phase space regions for Barr parameters can be found in [83]. 80 falls below a certain threshold. If the binning is too fine, the expected patterns predicted may be driven by statistical fluctuations rather than by physical features of the oscillation probabilities. In this case the hypersurface fits also fail to converge, failing goodness of fit tests, and the PISA fit parameters will be unreliable. Figure 5.1 shows 2D distributions in energy and zenith obtained with PISA for the nominal NO parameters in table 5.1. Instead of using reconstructed PID, the figures show distributions for each of the true neutrino event types. Notice oscillation patterns are largely visible for event counts related to muon neutrinos and oscillation patterns for electron neutrinos are subdominant, and thus, less visible. Finally, tau neutrinos are expected mainly from ๐œˆ ๐œ‡ โ†’ ๐œˆ๐œ oscillations, and are detected at a much lower rate. The main oscillation channel probed with the analysis sample is therefore muon neutrino disappearance. Muon neutrino oscillations can be recovered using the distributions in figure 5.1, although the resolutions of the energy and directional reconstructions reduce the power to resolve some features of the oscillation patterns. However the uncertainties in the predicted values per bin are high for the binning scheme of figure 5.1, and there are several empty bins making hypersurface fits not possible over all bins. The PISA binning used for this NMO analysis was imposed two minimum conditions: โ€ข Fits to PISA distributions generated with some set of injected parameters, ๐‘, ยฎ should converge on average to the injected values. โ€ข Hypersurface fits must be successfully performed over all analysis bins. The binning of the energy, zenith and the ๐œˆ/๐œˆยฏ scoring distributions was defined prior to, and remained fixed, for the analysis. A 4D binning including reconstructed PID is not explored in this work, but could be explored in a future iteration of the analysis. A 2D version of this binning, containing only energy and direction, was verified to satisfy the two minimal conditions such that the NMO sensitivities obtained with and without ๐œˆ/๐œˆยฏ scoring could be properly compared. In general, the goal of the analysis presented in this work was to show 81 (a) (b) Figure 5.1: Nominal Normal Ordering Monte Carlo distributions predicted by PISA using truth (a) and reconstructed (b) variables. 82 the potential for improvement of NMO sensitivities by adding ๐œˆ/๐œˆยฏ scoring in order to aid future NMO analyses with the IceCube Upgrade. Optimization of the current sensitivities to the NMO obtained with DeepCore was not explored. The following binning was ultimately chosen for each dimension considered: โ€ข log10 (๐ธ): 7 uniform bins between 6.31 GeV and 158.49 GeV, with the last two bins merged. โ€ข cos(๐œƒ): 6 uniform bins between -1 and 0.1 โ€ข ๐œˆ/๐œˆยฏ score: 2 bins, for values above and below a score of 0.605. Hypersurface fits were successfully obtained over each of the analysis bins, in both cases (see figures 5.2 and 5.3). Is important to add here that the fits were performed independently over different neutrino interaction types since their event rates had a different dependence on the detector systematics. Fits were therefore performed independently for purely cascade-like events (๐œˆ๐‘’ CC + ๐œˆ๐‘Ž๐‘™๐‘™ NC), ๐œˆ๐œ CC events (different from other cascade like events, since they can produce muons when decaying) and track-like events (๐œˆ ๐œ‡ ๐ถ๐ถ). Since there were 3 fits per bin (one per event type), and a total of 36 (72) bins for the 2D (3D) scheme, there were a total of 108 (216) fits per hypersurface. In addition, a total of 40 hypersurfaces were produced to cover different ฮ”๐‘š 32 2 values, and the results were interpolated separately over the NO and the IO. To test the quality of these fits, the ๐œ’2 /dof over each fit was calculated (following the method- ology in [2]). Distributions of the reduced ๐œ’2 values were generated for each neutrino event type and it was verified that < ๐œ’2 /dof> โ‰ˆ 1 in all cases. Figures 5.2 and 5.3 show examples of these distributions for the hypersurface fits obtained using the nominal NO value of ฮ”๐‘š 32 2 for the 2D and the 3D binning, respectively. The per bin reduced ๐œ’2 values were also plotted in order to evaluate whether there were particular regions in the reconstructed variable space with a higher ratio of "bad fits" (those at the tail of the reduced ๐œ’2 distribution). Examples of the distributions obtained using the nominal NO ฮ”๐‘š 32 2 value can be found in figures 5.4 and 5.5 for the 2D and 3D binning, respectively. No 83 Figure 5.2: Reduced ๐œ’2 values obtained for all hypersurface fits performed using the nominal 2 for a binning scheme without neutrino/antineutrino scores. It can be seen that NO value of ฮ”๐‘š 32 < ๐œ’2 /dof> โ‰ˆ 1 for all event types. 84 Figure 5.3: Reduced ๐œ’2 values obtained for all hypersurface fits performed using the nominal 2 , for a binning scheme with neutrino/antineutrino scores. It can be seen that NO value of ฮ”๐‘š 32 < ๐œ’2 /dof> โ‰ˆ 1 for all event types. 85 regions with a significantly higher rate of bad fits were found and, as a result, the hypersurface fitting process was concluded. Final distributions using the 2D and 3D analysis binning (applying the nominal NO parameters) can be found in figures 5.6 and 5.7, respectively. The distributions over each true neutrino event type is individually produced to properly apply hypersurface corrections. Muon distributions are also separately generated, since these did not include systematics effects when being calculated. Distributions over all neutrino types were ultimately merged, multiplied by an overall neutrino normalization factor (๐ด๐‘’ ๐‘“ ๐‘“ ), and added to the muon template (also including an appropriate muon normalization factor, the weight scale listed in table 5.1). The final (merged) distributions are shown in figures 5.6 and 5.7. 5.1.3 Likelihood function To perform fits with PISA, event counts in experimental data distributions are compared to event counts calculated with PISA according to some set of physics and nuisance parameters, ๐‘. ยฎ We define ๐‘˜ ๐‘– to be the number of event counts in data for the i-th bin, and ๐œ‡๐‘– ( ๐‘) ยฎ the number of events in the same bin predicted using PISA. Notice that the prediction by PISA is given by the sum of the weights from the ๐‘š๐‘– events contained within the i-th bin: ๐‘š๐‘– โˆ‘๏ธ ๐œ‡๐‘– = ๐œ‡๐‘– ( ๐‘) ยฎ = ๐‘ค ๐‘— ( ๐‘). ยฎ (5.4) ๐‘— The uncertainty in this prediction, ๐œŽ๐‘– , is given by the sum of squared weights: ๐‘š๐‘– โˆ‘๏ธ ๐œŽ๐‘– = ๐œŽ๐‘– ( ๐‘) ยฎ = ๐‘ค 2๐‘— ( ๐‘). ยฎ (5.5) ๐‘— A Poisson likelihood would return the probability of observing ๐‘˜ ๐‘– events in the i-th bin, given an expectation ๐œ‡๐‘– ( ๐‘): ยฎ ยฎ ๐‘˜ ๐‘– ๐‘’ โˆ’๐œ‡๐‘– ( ๐‘) ๐œ‡๐‘– ( ๐‘) ยฎ L ( ๐‘|๐‘˜ ยฎ ๐‘– ) = ๐‘ƒ(๐‘˜ ๐‘– |๐œ‡๐‘– ( ๐‘)) ยฎ = . (5.6) ๐‘˜! This expression, however, doesnโ€™t account for limitation of MC statistics used to infer the expected counts ๐œ‡๐‘– ( ๐‘). ยฎ The analysis instead uses the modified Poisson likelihood presented in [84] that includes corrections to address limited MC statistics. Specifically, one can return to 5.6, 86 Figure 5.4: Reduced ๐œ’2 values per analysis bin for all hypersurface fits obtained using a binning without neutrino/antineutrino scores. The example corresponds to the hypersurfaces obtained using the nominal NO value of ฮ”๐‘š 322 . 87 Figure 5.5: Reduced ๐œ’2 values per analysis bin for all hypersurface fits obtained using a binning with neutrino/antineutrino scores. The example corresponds to the hypersurfaces obtained using the nominal NO value of ฮ”๐‘š 322 . 88 Figure 5.6: Nominal normal ordering distributions per true event type obtained using the 2D analysis binning. The total distributions obtained by PISA (merging the per event type distributions) can be found in the lower right corner of the figure. starting with the general expression: ๐‘˜ ๐œ‡๐‘– ๐‘– ๐‘’ โˆ’๐œ‡๐‘– โˆซ L ( ๐‘|๐‘˜ ยฎ ๐‘–) = ๐‘ƒ(๐œ†๐‘– | ๐‘คยฎ๐‘– ( ๐‘)) ยฎ ๐‘‘๐œ†๐‘– (5.7) ๐‘˜! where ๐‘คยฎ๐‘– represents the set of all weights corresponding to events within the i-th bin and ๐‘ƒ(๐œ†๐‘– | ๐‘คยฎ๐‘– ( ๐‘)) ยฎ is given by a delta function:   ๐‘ƒ(๐œ†๐‘– | ๐‘คยฎ๐‘– ( ๐‘)) ยฎ = ๐›ฟ ๐œ† ๐‘– โˆ’ ๐œ‡๐‘– . (5.8) A more realistic expression for ๐‘ƒ(๐œ†๐‘– | ๐‘คยฎ๐‘– ( ๐‘))ยฎ based on limited MC statistics is given by [84]: ๐›ฝ๐›ผ ๐‘’ โˆ’๐œ†๐›ผ ๐œ†๐›ผโˆ’1 ๐‘ƒ(๐œ†๐‘– | ๐‘คยฎ๐‘– ( ๐‘)) ยฎ = , (5.9) ฮ“(๐›ผ) where ฮ“ represents the gamma function and: ๐œ‡๐‘–2 ๐›ผ= + 1, (5.10) ๐œŽ๐‘–2 ๐œ‡๐‘– ๐›ฝ= . (5.11) ๐œŽ๐‘–2 The effective Poisson likelihood used in this work follows from equation 5.7 together with equation 5.9. However, an additional term is added to the likelihood to include Gaussian priors for 89 (a) (b) Figure 5.7: Nominal normal ordering distributions per true event type obtained using the 3D analysis binning. The "nu-like" bins (a) and the "nubar-like" bins (b) are shown separately. The total distributions obtained by PISA (merging the per event type distributions) can be found in the lower right corner of (a) and (b) for the "nu-like" and "nubar-like" bins respectively. 90 the fit parameters, of the form: โˆ‘๏ธ ( ๐‘ โˆ’ ๐‘ ) 2 0 ๐‘†= , (5.12) ๐‘โˆˆ๐‘ ๐‘ฆ๐‘ ๐‘ก ๐œŽ๐‘2 where the sum is over all parameters with Gaussian priors, p is the tested value of the parameter, ๐‘ 0 is the nominal value of the parameter and ๐œŽ๐‘ is its Gaussian prior (estimated uncertainty). As a result, in order to fit physics and systematic parameters, bins contents (๐‘˜ ๐‘– ) in data were compared to MC expectations calculated with PISA (๐œ‡๐‘– ) via the following negative log-likelihood: "โˆซ ๐‘˜ ๐‘– โˆ’๐œ‡ # โˆ‘๏ธ ๐›ฝ๐›ผ ๐‘’ โˆ’๐œ†๐›ผ ๐œ†๐›ผโˆ’1 ๐œ‡๐‘– ๐‘’ ๐‘– 1 โˆ‘๏ธ ( ๐‘ โˆ’ ๐‘ 0 ) 2 ๐ฟ๐ฟ๐ป = โˆ’ log ๐‘‘๐œ†๐‘– + . (5.13) ๐‘– ฮ“(๐›ผ) ๐‘˜! 2 ๐‘โˆˆ๐‘ ๐‘ฆ๐‘ ๐‘ก ๐œŽ๐‘2 A minimizer algorithm determines the set of parameters, ๐‘, ยฎ that minimizes this expression. 5.1.4 Determination of the Neutrino Mass Ordering The determination of the NMO can ultimately be performed using a binary hypothesis test with two non-nested hypotheses. The method follows that applied to previous IceCube NMO analyses [3], consisting of the following recipe: โ€ข Finding best fits to experimental data under both mass orderings: use PISA and the log- likelihood from equation 5.13 to find the set of parameters ( ๐‘ยฎ๐‘๐‘‚ ) that best describes experi- mental data, under NO. Repeat for IO ( ๐‘ยฎ๐ผ๐‘‚ ). โ€ข Compare the log-likelihoods obtained when fitting for NO (LLHNO ) and IO (LLHIO ): the minimum likelihood corresponds to the preferred NMO according to the experimental data. โ€ข Generate pseudo-data: generate best fit NO distributions with PISA using the NO best fit parameters from the previous step. Modify the event counts in each bin drawn from a Poisson distribution with a mean given by the original event counts in the bin. This returns a pseudo- experiment that would be obtained if the best fit NO distribution was the true underlying distribution in data. Repeat this procedure N times to obtain a set of N pseudo-experiments under the NO hypothesis. Finally repeated the previous steps starting from the IO best fit parameters, to obtain N pseudo-experiments under the IO hypothesis. 91 โ€ข Calculate test statistic values using pseudo data: The pseudo-experiments obtained from the previous step are used in a similar fashion as experimental data, that is, they are fit for both the NO and the IO. These fits were then used to define a test statistic (TS) that compares the fit metrics under each ordering: TS = 2ฮ”LLHNO-IO = 2(LLHNO โˆ’ LLHIO ). (5.14) โ€ข Calculate p-value of the TS obtained from experimental data according to each TS distribu- tion: Using the TS distributions from the IO and the NO pseudo-experiments, it was possible to determine the p-value of the TS value obtained with experimental data, given each distri- bution. These p-values may then be used to determine the significance of the preference of the experimental data to the best fit ordering. The above steps were developed and tested on MC data only (blind analysis) to avoid introducing bias in the results (see Section 5.2) prior to performing a final measurement using experimental data (see Section 5.3). 5.2 Asimov sensitivities to the Neutrino Mass Ordering To test the methodology chosen for the measurement of the NMO, a MC-only study was first conducted with both 2D and 3D analysis binning. As noted below each case returned a different sensitivity to the NMO, with the 3D binning marginally improving sensitivity. The MC-only analyses presented here will be referred to as the Asimov sensitivity tests throughout the rest of the text. The Asimov sensitivity tests followed the method outlined subsection 5.1.4, with the exception that experimental data was replaced by the distribution obtained from PISA using the parameter values corresponding to nominal NO: ๐‘ยฎNO,nominal . This distribution was fitted to both the NO and the IO hypotheses. The returned best fit hypothesis was the NO with the retrieved values essentially ๐‘ยฎNO,nominal . There was, however, also a set of parameters ๐‘ยฎIO,best fit , corresponding to the best fit under IO. The distributions obtained using ๐‘ยฎNO,nominal and ๐‘ยฎIO,best fit were used to generate pseudo-experiments (via Poisson fluctuations). Each of these pseudo-experiments were fit to both 92 the NO and the IO using PISA. These fits were then used to perform an ensemble test, where it was tested whether the best fit parameters recovered from each pseudo-experiment fit returned, on average, the injected ๐‘ยฎNO,nominal (for NO pseudo-data) or ๐‘ยฎIO,best fit (for IO pseudo-data) values. The results of these tests can be found in figures 5.8, 5.9 and 5.10 for the 3D binning case. Similar results were obtained when repeating the test with 2D binning. Is important to understand prior ranges shown, as they represent Gaussian priors only for the systematics and physics parameters with defined Gaussian priors in table 5.1, whereas for all other parameters they represent fit ranges (hard cuts on the allowed values for the parameters in the fits). The injected systematic values were recovered in all cases (see figures 5.8 and 5.9). However, this was not the case for the physics parameters considered in the NMO analysis. The results for the fitted physics parameters are presented in figure 5.10 and injected parameters are not always recovered. However, this was expected, given the nature of the fits. For example, fits over ๐œƒ 23 (figure 5.10a) show three peaks, one near the injected values ๐œƒ 23,inj , another around (90-๐œƒ 23,inj ) and a peak at maximal mixing (๐œƒ 23 = 45โ—ฆ ). The two peaks at either side of ๐œƒ 23 = 45โ—ฆ are the result of the octant degeneracy in the vacuum oscillation probabilities. The peak at maximal mixing is a consequence of the broad binning. Since the oscillation features cannot be resolved fully, the width of the distributions about each octant is quite pronounced. For pseudo-experiments at the tails of the distributions the fits can only estimate more or less mixing of the neutrino flavors, and as discussed in [50], since maximal mixing corresponds to the ๐œƒ 23 = 45โ—ฆ , pseudo-experiments where the fits require more mixing to account for the observed distributions will remain at maximal mixing, and produce the observed peak. 2 parameter (see figure 5.10b), the fits cannot fully resolve the NMO. As a In the case of the ฮ”๐‘š 31 result, data injected with the NO will sometimes have a best fit corresponding to an IO hypothesis, and vice-versa. As a result, the peak of the NO (IO) distribution is shifted towards the IO (NO) injected value. Fits returning the injected systematic values, on average, was important since it verified that the defined 2D and 3D binning chosen for the analysis contained enough information to (roughly) 93 350 true NO true NO true NO true IO 300 true IO 600 true IO Number of pseudoexperiments Number of pseudoexperiments 300 Number of pseudoexperiments prior prior prior 250 500 250 200 400 200 150 150 300 100 100 200 50 50 100 0 1.0 0.5 0.0 0.5 1.0 0 1.0 0.5 0.0 0.5 1.0 0 0.6 0.8 1.0 1.2 1.4 MCCQE A MCCRES A NC norm true NO 300 true NO true NO 300 true IO true IO 1000 true IO Number of pseudoexperiments Number of pseudoexperiments Number of pseudoexperiments prior prior prior 250 250 800 200 200 600 150 150 400 100 100 50 50 200 0 2.0 1.5 1.0 0.5 0.0 0.5 1.0 0 0.2 0.1 0.0 0.1 0.2 0 0.7 0.8 0.9 1.0 1.1 1.2 1.3 hole ice p0 hole ice p1 DOM eff 300 true NO 800 true NO true IO 700 true IO Number of pseudoexperiments Number of pseudoexperiments 250 prior prior 600 200 500 150 400 300 100 200 50 100 0 0 1 2 3 0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 weight scale Aeff scale Figure 5.8: Ensemble tests using the 3D analysis binning for detector, neutrino interaction and overall normalization systematics. The three distributions in the upper row included Gaussian priors. A particular note is required for the weight scale. Notice that the muon normalization has a pronounced width, but since the minimum physical value for this parameter is zero, the tail of the distribution piles up at zero. No weight scale value was found below or above the allowed fit ranges. The apparent leakage beyond the fit range was a result of the binning used to produce the weight scale distribution. 94 300 true NO 300 true NO true NO true IO true IO 300 true IO 250 Number of pseudoexperiments Number of pseudoexperiments Number of pseudoexperiments prior prior prior 250 250 200 200 200 150 150 150 100 100 100 50 50 50 00.75 0.50 0.25 0.00 0.25 0.50 0 0.3 0.2 0.1 0.0 0.1 0.2 0.3 0 0.15 0.10 0.05 0.00 0.05 0.10 0.15 Barr af Pi Barr g Pi Barr h Pi 300 true NO 300 true NO 300 true NO true IO true IO true IO 250 Number of pseudoexperiments Number of pseudoexperiments Number of pseudoexperiments prior prior prior 250 250 200 200 200 150 150 150 100 100 100 50 50 50 0 0.4 0.2 0.0 0.2 0.4 0 0.4 0.2 0.0 0.2 0.4 0 0.4 0.2 0.0 0.2 0.4 Barr w K Barr y K Barr w antiK 300 true NO true IO Number of pseudoexperiments prior 250 200 150 100 50 0 0.10 0.05 0.00 0.05 0.10 Figure 5.9: Ensemble tests using the 3D analysis binning for atmospheric neutrino flux systematics. fit all systematic and physics parameters included. As a result, a finer binning condition was not explored. The NO and the IO best fits to each pseudo-experiment were used to calculate TS values following equation 5.14. The TS distributions obtained with the 2D and 3D binning are shown in figures 5.11a and 5.11b, respectively. Notice that the TS definition from equation 5.14 is defined such that pseudo-experiments with a preference for the NO (those generated using ๐‘ยฎNO,nominal ) will have on average a negative TS since they will tend to have NO as their preferred ordering. Similarly, pseudo-experiments with a preference for the IO, those generated using ๐‘ยฎIO,best fit , will 95 175 Number of pseudoexperiments Number of pseudoexperiments injected NO injected NO injected IO 250 injected IO 150 true NO true NO true IO true IO 125 200 100 150 75 100 50 50 25 0 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0 2.0 2.2 2.4 2.6 2.8 3.0 sin2( 23) 103 | 2 | 31 (a) (b) Figure 5.10: Ensemble tests using the 3D analysis binning for the neutrino oscillation parameters. The three peak structure in the ๐œƒ 23 fits (a) can be explained by the octant degeneracy in the vacuum oscillation probabilities, as well as the pronounce width of the distributions. The tails of the distributions at each octant pile up at maximum mixing. The shifted peaks in the ฮ”๐‘š 31 2 distributions (b) arise from the limited sensitivity to the NMO. In some cases the NO pseudo- experiments can fit best to an IO hypothesis, and vice-versa. have, on average, positive TS values. As the NO is the best fit hypothesis under the set of parameters ๐‘ยฎNO,nominal , the significance of the result is drawn from the ability to reject the non-preferred (IO) hypothesis. Although, there are infinite IO hypotheses to reject (one per set of physics and nuisance parameters under the IO), the IO hypothesis that most closely recovers the event expectations drawn from the true underlying NO hypothesis is that which corresponds to the parameter values ๐‘ยฎIO,best fit . As a result, if the IO hypothesis given by the set of parameters ๐‘ยฎIO,best fit is rejected at a significance ๐›ผ, all IO hypotheses will be rejected at this significance or lower. Since the TS definition for the NMO analysis generates a NO TS distribution centered at negative values, whereas for IO the distribution is centered at positive values, p-values for rejection of the NO (IO) hypothesis are calculated using the fraction of the distribution to the right (left) of the tested TS value. The p-values obtained using this definition are also illustrated in figure 5.11. To calculate 96 2001 Pseudoexperiments (per hypothesis) 2001 Pseudoexperiments (per hypothesis) Pseudotrial sensitivity = 0.2260.068 0.072 Pseudotrial sensitivity = 0.3580.060 0.057 gaus fit (NO) IO injected (pIO = 0.410537, 0.23 ) gaus fit (NO) IO injected (pIO = 0.360013, 0.36 ) gaus fit (IO) NO injected (pNO = 0.499911) gaus fit (IO) NO injected (pNO = 0.500759) 100 100 350 300 300 250 Pseudo-Experiments Pseudo-Experiments 250 10 1 10 1 200 p-value p-value 200 150 150 10 2 10 2 100 100 50 50 0 10 3 0 10 3 2 1 0 1 2 4 3 2 1 0 1 2 3 TS TS (a) (b) Figure 5.11: TS distributions obtained with NO (orange) and IO (blue) pseudo-experiments for the Asimov test performed with the nominal NO hypothesis as the underlying true hypothesis. The solid orange (blue) curve represents the p-values to the right (left) of the NO (IO) distributions. The vertical dashed black line represents the median of the NO TS distribution. The test was performed twice using the 2D (a) and the 3D (b) analysis binning. Asimov sensitivities to the NMO under the nominal NO hypothesis were derived by finding the p-values of the WO hypothesis (IO) at the median TS for the TO (NO) hypothesis. the median sensitivity to the preferred (NO) order, the p-value of the IO TS distribution at the median of the NO TS distribution was found and converted to its sigma equivalent (using a 1-sided Gaussian z-score). Uncertainties in these sensitivities were calculated by performing Gaussian fits to the TS distributions, using a least squares algorithm. The fits returned the following parameters: โ€ข The mean of the NO and IO TS distributions (and their fit uncertainties): ๐œ‡ ๐‘๐‘‚ ยฑ ๐›ฟ(๐œ‡ ๐‘๐‘‚ ) and ๐œ‡ ๐ผ๐‘‚ ยฑ ๐›ฟ(๐œ‡ ๐ผ๐‘‚ ). โ€ข The std deviation of the NO and IO TS distributions (and their fit uncertainties): ๐œŽ๐‘๐‘‚ ยฑ ๐›ฟ(๐œŽ๐‘๐‘‚ ) and ๐œŽ๐ผ๐‘‚ ยฑ ๐›ฟ(๐œŽ๐ผ๐‘‚ ). These values could be used to estimate uncertainties in the calculated sensitivities as follows: 97 โ€ข The sensitivity uncertainty upper limit is obtained by shifting the IO TS Gaussian to the right by ๐›ฟ(๐œ‡ ๐ผ๐‘‚ ), the NO TS Gaussian to the left by ๐›ฟ(๐œ‡ ๐ผ๐‘‚ ), and reducing the width of the IO Gaussian by ๐›ฟ(๐œŽ๐ผ๐‘‚ ). These Gaussian distributions are then used to find the p-value of the IO TS distribution at the median of the NO TS distribution. The p-value is converted to its 1-sided Gaussian sigma equivalent. โ€ข The sensitivity uncertainty lower limit is obtained by shifting the IO TS Gaussian to the left by ๐›ฟ(๐œ‡ ๐ผ๐‘‚ ), the NO TS Gaussian to the right by ๐›ฟ(๐œ‡ ๐ผ๐‘‚ ), and increasing the width of the IO Gaussian by ๐›ฟ(๐œŽ๐ผ๐‘‚ ). The Gaussian distributions are then used to find the p-value of the IO TS distribution at the median of the NO TS distribution. The p-value is converted to its 1-sided Gaussian sigma equivalent. The sensitivities obtained using the above method are: โ€ข For 2D binning (no ๐œˆ/๐œˆยฏ scores): 0.226+0.072 โˆ’0.068 ๐œŽ. โ€ข For 3D binning (with ๐œˆ/๐œˆยฏ scores): 0.358+0.060 โˆ’0.057 ๐œŽ Evidently, sensitivities were modestly improved when ๐œˆ/๐œˆยฏ scores were included. Ideally, a scan over different injected values of the nominal NO hypotheses would help to understand whether the observed improvement was still present for different values of the underlying physics parameters. In addition, scans over the injected IO hypotheses were also desired, but given time constrains, these are planned for future studies. 5.3 Measurement of the Neutrino Mass Ordering using 7.5 years of Deep- Core data The measurement of the NMO was performed using 7.5 years of DeepCore experimental data, processed through the analysis neutrino selection criteria. Experimental data showed a preference for the NO using both 2D and 3D binning to perform the measurement; that is, the TS value obtained with data was negative in both instances. The TS distributions obtained from the best fit NO and IO hypotheses to experimental data for both binning scenarios are shown in figure 5.12, 98 2001 Pseudoexperiments (per hypothesis) 2001 Pseudoexperiments (per hypothesis) Data NO injected (pNO = 0.648403) Data NO injected (pNO = 0.720680) IO injected (pIO = 0.277634 IO injected (pIO = 0.170812 100 300 100 350 250 300 Pseudo-Experiments Pseudo-Experiments 250 10 1 200 10 1 p-value p-value 200 150 150 10 2 100 10 2 100 50 50 0 10 3 0 10 3 3 2 1 0 1 2 3 2 1 0 1 2 3 TS TS (a) (b) Figure 5.12: Measurement of the neutrino mass ordering without (a) and with (b) neu- trino/antineutrino binning, derived using 7.5 years of experimental data. The best fit NO and IO hypotheses to the experimental data were used to generate TS distributions. The experimental TS value is shown as a vertical solid black line and although a small preference to the NO is observed, the TS is consistent with both orderings (for both binning schemes). There is a small improvement in the results when using neutrino/antineutrino binning, apparent by a increased sep- aration between the NO and the IO TS distributions. together with the TS values obtained using experimental data. Notice that the distributions are clearly more separated when adding neutrino/antineutrino binning (as expected from the Asimov sensitivity tests). The plots also show the p-values to the right (left) of the NO (IO) TS distributions. Notice that, given the NO was the best ordering hypothesis, we expect the TS value in data to have been drawn from the NO TS distribution. A goodness of fit test can therefore be implemented where (1 โˆ’ ๐‘ ๐‘๐‘‚ ) > 0.05. This is indeed satisfied for both the 2D (figure 5.12a) and the 3D (figure 5.12b) binning. The obtained p-values at the experimental data TS under each binning are: โ€ข For 2D binning (no ๐œˆ/๐œˆยฏ scores): ๐‘ NO = 0.6484 and ๐‘ IO = 0.2776. โ€ข For 3D binning (with ๐œˆ/๐œˆยฏ scores): ๐‘ NO = 0.7207 and ๐‘ IO = 0.1708. 99 Although median sensitivities are marginally improved using 3D binning, notice that the TS in experimental data is still compatible with both the NO and IO TS distributions. A quantitative measure can be obtained to describe how much less probable it is to obtain the experimental TS value under the IO (less preferred) hypothesis. This measure is given by the ๐ถ ๐ฟ ๐‘  value, as discussed in [3]. The ๐ถ ๐ฟ ๐‘  value of interest for the results obtained in this analysis is: ๐‘ ๐ผ๐‘‚ ๐ถ ๐ฟ ๐‘  (๐ผ๐‘‚) = . (5.15) 1 โˆ’ ๐‘ ๐‘๐‘‚ This expression is derived and discussed in detail in [85]. A ๐ถ ๐ฟ ๐‘  โ‰ˆ 0 indicates a p-value under the IO much smaller than that under the NO, and thus a strong preference to the NO, whereas ๐ถ ๐ฟ ๐‘  โ‰ˆ 1 is indicative that there isnโ€™t a strongly preferred ordering. Notice the expression accounts for the fact that a very low TS value rejects both hypotheses and should therefore not be interpreted as a significant preference to the NO. The ๐ถ ๐ฟ ๐‘  values obtained for the 2D and 3D binning were: โ€ข For 2D binning (no ๐œˆ/๐œˆยฏ scores): ๐ถ ๐ฟ ๐‘  (๐ผ๐‘‚) โ‰ˆ 0.79. โ€ข For 3D binning (with ๐œˆ/๐œˆยฏ scores): ๐ถ ๐ฟ ๐‘  (๐ผ๐‘‚) โ‰ˆ 0.61. Therefore, in both cases there was not a significant preference to the NO, although the results were marginally better when including ๐œˆ/๐œˆยฏ scores. 100 CHAPTER 6 CONCLUSION Solving the NMO problem is key to completing the picture of neutrino oscillations. Sensitivity to this parameter can be achieved in the DeepCore detector via matter effects in the oscillation probabilities of atmospheric neutrinos below โˆผ 20 GeV energies, arising for neutrinos in the NO and antineutrinos in the IO. Although the optimal method for resolving the NMO with atmospheric neutrino experiments involves searching for deviations from vacuum neutrino oscillations in pure neutrino or antineutrino samples, previous NMO measurements in DeepCore [3] utilized combined samples of neutrinos and antineutrinos since separation of these event types had not been explored in IceCube. Sensitiv- ity to the NMO was still achieved since neutrinos are detected at a higher rate than antineutrinos, given their higher cross section and atmospheric flux. As a result, the combined sample pre- sented oscillation probabilities with a larger compatibility to matter effects in the NO and a lower compatibility to these effects in the IO. In the work described here, the prospect of an improved NMO analysis with DeepCore that included ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ separation using muon decay signals was studied. It was performed on a clean sample of track-like atmospheric neutrino events in DeepCore with a total livetime of 7.5 years. The novel method proposed in this analysis to separate ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ was based on the fact that the ๐œ‡+ and ๐œ‡โˆ’ produced in ๐œˆ ๐œ‡ CC interactions have different effective livetimes in DeepCore (on the order of 400 ns), resulting from ๐œ‡โˆ’ capturing in oxygen within the detector volume. Cherenkov light produced by the decay secondaries may be used to infer these decay times with the expected detection time of these signals delayed by several hundreds of nanoseconds with respect to the Cherenkov light produced by the ๐œ‡ยฑ . The envisioned strategy to extract these signals was therefore to extract signatures of delayed photons in the DeepCore DOMs closest to the muon stopping point. Through the course of the analysis, it was found that muon decay signals couldnโ€™t be extracted with sufficient efficiency in the DeepCore detector, given the relatively low light yield and the high level of background signatures in the delayed time window. 101 The primary background identified to the muon decay signature was afterpulses; delayed pulses produced in a PMT resulting from the ionization of residual gas by signal photons. Additional backgrounds included dark noise signals in the DeepCore optical modules, and scattered muon (and hadronic cascade) Cherenkov light. A BDT classifier was developed with all variables found to enhance the decay signals or mitigate the above backgrounds. While the method demonstrated potential, limited progress was made in significantly extracting the decay signals. Nevertheless, this classifier serves as an important proof of concept for future analyses in the IceCube Upgrade. In this case the decay signal extraction will be enhanced by an increased optical module density (increasing the fraction of events with detected decay signals) and improved afterpulse rejection, by leveraging the fact that decays close to multi-PMT modules produce coincident decay hits, whereas afterpulses are random processes. These studies remain ongoing utilizing simulated Upgrade detector events. An NMO analysis with ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ discrimination was performed. Inelasticity distributions also showed differences for neutrino and antineutrino events at low energies, and a ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT classifier that used both inelasticity and the most decay-like signals per event was subsequently trained. The classifier provided a small power of discrimination between ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ events. Applying the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT classifier separated the NMO analysis sample into ๐œˆ ๐œ‡ -like and ๐œˆยฏ ๐œ‡ -like events. These two classes contained different neutrino to antineutrino ratios, although both classes had a higher contribution from neutrino events given their higher detection rate in the detector. A binned maximum likelihood analysis was implemented where event distributions in the detector as a function of energy, zenith and ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT scores were used to probe into the oscillation probabilities of atmospheric neutrinos. A wide binning was used for this analysis, to ensure sufficient MC statistics were applied in each bin. In addition, an effective Poisson likelihood function, with corrections that addressed the use of limited MC statistics, was also used. Optimization of this binning to obtain maximum sensitivity to the oscillation parameters was deemed unnecessary given the low sensitivity of DeepCore to the NMO. The likelihood fits performed in the analysis also included several nuisance parameters that were used to assess the impact on the fits of uncertainties in the atmospheric neutrino and muon fluxes, the neutrino-nucleon interactions and the detector modeling. 102 The sensitivity of the NMO was first tested in a MC-only (Asimov) study, where a distribution corresponding to a nominal set of physics and nuisance parameters under the NO was used in place of experimental data. The likelihood fits found the NO as the best fit ordering, as would be expected. Pseudo-experiments were then draw from the best fit NO and IO distributions to perform a binary hypothesis test, that returned a median sensitivity for rejection of the IO. This analysis was performed using binning with and without ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT scores. The median sensitivities obtained +0.060 ๐œŽ and 0.226+0.072 ๐œŽ respectively, showing that there was a marginal increase in were 0.358โˆ’0.057 โˆ’0.068 the sensitivities in the first case. Following the development of the blind analysis, experimental data was applied in the study. Both fits with and without ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT scores showed a preference for the NO hypothesis. The fits returned ๐ถ ๐ฟ ๐‘  (๐ผ๐‘‚) values of 0.61 and 0.79, respectively. Since values close to ๐ถ ๐ฟ ๐‘  (๐ผ๐‘‚) = 0 correspond to an scenario where the IO is rejected at a significant level, including the ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ BDT scores again provided improvement in the sensitivity, but no significant determination of the NMO was obtained with the results being fully compatible with both orderings. The analysis performed in this thesis serves as a strong proof of concept for future NMO measurements in the IceCube Upgrade employing ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ discrimination. Future studies are required to test the potential for an improved ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ classifier in the IceCube Upgrade and the impact it could have on next generation NMO measurements. 103 BIBLIOGRAPHY [1] S.Blot et al. The IceCube Collaboration (internal website), oscnext (v00.06). https: //wiki.icecube.wisc.edu/index.php/OscNext. Accessed: 2022-07-24. [2] S.Blot, L. Fischer, J. Hignight, A. Trettin and W. Yan Ma. The IceCube Collaboration (internal website), oscnext verification sample analysis (v1.3). https://wiki.icecube. wisc.edu/index.php/OscNext. Accessed: 2022-07-24. [3] M. G. Aartsen et al. Development of an analysis to probe the neutrino mass ordering with atmospheric neutrinos using three years of IceCube DeepCore data. Eur. Phys. J. C, 80(1):9, 2020. [4] M. G. Aartsen et al. Computational techniques for the analysis of small signals in high- statistics neutrino oscillation experiments. Nucl. Instrum. Meth. A, 977:164332, 2020. [5] David J Griffiths. Introduction to elementary particles; 2nd rev. version. Physics textbook. Wiley, New York, NY, 2008. [6] C. L. Cowan, F. Reines, F. B. Harrison, H. W. Kruse, and A. D. McGuire. Detection of the free neutrino: a confirmation. Science, 124(3212):103โ€“104, 1956. [7] G. Danby, J. M. Gaillard, Konstantin A. Goulianos, L. M. Lederman, Nari B. Mistry, M. Schwartz, and J. Steinberger. Observation of High-Energy Neutrino Reactions and the Existence of Two Kinds of Neutrinos. Phys. Rev. Lett., 9:36โ€“44, 1962. [8] K. Kodama et al. Observation of tau neutrino interactions. Phys. Lett. B, 504:218โ€“224, 2001. [9] Q. R. et al. Ahmad. Direct evidence for neutrino flavor transformation from neutral-current interactions in the sudbury neutrino observatory. Phys. Rev. Lett., 89:011301, Jun 2002. [10] Takaaki Kajita. Atmospheric neutrino results from Super-Kamiokande and Kamiokande: Evidence for neutrino(mu) oscillations. Nucl. Phys. B Proc. Suppl., 77:123โ€“132, 1999. [11] S. Boyd. Neutrino Oscillations (Lecture Notes). University of Warwick, UK, 2016. [12] Stephen F. King and Christoph Luhn. Neutrino Mass and Mixing with Discrete Symmetry. Rept. Prog. Phys., 76:056201, 2013. [13] Ivan Esteban, M. C. Gonzalez-Garcia, Michele Maltoni, Thomas Schwetz, and Albert Zhou. The fate of hints: updated global analysis of three-flavor neutrino oscillations. Journal of High Energy Physics, 2020(9):178, 2020. [14] Standard Model. https://www.physik.uzh.ch/en/researcharea/lhcb/outreach/ StandardModel.html. Accessed: 2022-02-12. [15] T. Tantau. The TikZ and PGF Packages, Manual for version 3.0.0, 2013. http://sourceforge.net/projects/pgf/. 104 [16] P.A. Zyla et al. Review of Particle Physics. PTEP, 2020(8):083C01, 2020. [17] J. Kopp. Phenomenology of Three-Flavour Neutrino Oscillations. 2006. [18] Michael James Larson. A Search for Tau Neutrino Appearance with IceCube-DeepCore. PhD thesis, Bohr Inst., 2018. [19] Ivan Esteban, M. C. Gonzalez-Garcia, Alvaro Hernandez-Cabezudo, Michele Maltoni, and Thomas Schwetz. Global analysis of three-flavour neutrino oscillations: synergies and tensions in the determination of ๐œƒ 23 , ๐›ฟ๐ถ๐‘ƒ , and the mass ordering. Journal of High Energy Physics, 2019(1):106, 2019. [20] Nu-fit.org. v4.0: Three-neutrino fit based on data available in November 2018 | NuFIT. http://www.nu-fit.org/?q=node/177. Accessed: 2022-02-25. [21] B. T. Cleveland, Timothy Daily, Raymond Davis, Jr., James R. Distel, Kenneth Lande, C. K. Lee, Paul S. Wildenhain, and Jack Ullman. Measurement of the solar electron neutrino flux with the Homestake chlorine detector. Astrophys. J., 496:505โ€“526, 1998. [22] J. Hosaka et al. Solar neutrino measurements in super-Kamiokande-I. Phys. Rev. D, 73:112001, 2006. [23] J. P. Cravens et al. Solar neutrino measurements in Super-Kamiokande-II. Phys. Rev. D, 78:032002, 2008. [24] K. Abe et al. Solar neutrino results in Super-Kamiokande-III. Phys. Rev. D, 83:052010, 2011. [25] J. N. Abdurashitov et al. Measurement of the solar neutrino capture rate with gallium metal. III: Results for the 2002โ€“2007 data-taking period. Phys. Rev. C, 80:015807, 2009. [26] F. Kaether, W. Hampel, G. Heusser, J. Kiko, and T. Kirsten. Reanalysis of the GALLEX solar neutrino flux and source experiments. Phys. Lett. B, 685:47โ€“54, 2010. [27] B. Aharmim et al. Combined Analysis of all Three Phases of Solar Neutrino Data from the Sudbury Neutrino Observatory. Phys. Rev. C, 88:025501, 2013. [28] G. Bellini et al. Precision measurement of the 7Be solar neutrino interaction rate in Borexino. Phys. Rev. Lett., 107:141302, 2011. [29] G. Bellini et al. Measurement of the solar 8B neutrino rate with a liquid scintillator target and 3 MeV energy threshold in the Borexino detector. Phys. Rev. D, 82:033006, 2010. [30] G. Bellini et al. Neutrinos from the primary protonโ€“proton fusion process in the Sun. Nature, 512(7515):383โ€“386, 2014. [31] A. Gando et al. Reactor On-Off Antineutrino Measurement with KamLAND. Phys. Rev. D, 88(3):033001, 2013. [32] Michael Wurm. Solar Neutrino Spectroscopy. Phys. Rept., 685:1โ€“52, 2017. 105 [33] D. Adey et al. Measurement of the Electron Antineutrino Oscillation with 1958 Days of Operation at Daya Bay. Phys. Rev. Lett., 121(24):241805, 2018. [34] G. Bak et al. Measurement of Reactor Antineutrino Oscillation Amplitude and Frequency at RENO. Phys. Rev. Lett., 121(20):201801, 2018. [35] Y. Abe et al. Improved measurements of the neutrino mixing angle ๐œƒ 13 with the Double Chooz detector. JHEP, 10:086, 2014. [Erratum: JHEP 02, 074 (2015)]. [36] Xin Qian and Jen-Chieh Peng. Physics with Reactor Neutrinos. Rept. Prog. Phys., 82(3):036201, 2019. [37] P. Adamson et al. Measurement of Neutrino and Antineutrino Oscillations Using Beam and Atmospheric Data in MINOS. Phys. Rev. Lett., 110(25):251801, 2013. [38] P. Adamson et al. Electron neutrino and antineutrino appearance in the full MINOS data sample. Phys. Rev. Lett., 110(17):171801, 2013. [39] Mayly Sanchez. Nova results and prospects, June 2018. [40] Taichiro Koga. Measurement of neutrino interactions on water and search for electron anti-neutrino appearance in the T2K experiment. PhD thesis, Tokyo U.โ€ž 2018. [41] K. Abe et al. Atmospheric neutrino oscillation analysis with external constraints in Super- Kamiokande I-IV. Phys. Rev. D, 97(7):072001, 2018. [42] M. G. Aartsen et al. Determining neutrino oscillation parameters from atmospheric muon neutrino disappearance with three years of IceCube DeepCore data. Phys. Rev. D, 91(7):072004, 2015. [43] S. Euler. Observation of oscillations of atmospheric neutrinos with the IceCube Neutrino Observatory. PhD thesis, Aachen, 2014. Aachen, Techn. Hochsch., Diss., 2014. [44] Evgeny K. Akhmedov. Parametric resonance of neutrino oscillations and passage of solar and atmospheric neutrinos through the earth. Nucl. Phys. B, 538:25โ€“51, 1999. [45] T. K. Gaisser and M. Honda. Flux of atmospheric neutrinos. Ann. Rev. Nucl. Part. Sci., 52:153โ€“199, 2002. [46] M G Aartsen, R Abbasi, M Ackermann, J Adams, J A Aguilar, M Ahlers, M Ahrens, C Alispach, P Allison, N M Amin, and et al. Icecube-gen2: the window to the extreme universe. Journal of Physics G: Nuclear and Particle Physics, 48(6):060501, Apr 2021. [47] K.A. Olive. Review of particle physics. Chinese Physics C, 38(9):090001, aug 2014. [48] Anatoli Fedynitch, Ralph Engel, Thomas K. Gaisser, Felix Riehn, and Todor Stanev. Cal- culation of conventional and prompt lepton fluxes at very high energy. EPJ Web Conf., 99:08001, 2015. 106 [49] M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa. Atmospheric neutrino flux calculation using the nrlmsise-00 atmospheric model. Phys. Rev. D, 92:023004, Jul 2015. [50] M. Leuermann. Testing the Neutrino Mass Ordering with IceCube DeepCore. PhD thesis, RWTH Aachen U., 2018. [51] Adam M. Dziewonski and Don L. Anderson. Preliminary reference earth model. Physics of the Earth and Planetary Interiors, 25(4):297โ€“356, 1981. [52] S. Wren. Neutrino mass ordering studies with IceCube-DeepCore. PhD thesis, University of Manchester, 2018. [53] J. A. Formaggio and G. P. Zeller. From eV to EeV: Neutrino Cross Sections Across Energy Scales. Rev. Mod. Phys., 84:1307โ€“1341, 2012. [54] Ian J. R. Aitchison and Anthony J. G. Hey. Gauge Theories in Particle Physics: A Practical Introduction, Volume 2: Non-Abelian Gauge Theories : QCD and The Electroweak Theory, Fourth Edition. Taylor & Francis, 2013. [55] Raj Gandhi, Chris Quigg, Mary Hall Reno, and Ina Sarcevic. Ultrahigh-energy neutrino interactions. Astropart. Phys., 5:81โ€“110, 1996. [56] PDG, atomic and nuclear properties of materials for more than 350 materials. https: //pdg.lbl.gov/2021/AtomicNuclearProperties/MUE/muE_water_ice.pdf. Ac- cessed: 2022-03-11. [57] T.D.S. Stanislaus. Atomic Capture of Negative Muons in Oxides [microform]. Thesis: M. Sc. Thesis (M.Sc.)โ€“University of British Columbia, 1983. [58] Nimai C. Mukhopadhyay. NUCLEAR MUON CAPTURE. Nucl. Phys. A, 335:111โ€“135, 1980. [59] S. Agostinelli et al. GEANT4โ€“a simulation toolkit. Nucl. Instrum. Meth. A, 506:250โ€“303, 2003. [60] E. Fermi and E. Teller. The capture of negative mesotrons in matter. Phys. Rev., 72:399โ€“408, Sep 1947. [61] D. F. Measday. The nuclear physics of muon capture. Phys. Rept., 354:243โ€“409, 2001. [62] Shirvel Stanislaus, Frida Entezami, and David F. Measday. Atomic capture of muons and density. Nuclear Physics, 475:642โ€“656, 1987. [63] Nimai C. Mukhopadhyay. Nuclear muon capture. Physics Reports, 30(1):1โ€“144, 1977. [64] John David Jackson. Classical electrodynamics. Wiley, New York, NY, 3rd ed. edition, 1999. [65] J. Feintzeig. Searches for Point-like Sources of Astrophysical Neutrinos with the IceCube Neutrino Observatory. PhD thesis, The University of Wisconsin - Madison, January 2014. 107 [66] The IceCube Collaboration (internal website), refractive index of ice. https://wiki. icecube.wisc.edu/index.php/Refractive_index_of_ice. Accessed: 2021-08-23. [67] P Price. Role of group and phase velocity in high-energy neutrino observatories. Astropar- ticle Physics, 15(1):97โ€“100, Mar 2001. [68] M.G. Aartsen, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers, M. Ahrens, D. Altmann, K. Andeen, T. Anderson, I. Ansseau, and et al. The icecube neutrino observatory: instru- mentation and online systems. Journal of Instrumentation, 12(03):P03012โ€“P03012, Mar 2017. [69] The IceCube Collaboration, icecube internal gallery. https://gallery.icecube.wisc. edu/. Accessed: 2021-08-03. [70] M.G. Aartsen, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers, M. Ahrens, D. Altmann, K. Andeen, T. Anderson, I. Ansseau, and et al. Measurement of atmospheric tau neutrino appearance with icecube deepcore. Physical Review D, 99(3), Feb 2019. [71] R. Abbasi, Y. Abdou, T. Abu-Zayyad, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers, M.M. Allen, D. Altmann, K. Andeen, and et al. The design and performance of icecube deepcore. Astroparticle Physics, 35(10):615โ€“624, May 2012. [72] R. Abbasi, Y. Abdou, M. Ackermann, J. Adams, J.A. Aguilar, M. Ahlers, D. Altmann, K. Andeen, J. Auffenberg, X. Bai, and et al. Icetop: The surface component of icecube. Nu- clear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 700:188โ€“220, Feb 2013. [73] The IceCube Collaboration (internal website), coordinate system. https://wiki. icecube.wisc.edu/index.php/Coordinate_system. Accessed: 2022-03-18. [74] C. Kopper. Performance Studies for the KM3NeT Neutrino Telescope. PhD thesis, Friedrich- Alexander-Universitรคt Erlangen-Nรผrnberg, 2010. [75] M. Ackermann et al. Optical properties of deep glacial ice at the South Pole. J. Geophys. Res., 111(D13):D13203, 2006. [76] The IceCube Collaboration (internal website), ice flow. https://wiki.icecube.wisc. edu/index.php/Ice_flow. Accessed: 2022-03-22. [77] The IceCube Collaboration (internal website), ice model plots for public release. https://wiki.icecube.wisc.edu/index.php/Ice_model_plots_for_ public_release. Accessed: 2022-03-23. [78] P. Eller. The IceCube Collaboration (internal website), unified angular acceptance model. https://wiki.icecube.wisc.edu/index.php/Ice_models. Accessed: 2022-07-24. [79] Michael James Larson. Simulation and identification of non-Poissonian noise triggers in the IceCube neutrino detector. PhD thesis, Alabama U., 2013. 108 [80] The IceCube Collaboration, PMTResponseSimulator. https://docs.icecube.aq/ icetray/main/projects/DOMLauncher/PMTRes.html. Accessed: 2022-03-18. [81] Aya Ishihara. The IceCube Upgrade - Design and Science Goals. PoS, ICRC2019:1031, 2021. [82] D. Kraft. A Software Package for Sequential Quadratic Programming. Deutsche Forschungs- und Versuchsanstalt fรผr Luft- und Raumfahrt Kรถln: Forschungsbericht. Wiss. Berichtswesen d. DFVLR, 1988. [83] G. D. Barr, T. K. Gaisser, S. Robbins, and Todor Stanev. Uncertainties in Atmospheric Neutrino Fluxes. Phys. Rev. D, 74:094009, 2006. [84] Carlos A. Argรผelles, Austin Schneider, and Tianlu Yuan. A binned likelihood for stochastic models. JHEP, 06:030, 2019. [85] X. Qian, A. Tan, J.J. Ling, Y. Nakajima, and C. Zhang. The gaussian cls method for searches of new physics. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 827:63โ€“78, 2016. [86] Hamamatsu, data sheet. photomultiplier tube r7081-02 for icecube experiment. https://user-web.icecube.wisc.edu/~kitamura/NK/PMT/031112%20R7081-02% 20data%20sheet.pdf. Accessed: 2021-08-23. [87] R. Abbasi, Y. Abdou, T. Abu-Zayyad, J. Adams, J.A. Aguilar, M. Ahlers, K. Andeen, J. Auffenberg, X. Bai, M. Baker, and et al. Calibration and characterization of the icecube photomultiplier tube. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 618(1-3):139โ€“152, Jun 2010. [88] R. Abbasi, M. Ackermann, J. Adams, M. Ahlers, J. Ahrens, K. Andeen, J. Auffenberg, X. Bai, M. Baker, S.W. Barwick, and et al. The icecube data acquisition system: Signal capture, digitization, and timestamping. Nuclear Instruments and Methods in Physics Research Sec- tion A: Accelerators, Spectrometers, Detectors and Associated Equipment, 601(3):294โ€“316, Apr 2009. [89] Rasha Abbasi et al. Design and performance of the multi-PMT optical module for IceCube Upgrade. PoS, ICRC2021:1070, 2021. [90] M. G. Aartsen et al. IceCube-Gen2: the window to the extreme Universe. J. Phys. G, 48(6):060501, 2021. [91] Ryo Nagai and Aya Ishihara. Electronics Development for the New Photo-Detectors (PDOM and D-Egg) for IceCube-Upgrade. PoS, ICRC2019:966, 2020. [92] K.J. Ma, W.G. Kang, J.K. Ahn, S. Choi, Y. Choi, M.J. Hwang, J.S. Jang, E.J. Jeon, K.K. Joo, H.S. Kim, and et al. Time and amplitude of afterpulse measured with a large size photomul- tiplier tube. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 629(1):93โ€“100, Feb 2011. 109 [93] The IceCube Collaboration (internal website), afterpulse measurement. https://wiki. icecube.wisc.edu/index.php/Afterpulse_measurement. Accessed: 2022-03-18. [94] J Haser, F Kaether, C Langbrandtner, M Lindner, S Lucht, S Roth, M Schumann, A Stahl, A Stรผken, and C Wiebusch. Afterpulse measurements of r7081 photomultipliers for the double chooz experiment. Journal of Instrumentation, 8(04):P04029โ€“P04029, Apr 2013. [95] Hamamatsu, data sheet. photomultiplier tube r7081. https://www.hamamatsu.com/ resources/pdf/etd/LARGE_AREA_PMT_TPMH1376E.pdf. Accessed: 2022-01-04. [96] The IceCube Collaboration (internal website), afterpulse data. https://wiki.icecube. wisc.edu/index.php/Afterpulse_Data. Accessed: 2022-03-18. [97] D. Liu and C. Wendt. The IceCube Collaboration (internal website), Dahai Liuโ€™s mea- surements of large pulse response. https://user-web.icecube.wisc.edu/~chwendt/ afterpulse-notes/. Accessed: 2022-07-21. [98] A Rohatgi. Webplotdigitizer:version 4.5. https://automeris.io/WebPlotDigitizer, 2021. [99] The IceCube Collaboration (internal website), weights in nugen. https://docs.icecube. aq/icetray/main/projects/neutrino-generator/weighting.html. Accessed: 2022-06-03. [100] The IceCube Collaboration, weighting of genie-reader simulation. https://docs. icecube.aq/icetray/main/projects/genie-reader/event_weighting.html. Accessed: 2022-06-03. [101] The IceCube Collaboration (internal report), effective livetime for weighted cor- sika. https://internal-apps.icecube.wisc.edu/reports/data/icecube/2009/ 02/001/icecube_200902001_v2.pdf. Accessed: 2022-03-26. [102] Costas Andreopoulos, Christopher Barry, Steve Dytman, Hugh Gallagher, Tomasz Golan, Robert Hatcher, Gabriel Perdue, and Julia Yarba. The GENIE Neutrino Monte Carlo Gen- erator: Physics and User Manual. 10 2015. [103] The IceCube Collaboration, MuonGun. https://docs.icecube.aq/icetray/main/ projects/MuonGun/index.html. Accessed: 2022-07-21. [104] Leif Radel and Christopher Wiebusch. Calculation of the Cherenkov light yield from electromagnetic cascades in ice with Geant4. Astropart. Phys., 44:102โ€“113, 2013. [105] J.-H. Koehne, K. Frantzen, M. Schmitz, T. Fuchs, W. Rhode, D. Chirkin, and J. Becker Tjus. Proposal: A tool for propagation of charged leptons. Computer Physics Communications, 184(9):2070โ€“2090, 2013. [106] Leif Rรคdel and Christopher Wiebusch. Calculation of the cherenkov light yield from low energetic secondary particles accompanying high-energy muons in ice and water with geant4 simulations. Astroparticle Physics, 38:53โ€“67, 2012. 110 [107] The IceCube Collaboration (internal website), slc hit cleaning. https://wiki.icecube. wisc.edu/index.php/SLC_hit_cleaning. Accessed: 2022-03-26. [108] J. P. Yanez. Status of standard oscillation physics with IceCube DeepCore. 1468:012122, feb 2020. [109] The IceCube Collaboration (internal website), santa algorithm. https://wiki.icecube. wisc.edu/index.php/SANTA_algorithm. Accessed: 2022-03-27. [110] J. Ahrens, X. Bai, R. Bay, S.W. Barwick, T. Becka, J.K. Becker, K.-H. Becker, E. Bernardini, D. Bertrand, A. Biron, and et al. Muon track reconstruction and data selection techniques in amanda. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 524(1-3):169โ€“194, May 2004. [111] J.A. Aguilar et al. A fast algorithm for muon track reconstruction and its application to the antares neutrino telescope. Astroparticle Physics, 34(9):652โ€“662, 2011. [112] Jan-Patrick HรผlรŸ and Christopher Wiebusch. Search for neutrinos from the direction of the galactic center with the icecube neutrino telescope. 2010. [113] The IceCube Collaboration, first guess of a start/stop point. https://docs.icecube.aq/ icetray/main/projects/finiteReco/I3StartStopPoint.html. Accessed: 2022- 07-21. [114] J. Lundberg, P. Mioฤinoviฤ‡, K. Woschnagg, T. Burgess, J. Adams, S. Hundertmark, P. Desiati, and P. Niessen. Light tracking through ice and waterโ€”scattering and absorption in heteroge- neous media with photonics. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 581(3):619โ€“631, Nov 2007. [115] Nathan Whitehorn, Jakob van Santen, and Sven Lafebre. Penalized splines for smooth rep- resentation of high-dimensional monte carlo datasets. Computer Physics Communications, 184(9):2214โ€“2220, Sep 2013. 111 APPENDIX A DATA ACQUISITION AND ONLINE PROCESSING IN ICECUBE In order to identify single neutrino interactions within the IceCube detector, from here onwards identified as neutrino events, light recorded by the optical modules has to be analyzed to determine time windows of interest, where physics-like light deposition patterns can be identified. This identification process has two distinct stages: โ€ข Online stage: raw information produced by the data acquisition system within the DOMs is collected at the ICL computer hub, time ordered and processed through a series of triggers to determine tentative physics event data blocks. Furthermore, a series of simple filters are run on each event, identifying what physics analyses they could be useful for. This stage takes place at the SPS. โ€ข Offline stage: in the Northern Hemisphere, computationally expensive filters are run over the events identified during the online processing. In addition, within the framework of each specific physics analysis, a series of increasingly complex cuts and reconstruction algorithms to identify neutrino event samples with specific characteristics are run. In this appendix both the data acquisition system within the IceCube DOMs and the online candidate event identification and processing at the south pole are presented. On the other hand, the specific offline event selection criteria for the neutrino sample used in the NMO analysis with neutrino/antineutrino separation is presented in Appendix E. A.1 The IceCube DOM: in-ice data acquisition The fundamental detection component in IceCube is the DOM (see figure A.1). These optical modules detect Cherenkov light produced by charged particles in the ice, digitize the signal and send the digitized waveforms, together with a time stamp, towards the ICL [68]. There are a total of 5160 DOMs within the Antarctic ice and another 324 DOMs within the IceTop surface array. The main components of these optical modules are: 112 โ€ข A โˆผ13 cm radius R7081-02 PMT [86], facing downwards towards the northern sky (R7081- 02MOD for HQE DOMs); that is, with its sensitive surface located at the lower hemisphere of the DOM [68]. These PMTs have a peak quantum efficiency of โˆผ 25% at 390 nm (or โˆผ34% in HQE DOMs) and are run at a 107 gain, a value chosen to measure individual photon pulses with amplitudes much higher than the electronic noise level and the precision of the digitizers (โˆผ8 mV for a single photo-electron, โˆผ0.1 mV for both noise components) [87]. โ€ข A glass sphere, protecting the PMT and electronics from deep ice pressure. The sphere has a โˆผ 16.5 cm radius and 1.27 cm thickness [68]. โ€ข A Silicone gel filling in the lower hemisphere of the DOM, providing both optical coupling between the PMT surface and the DOM glass sphere and structural support to the PMT overall [68]. โ€ข A mu-metal grid surrounding the PMT to provide shielding against external magnetic fields (ambient flux density of 550 mG at South Pole) [68]. โ€ข A Main Board, the data acquisition (DAQ) circuit within the DOM. The main board supplies power to all components and in addition, it contains all elements involved in waveform digi- tization and processing, manages communication with neighboring DOMs and the IceCube Lab and allows calibration of DOM components [68]. โ€ข A Flasher Board, containing 12 LEDs that provide an in-situ light source for calibration purposes [68]. โ€ข A Cable system, providing communication with neighboring DOMs and the ICL. This cable system, consisting of three wire pairs, enters the DOM through a 16.3 mm hole in the glass and ultimately couples to the main string cable. The first wire pair carries power and allows bidirectional communication with a computer in the IceCube Laboratory. The other two communicate to neighboring DOMs and are used to find local coincidence (LC) of light signals, aiding on the identification of physical events [68]. 113 Figure A.1: An overall schematic of the IceCube DOM showing its main components [69]. Signal processing in the DOM starts at the PMT anode. Voltage signals are compared to a threshold voltage equivalent to a 0.25 PE signal, typically called the discriminator threshold. Signals with voltages above the threshold trigger a data taking instance or launch, in which the signal waveform is digitized using both Analog Transient Waveform Digitizers (ATWDs) and a fast analog-to-digital converter (fADC). An ATWD chip contains three distinct channels with amplifier gains of x16, x2 and x0.5 for the readout and digitization of waveforms with a wide range of amplitudes, as would be expected from neutrino signals ranging over several orders of magnitude in energy, and a fourth channel for calibration purposes. Information in a low gain channel is digitized only when the channels with higher amplification are at a level of 75% saturation or above [68]. There are two ATWD chips per DOM. This allows recording signals even when one of the chips is in the process of digitization, which takes about 29 ๐œ‡s per readout ATWD channel. A DOM will only have significant dead-time after two successive launches [88]. 114 The ATWDs collect up to a total of 128 voltage samples per launch, recording at a rate of 300 Msps. This is equivalent to 427 ns of data and covers the time range expected from physical signals originated within tens of meters from a DOM [68]. Events further away from the optical modules will nevertheless produce lower amplitude signals at a delayed time window, as a consequence of light scattering in the ice. The fADC is responsible for recording these signals, collecting up to 256 samples per launch and sampling continually at 25 ns intervals [68]. The full digital record collected during a launch is typically referred to as a "hit" and is the basic unit datum produced by the IceCube DOMs. Although both ATWDs and fADC start recording data in a launch, fully digitized waveforms are only saved and subsequently sent to the IceCube laboratory if the LC condition is satisfied. When a launch starts in any DOM, a LC signal is sent to its neighboring DOMs and the DOMs immediately above and below those neighbors [68]. The LC condition requires a launched DOM to receive a LC signal sent within 1 ๐œ‡s of its launch time, accounting for the travel time of this signal (takes up to 1325 ns to reach all relevant DOMs). When the LC condition fails, ATWD digitization is aborted and only a time stamp is saved, together with a course charge stamp from the highest recording in the fADC since the start of the launch (together with the fADC samples before and after it). This is known as a soft local coincidence (SLC) hit. Otherwise fully digitized information including all relevant ATWD channels, 6.4 ๐œ‡s of fADC data and both time and course charge stamps are recorded and eventually transferred to the IceCube laboratory; this is called a hard local coincidence (HLC) hit [88]. Individual pulses retrieved from these hits are then labelled HLC or SLC as well, an important fact when discussing neutrino event selection for the NMO analysis (see Appendix E). A.2 Online processing The online processing in the ICL includes data acquisition, identification of select time windows that are likely to contain physics events, data storage and transference of data via satellite to the northern hemisphere. At the SPS, there is one computer, or DOMHub, assigned to each in-ice string (and 11 additional computers assigned to IceTop). DOM calibration and monitoring data, together with pulse data over 115 the full string, is read each second at the DOMHubs and processed by the so called "StringHub" Software. StringHub time orders the pulse stream and generates "condensed" HLC launches, containing only the time and a DOM ID. The condensed launches are then transferred to appropriate trigger algorithms, that search for time windows of interest [68]. There are several triggers within the online data acquisition system. They mostly search for clusters of light within the detector, since these are typical of physical events. Uncorrelated launches in the detector, on the other hand, are mainly related to dark noise. A comprehensive list with the main IceCube triggers is presented next: โ€ข The Simple Multiplicity Trigger (SMT): SMTN triggers requires a total of N HLC hits within a given time window. There are several flavors for this Trigger. The most common are the SMT8, applied to in-ice DOMs and with a time window of 5 ๐œ‡s and the SMT3, applied over DeepCore DOMs, with a time window of 2.5 ๐œ‡s. Once the SMTN condition is satisfied, the trigger window is extended until no HLC hits are found for a time period as long as the original trigger time window [68]. โ€ข Volume Trigger: Developed to trigger on localized light depositions with a low number of hits, typical of low energy events. The most common version of the trigger searches for 4 in-ice HLC hits within a time window of 1 ๐œ‡s and inside a cylinder of radius 175 m and height 75 m centered at each hit DOM (at each DOM containing a hit) [68]. โ€ข String Trigger: Developed to trigger on low energy muons traveling vertically. It identifies each possible set of 7 adjacent DOMs within a string and for each set, requires 5 HLC hits within a time window of 1.5 ๐œ‡s [68]. โ€ข The Slow Particle (SLOP) trigger: The trigger was developed to search for light deposition patterns expected from slow particles crossing the detector (i.e. magnetic monopoles), at speeds lower than 0.01c. Characteristic time windows for this trigger are in the order of hundreds of ๐œ‡s [68]. 116 Figure A.2: An example of a global trigger built from several individual trigger readout windows [68]. An event is formed from the collection of all HLC and SLC hits within the merged time window, together with all the trigger data used to define the global trigger. Trigger windows that are long enough to contain several independent physics events (i.e. the SLOP trigger) are typically re-split to make them useful to most physics analyses. Each trigger has a time window defined around it (typically 10 ๐œ‡s), to cover scattered light and early physical hits. The windows are asymmetric for triggers related to in-ice DOMs, but symmetric for IceTop triggers. Since a physical event can contain several triggers, all overlapping triggers are merged into a single one, called the global trigger (see figure A.2). The global trigger is usually re-split for readout windows long enough to contain several independent physics events, which is possible for some special triggers like the SLOP trigger [68]. Once a global trigger has been identified, a special software component, called the Event Builder, retrieves all hit data within its readout time window (both HLC and SLC), together with trigger information, to form an individual event. The next stage of the online processing pipeline, the Processing and Filtering system (PnF) [68], reduces amount of data produced by the Event Builder to a level where it can be transferred via satellite (bandwidth โˆผ 100 GB/day) to the Northern hemisphere. Events that do not satisfy the PnF requirements are still saved in local memory disks that are transferred physically to the northern hemisphere once per year [68]. The PnF processing starts with the conversion of hit digitized waveforms into individual photon signals (reconstructed pulses), characterized by a charge and time. The process starts by trans- forming counts read by the digitizer to voltage. Times are also corrected for photon transit times 117 Figure A.3: Schematic of the online processing at the SPS [68]. The main steps indicated in this work are: data acquisition and basic processing by the StringHubs, identification of trigger windows, event building, the PnF and data movement to the Northern hemisphere via satellite. Additional paths in the processing, not accounted for in this description, are irrelevant for the work presented in this thesis. through the PMT dynode chain. Pulse profiles corresponding to the known responses of each DOM to individual photons are then fit to the digitized data to unfold charge and time. IceCube events therefore contain a list of all pulses identified during unfolding, each including charge, time and width information, as well as LC condition, digitizer (ATWD and/or FADC) and DOM ID [68]. Afterwards, a set of filter algorithms are run on each event, each designed to identify events that match a specific physics signal. Basic reconstructions are run on the events to aid their filtering, obtaining first estimates on energy, direction, interaction vertex, etc, by matching the light depositions in each event to those expected from track or shower hypotheses. There is a vast list of 118 filters within the online processing, targeting a wide range of physics analyses [68]. Finally, all events satisfying at least one filter condition are listed for transfer via satellite to the Northern hemisphere, these events include trigger and filter result information, basic reconstruc- tions, calibration data and hit information. A overview of the online processing presented in this section is depicted in figure A.3. 119 APPENDIX B THE ICECUBE UPGRADE The IceCube Upgrade is an upcoming improvement to the detector that consists of the addition of seven densely instrumented strings within the DeepCore sub-array. Several analyses currently performed using the IceCube detector will benefit from it, including neutrino astrophysics and neutrino oscillation analyses [81]. Upgrade strings will be deployed as shown in figure B.1, to form a grid with average inter-string distances of 20 m [81]. Each string will contain about 700 optical sensors, organized in two distinct regions: a physics region, located at depths between 2150 m and 2425 m, where the ice is cleanest, and a calibration region, located above and below the physics region and used for muon veto, calibration and testing of the Upgrade optical modules [81]. Inter optical-module distances in the physics and calibration region will be 3 m and โˆผ 25 m respectively [81]. In addition to the increased optical module density in the physics region of the Upgrade array, the new strings will also showcase three newly developed optical modules that represent an improvement over current IceCube DOMs, multi-PMT digital optical modules (mDOMs), Dual optical sensors in an Ellipsoid Glass for Gen2 (D-Eggs) and PDOMs [81]. The mDOMs will feature a total of 24 PMTs (โˆผ 3 inch diameter), arranged within a 3D printed support structure [89]. A schematic of an mDOM can be found in figure B.2. The inclusion of several PMTs pointing outwards will provide mDOMs with a larger photocathode area and a homogeneous angular acceptance curve [89], which is an improvement over traditional DOMs, where photons incoming from above were severely suppressed. In addition, mDOMs will open the possibility to recover photon arrival angles using a single optical module, which will enhance the detectorโ€™s angular resolution [89]. D-Eggs, on the other hand, will contain only two PMTs, pointing upwards and downwards [90], as shown in figure B.2. They have an ellipsoidal shape and offer better angular acceptance and photocathode area than traditional DOMs, but lower than that of mDOMs [90]. Lastly, PDOMs are simply improved DeepCore DOMs with updated electronics [81]. These 120 Figure B.1: Schematic of the IceCube detector that includes the Upgrade strings [81]. The figure shows the depth range for the physics region of the Upgrade and the geometry that was chosen for the instrumentation of the seven new strings (calibration regions have yet to be finalized). The Upgrade will also showcase two main new optical module types, as illustrated in the lower right diagram. Figure B.2: Showing an mDOM (left) and a D-Egg (right), the two main optical modules to be used in the IceCube Upgrade [90]. 121 Figure B.3: Improvement to low energy neutrino event rates in the IceCube Upgrade for an example neutrino selection criteria looking for fully contained events [81]. The inner (outer) fiducial volume is defined as a cylindrical region of height 275 m and radius 50 m (145 m), centered at the center of the DeepCore/Upgrade array. optical modules will mostly be found in the calibration region of the Upgrade, whereas mDOMs and D-Eggs will be the main optical modules deployed in the Upgrade strings and will be found in all regions [81]. Similarly to regular DOMs, the new modules will be enclosed within a glass pressure vessel and will include Silicone gel in the region between the photocathodes (and support structure for mDOMs) and the pressure vessel to improve optical coupling [90]. They will also contain their own DAQ and will be able to digitize continuously, removing dead-time effects [91, 89]. Amongst the upgrade modules, several calibration devices will be included, which will be useful for improving understanding on the properties of the Antarctic Ice [81]. There are two main reason why the Upgrade will offer improvements to future oscillations analyses. The first reason is that it will offer an improvement to low energy neutrino detection (since the detected photons per event will increase) and the probability for these neutrinos to pass simple selection criteria. To select neutrino events at low energies we need enough light to pass simple event quality selection cuts and to properly separate them from noise triggered events using directionality criteria. As shown in figure B.3, there is an observed increase in the rate of low energy neutrinos in the Upgrade, when compared to DeepCore. The second reason is that reconstruction algorithms will see a considerable improvement in resolution with the Upgrade, given the increased amount of photons [81] and (potentially) the 122 Figure B.4: Improvement to reconstruction resolutions for low energy neutrino events with the IceCube Upgrade [81]. directionality information per module. This is exemplified in figure B.4 The IceCube upgrade will have an important role in the future of the analysis presented in this work since it will incorporate an array design that will enhance the detectability of the muon decay signature photons while significantly reducing the impact of internal optical module backgrounds, as discussed in section 4.2. 123 APPENDIX C AFTERPULSES When a photoelectron travels through a PMT dynode chain, it can ionize residual gas in the space between the first dynode and the photocathode. The resulting ions will have a positive charge, and as a result, will travel back and hit the photocathode surface, releasing secondary photoelectrons and producing a secondary delayed signal, referred to as an afterpulse. These pulses can be found in a time window ranging between 300 to about 16000 ns after the time of the primary signal. Their large time delays are the result of the slower transit time of ions through the dynode chain, when compared to that of electrons, given their higher mass. As described in [92], the transit time of an ion through the space between the first dynode and the photocathode will be given by the following relation, assuming a linear voltage profile: โˆš๏ธ„ 2๐‘š๐ฟ (๐ฟ โˆ’ ๐‘ ) ๐‘ก= , (C.1) ๐‘ž๐‘‰0 where ๐‘‰0 and L are respectively the voltage difference and distance between photocathode and first dynode, q and m are the ion charge and mass, and s is the distance between the ion production point and the photocathode. As a result, the afterpulse profile within a PMT will include several time components, each corresponding to a different ion (element) within the residual gas in the dynode chain. C.1 Afterpulse model in the IceCube simulation chain In the past, there has has been a limited interest in characterizing the afterpulse profile in IceCube data, mainly since most physics analyses do not require pulses in the time window beyond a couple hundred ns and since the total charge contribution from afterpulses is significantly smaller than that from primary pulses, given their low generation probability [87]. As a result, the main afterpulse study for IceCube PMTs, performed in 2007 and described in [87], was carried out with only a few PMTs under laboratory conditions (room temperature) and only the results from a single PMT were introduced into the standard IceCube simulation chain. Although a more recent afterpulse measurement was performed by the collaboration in 2011 [93], it is generally 124 not included in the main simulation chain and does not introduce a noticeable improvement in the afterpulse modeling, measuring afterpulse rates for a small sample of PMTs, again under laboratory conditions. Afterpulse simulation has remained mainly untouched ever since and as a result, IceCube simulation currently lacks a proper characterization of afterpulse rate variations across the detector. Furthermore, more recent afterpulse studies, found in [92, 94], have shown different charge distributions for the several afterpulse components in the PMT model R7081 [95] (this model should have a similar afterpulse profile as the modified model for IceCube, since they differ mainly by the operation temperature range). As an example, figure C.1a shows the afterpulse profile for single photoelectron (SPE) and multi-photoelectron (MPE) afterpulses, from the measurements performed in [94]. MPE afterpulse components were found around the 600, 2000 and 8000 ns time windows. The current afterpulse model in IceCube does not account for this, instead all afterpulse components were assumed to produce charge according to a SPE distribution, with the exception of early afterpulses (around 600 ns) that were identified as multi-electron components with an average charge of โˆผ 13 PE, in the 2007 measurement, given the high charge peak that is found in this time range. At the time, it was determined that the charge profile for these afterpulses could be modelled using a Fisher-Tippett distribution: 1 ๐‘“ (๐‘ž|๐œ‡, ๐œŽ) = exp(โˆ’๐‘ ) exp(โˆ’ exp(โˆ’๐‘ )), (C.2) ๐œŽ where ๐œ‡, ๐œŽ are parameters of the distribution, q is charge and ๐‘  = (๐‘ž โˆ’ ๐œ‡)/๐œŽ. The charge fit results to afterpulse data between 500 and 800 ns were: ๐œ‡ = 13.31 and ๐œŽ = โˆ’3.386. To derive the current afterpulse model in IceCube, reviewed in [96, 87], bright LED pulses of varying intensity were fired towards the PMT surface and the resulting PMT waveforms, in a time window between 300 ns 16.8 ๐œ‡๐‘  after the primary peak from the LED signal, were measured. These waveforms were generated from the added afterpulse signals from all photoelectrons within the original pulse. Eleven distinct peaks in the afterpulse profile were then identified by fitting 125 IC AP components 103 All APs [0,2] PE [2,10] PE [10,25] PE 102 events 101 100 0 2 4 6 8 10 12 t [ s] (a) (b) Figure C.1: The figures above show a comparison between the afterpulse profiles reported in [94] for R7081 PMTs and those obtained with the MPE afterpulse model introduced into the IceCube simulation in this work. The profiles represent the cumulative distributions of the PMT response to SPE primary signals, over a large sample of SPE primary pulses. Each colored curve corresponds to afterpulses with different charge (Q) ranges. Green for Q < 2 PE, blue for 2 < Q < 10 PE, red for 10 < Q < 25 PE and black for the full profile. a) Distribution found by the Double Chooz collaboration in [94] using SPE LED pulses. The exponential profile at shorter times was identified as afterglowing in the LED and was not present in the afterpulse measurements by IceCube. b) Distribution obtained from afterpulse simulation using several simulated SPE primary pulses and the IceCube MPE afterpulse model produce during this work. The exponential component from afterglowing is missing. Vertical lines correspond to the 11 afterpulse components in IceCube simulation. the waveforms with a set of Gaussian components. Table C.1 shows the Gaussian fit that was performed for afterpulses produced by a 1.9 ร— 106 PE LED pulse [96]. On the other hand, afterpulse rate was calculated using the ratio between primary pulse charge and derived afterpulse charge. Within this process, it was noticed that the current readout at times larger than 11.3 ๐œ‡๐‘  after the original LED pulse corresponded to a constant baseline response from the PMT in the absence of light. As a result, the total afterpulse integrated charge between 300 ns and 11.3 ๐œ‡๐‘  was found, and then the integrated charge at later times was subtracted from it. This procedure was repeated for several different LED pulses, and then figure C.2, showing the relation between primary pulse charge and integrated afterpulse charge was generated. A linear fit to the plot in figure C.2 demonstrated that for charges below the saturation charge 126 Component Amplitude [mA] Peak time [ns] Peak width [ns] Average Photoelectrons 1 -10 500 200 1 2 -23 540 20 13 3 -45 660 20 13 4 -6.5 1100 100 1 5 -4.5 1300 200 1 6 -8.5 1650 225 1 7 -8.95 2075 300 1 8 -5.15 2850 500 1 9 -3 4750 700 1 10 -5 6150.3 400 1 11 -17.44 7833.4 944.6 1 Table C.1: Table indicating all afterpulse components obtained when fitting a multiple Gaussian model to the waveform obtained using a bright LED pulse of 1.9 ร— 106 PE. This is the model currently available in IceCube simulation (adapted from [96]). Figure C.2: Showing the method used in the 2007 afterpulse measurement by the IceCube collab- oration [97] to infer afterpulse rate. 127 of the PMT, there was a linear relation between primary pulse charge and afterpulse charge, that corresponded to about 0.06 PEs of afterpulsing per primary SPE. C.2 MPE Afterpulse model developed for future decay signal extraction studies The NMO analysis performed in this work uses ๐œˆ ๐œ‡ /๐œˆยฏ ๐œ‡ discrimination based on decay signals and inelasticity. Decay signals are currently not extracted at a significant level due, partially, to their low level of separation from afterpulse signals. MPE afterpulses have the potential for improving the extraction of decays signals by introducing differences in the charge distributions of afterpulses and decay pulses. A MPE model for afterpulsing compatible with the IceCube simulation chain was developed for the analysis, but the effect of MPE afterpulsing was ultimately not explored. A future study is required to test the effect of MPE afterpulses in decay signal extraction. Below, the modeling process is outlined. Given the lack of a measurement on the charge profile of afterpulses performed by IceCube, an intermediate approach had to be taken, where the charge distributions for MPE afterpulse com- ponents observed in [92, 94] were taken from literature and matched to the afterpulse components measured by IceCube in 2007. The procedure started with extracting charge data from the plots presented in [92] for the MPE afterpulses around 2 and 8 ๐œ‡๐‘ , using the automeris web plot digitizer [98]. The data points were then fitted using a Fisher-Tippett distribution (equation C.2), to mimic the method used in the 2007 measurement when modeling MPE afterpulses around 600 ns. The results are shown in figure C.3. Although the Fisher-Tippett distribution were not the most optimal choice for fitting these charge distributions, the lack of a dedicated IceCube charge measurement meant that any description of MPE afterpulses would be more correct than the current model in simulation, but would always be a rough approximation. As a result, the fits were deemed sufficient. After fitting the charge distributions, the next step was to match the afterpulse component in simulation to the right distributions. The distribution in figure C.3a was matched to all afterpulse components between 1000 ns and 3000 ns, whereas the one in figure C.3b was only matched to the latest afterpulse component in IceCube simulation (at โˆผ 8 ๐œ‡s ). The afterpulsing probabilities 128 (a) (b) Figure C.3: The figures above show fits to the charge distributions for MPE afterpulses presented in [92] for a) MPE afterpulses around 2๐œ‡ s and b) MPE afterpulses around 8๐œ‡ s, obtained for a PMT high-voltage of 1450 V. Fits were performed using a simple least squares to a Fisher Tippett distribution. Obtained p-values are low, implying this is not an optimal modeling for the underlying charge distributions. Nevertheless, the goal of the MPE afterpulse study in this work was restricted to investigate if the existence of MPE afterpulses could aid in the extraction of decay signals. A precise measurement was beyond the scope of the work presented in this thesis, given the lack of a proper measurement with IceCube DOMs and as a result, the approximated charge distributions were deemed sufficient. per primary PE for MPE afterpulses were then modified in simulation, in order to retain the overall amplitude ratios between afterpulse components in figure C.1 and ultimately an average charge ratio of 0.06 afterpulse PEs per primary SPE. Figure C.1b shows the resulting aggregation of all afterpulses produced from a set of 1.5 ร— 107 SPE simulated primary pulses, using the modified simulation scripts. The simulation also included an internal loop that allowed generation of afterpulses from prior afterpulses. The results shows similarities with those found in [94], although a perfect fit was not possible for several reasons; the charge distributions used corresponded to a different article, mainly due to no access to the raw data used in [94], the experimental results in C.1a include an exponential profile at lower charges from afterglow in the LED that was used to produce the primary pulses, which is not present in IceCube simulation and finally, our charge model itself is only approximate. The afterpulse spectrum shown 129 on figure C.1b corresponds to the simulated output at the PMT anode and is further modified once simulation of DOM readout and pulse reconstruction is included. Although an improved afterpulse model was generated for the purpose of this work, the potential improvement in decay signal extraction provided by the presence of MPE afterpulse components has not yet been studied. 130 APPENDIX D DEEPCORE NMO ANALYSIS - EVENT SIMULATIONS AND WEIGHTING The following appendix summarizes the steps to generate the MC simulation used in this work. This sample was developed independently by the IceCube collaboration and tested on a simple neutrino oscillation analysis where both ๐œƒ 23 and ฮ”๐‘š 32 2 were fit assuming NO [2]. Note that the analysis sample used on this work only has a sub-percent contribution from atmospheric muons (see Appendix E). As a result, it was deemed unnecessary to reprocess these events to include muon capture physics as the corrections to the fit oscillation parameters from separation of atmospheric muons and antimuons should be negligible, especially since depositions of decay light in the detector are expected to be extremely rare (see section 4.1). In addition, this appendix introduces the weighting scheme used in oscillation analyses with IceCube simulation. This scheme allows simulation production to be independent of atmospheric neutrino fluxes, neutrino interaction models, neutrino oscillation parameters, etc. Instead of generating multiple MC sets for each model available, expectations according to different models can be obtained from a unique MC set by applying appropriate weights. D.1 Particle generation of atmospheric neutrino events Particle generation is the first step in the simulation chain for atmospheric muon and atmospheric neutrino events. Atmospheric neutrinos were injected close to and within the detector, without a dedicated simulation of the cosmic ray air showers. Specifically, each individual neutrino was injected at the end cap of a cylinder that could be rotated around an specific pivot, located roughly at the center of the DeepCore fiducial volume. The cylinder axis was chosen to match the direction of propagation of the incoming neutrinos, and the neutrinos were forced to interact within the cylinder volume (uniformly distributed). The direction of propagation of each neutrino event was drawn from a flat solid angle distribution; which in turn rotated the cylinders about their common pivot, allowing for simulation of events with all zenith and azimuth directions. The lengths and radii of the injection cylinders were energy dependent, chosen to efficiently sample all events that could generate light signals within the detector including those where the neutrino interaction itself didnโ€™t 131 happen inside the detector. Events were generated following a ๐ธ โˆ’2 spectrum and a flat directional spectrum. To convert the detection rates of the events generated using this method to those expected in real data each event was assigned a weight that accounted for its interaction probability and for a proper atmospheric neutrino flux spectrum. The methodology used to define these weights is summarized as: โ€ข Each neutrino was injected arriving from some zenith (๐œƒ) and azimuth (๐œ™) direction with an energy drawn from an ๐ธ โˆ’2 spectrum and was forced to interact at some position ๐‘ฅยฎ๐‘–๐‘›๐‘ก inside its injection cylinder. โ€ข The probability to inject an event with energy E was given by: ๐ธ โˆ’2 ๐‘ƒ(๐ธ) = โˆซ ๐ธ , (D.1) ๐‘š๐‘Ž๐‘ฅ ๐ธ ๐‘š๐‘–๐‘› ๐ธ โˆ’2 ๐‘‘๐ธ where ๐ธ ๐‘š๐‘Ž๐‘ฅ and ๐ธ ๐‘š๐‘–๐‘› were the maximum and minimum injected neutrino energies. โ€ข Given the direction of propagation of a neutrino, its path length from the end-cap of the cylinder to its interaction point can be used to calculate the probability that the neutrino interacted at that position, given by:  ๐‘ƒ๐‘–๐‘›๐‘ก (๐ธ) โ‰ˆ 1 โˆ’ exp โˆ’๐œŽ๐‘ก๐‘œ๐‘ก (๐ธ) โˆ— ๐ฟ ๐‘‘ /๐‘€ ๐‘ , (D.2) โˆซ where ๐‘€ ๐‘ is the proton mass, ๐ฟ ๐‘‘ is total column depth (๐ฟ ๐‘‘ = ๐œŒ(๐‘ฅ)๐‘‘๐‘ฅ along the path length), and ๐œŽ๐‘ก๐‘œ๐‘ก (๐ธ) is the energy dependent total neutrino cross section for neutrino- nucleon interactions [99]. This probability would typically then be multiplied by the survival probability of the neutrinos from production at the atmosphere to the cylinder end-cap, although this probability is essentially 1, and therefore not accounted for in the interaction weights. โ€ข The weights for each event were defined as follows [100]: " # ๐‘ƒ (๐ฟ ๐‘–๐‘›๐‘ก ๐‘‘ , ๐ธ, ๐œˆ๐‘– ) ๐‘Š๐‘–๐‘ ๐‘–๐‘š = ๐ด๐‘”๐‘’๐‘› ยท ฮฉ๐‘”๐‘’๐‘› ยท , (D.3) ๐‘ƒ(๐ธ) ยท ๐‘๐‘ก๐‘œ๐‘ก (๐œˆ๐‘– ) 132 where ๐ด๐‘”๐‘’๐‘› is the area of the generation cylinder end-cap, ฮฉ๐‘”๐‘’๐‘› is the generation solid angle, ๐œˆ๐‘– is the neutrino type interacting in the event, and ๐‘๐‘ก๐‘œ๐‘ก (๐œˆ๐‘– ) is the total number of events generated of the neutrino type ๐œˆ๐‘– . โ€ข These weights could be used after the fact (once simulations were finalized and processed through a neutrino selection criteria) to define flux model dependent weights of the form: ๐‘‘ฮฆ๐œˆ๐‘– (๐ธ, ๐œƒ, ๐œ™) ๐‘’๐‘ฅ ๐‘ ๐‘ค๐‘– = ๐‘Š๐‘–๐‘ ๐‘–๐‘š ยท , (D.4) ๐‘‘๐ธ ๐‘‘ฮฆ(๐ธ,๐œƒ,๐œ™) ๐‘’๐‘ฅ ๐‘ where ๐‘‘๐ธ is the expected neutrino flux per unit time (derived using an atmospheric neutrino flux model). The weights in equation D.4 are defined such that the sum of the weights over all events of type ๐œˆ๐‘– , generated with parameters E, ๐œƒ, ๐œ™ and interacting at a given position ๐‘ฅยฎ๐‘–๐‘›๐‘ก returns the expected rates of event interactions if the simulation had been performed using realistic atmospheric neutrino fluxes and neutrino interaction models. The weights may be rewritten as: expected flux per unit time ๐‘ค๐‘– = ๐‘ƒ๐‘–๐‘›๐‘ก (๐ธ, ๐œƒ, ๐œ™, ๐œˆ๐‘– , ๐‘ฅยฎ๐‘–๐‘›๐‘ก ) ยท , (D.5) generated flux ๐ธ,๐œƒ,๐œ™,๐œˆ๐‘– ,ยฎ ๐‘ฅ๐‘–๐‘›๐‘ก since ๐ฟ ๐‘‘ can be inferred from ๐œƒ, ๐œ™ and ๐‘ฅยฎ๐‘–๐‘›๐‘ก . Weights could therefore be used to trans- form distributions drawn from raw counts of simulated events into real rate (in ๐‘ โˆ’1 units) expectations in the detector given current best models of atmospheric neutrino fluxes and neutrino-nucleon interactions. To obtain expected event counts the rates calculated could be further scaled by the desired detector livetime (7.5 years for this analysis). โ€ข The weights in equation D.4 are valid per neutrino type, where the cross sections and expected fluxes are different for each type. As a result, weights were defined differently for neutrinos and antineutrinos. In the simulation, these particle types were produced in a 0.7:0.3 ratio. As a result, the weights are redefined as follows: ๏ฃด ๐‘‘ฮฆ๐œˆ (๐ธ,๐œƒ,๐œ™)๐‘’๐‘ฅ ๐‘ ๐‘ƒ๐‘–๐‘›๐‘ก,๐œˆ (๐ฟ ๐‘‘ ,๐ธ) ๏ฃฑ ยท ๐ด๐‘”๐‘’๐‘› ยท ฮฉ๐‘”๐‘’๐‘› ร— ๐‘ƒ(๐ธ)ยท0.7ยท๐‘ , for neutrinos ๏ฃด ๏ฃด ๐‘‘๐ธ ๐‘ก๐‘œ๐‘ก (๐œˆ) ๏ฃฒ ๏ฃด ๐‘ค๐‘– = , (D.6) ๐‘ƒ๐‘–๐‘›๐‘ก, ๐œˆยฏ (๐ฟ ๐‘‘ ,๐ธ) ๏ฃด ๐‘‘ฮฆ๐œˆยฏ (๐ธ,๐œƒ,๐œ™)๐‘’๐‘ฅ ๐‘ ยท ๐ด๐‘”๐‘’๐‘› ยท ฮฉ๐‘”๐‘’๐‘› ร— ๏ฃด ยฏ , for antineutrinos ๏ฃด ๐‘‘๐ธ ๐‘ƒ(๐ธ)ยท0.3ยท๐‘ ๐‘ก๐‘œ๐‘ก ( ๐œˆ) ๏ฃด ๏ฃณ where weights are also to be defined independently for ๐œˆ๐‘’ , ๐œˆ ๐œ‡ and ๐œˆ๐œ events. 133 The statistical uncertainties for counts predicted using unweighted events follow Poisson statistics (for large counts). Those predicted with weighted events are calculated using the sum of squared weights ( ๐‘ค๐‘–2 ). The 0.7:0.3 generation ratio used in the simulation ensured comparable magnitudes ร for neutrino and antineutrino weights. The effective livetimes (number of years of data that would be required to obtain the same statistical fluctuations in the number of events as those obtained with the weighted MC) of weighted simulation are given by [101]: ร ๐‘ค๐‘– ๐‘‡๐‘’ ๐‘“ ๐‘“ = ร , (D.7) ๐‘ค๐‘–2 and thus similar magnitudes for neutrino and antineutrino weights ensured the MC sample had comparable effective livetimes for the two particle types. Cross-sections for all neutrino interactions and flavors were calculated using a neutrino gener- ation software package called GENIE [102], that was also used to simulate the specific neutrino- nucleon interactions (kinematics) in each event. An important detail to mention here is that the neutrino interaction models in GENIE contained several parameters that could be tuned to match expectations found in experimental data. These included form factors that describe the internal structure of the nucleons and thus affected the final cross sections calculated for all processes. There are two particular form factors that were identified by previous oscillation analyses in Ice- Cube [2] as having a sizable effect in fitting neutrino oscillations parameters: those related to CC quasi-elastic and CC resonant scattering. These form factors are functions of the invariant quantity ๐‘„ 2 = โˆ’(๐‘˜ โ€ฒ โˆ’ ๐‘˜) 2 (where k and ๐‘˜ โ€ฒ are the 4-momenta of the incoming neutrino and outgoing charged lepton respectively) and are parametrized as: 1 ๐‘“ (๐‘„ 2 = โˆ’๐‘ž 2 ) โ‰ˆ (D.8) (1 โˆ’ (๐‘„ 2 /๐‘€ ๐ด2 )) 2 ๐ถ๐ถ ๐‘…๐ธ ๐‘† and ๐‘€ ๐ถ๐ถ ๐‘„๐ธ where ๐‘€ ๐ด is the so called axial mass [2]. The parameters ๐‘€ ๐ด ๐ด are included in the NMO analysis as nuisance parameters when fitting for the best fit mass ordering. In simulation, GENIE is used to recover the interaction weights that would be expected from the same event if the ๐ถ๐ถ ๐‘…๐ธ ๐‘† and ๐‘€ ๐ถ๐ถ ๐‘„๐ธ ๐‘€๐ด ๐ด parameters were shifted around their nominal value. An interpolation is ๐ถ๐ถ ๐‘…๐ธ ๐‘† and ๐‘€ ๐ถ๐ถ ๐‘„๐ธ performed to retrieve ๐‘€ ๐ด ๐ด dependent weights that can later be used in oscillation 134 analyses to predict how event rate expectations change as the ๐‘€ ๐ด parameters are shifted. Simulation weights could also be varied in the NMO analysis by shifting parameters related to the atmospheric neutrino flux model, as described in section 5.1.1. The total effective livetime of the generated atmospheric neutrino sample generated averaged over all energies, per neutrino flavor was of โˆผ 70 years, in an energy range between 1 GeV and 10 TeV. At the end of this simulation step, each simulated event contained detailed information on the energies and propagation directions of all neutrino interaction secondaries (cascades and outgoing leptons for CC events). D.2 Particle generation of atmospheric muon events Atmospheric muon simulations were generated using a dedicated software for this purpose, called MuonGun [103]. A description of the internal workings of this software is beyond the scope of this work. However, the final result was a simulation of atmospheric muon bundle events, weighted to an appropriate flux model. The effective livetime of this simulation was โˆผ 7.5 years [1] since these simulations are far more intensive, given the large atmospheric muon event rates prior to employing a dedicated neutrino selection criteria. D.3 Particle propagation The next step in the simulation chain is the propagation through the detector ice of the neutrino interaction secondaries and atmospheric muon bundles generated in the previous step. The development of the hadronic and electromagnetic cascades produced at the neutrino inter- action point were performed using GEANT4 [59], but only for hadronic cascades with energies below 30 GeV and electromagnetic cascades with energies below 100 MeV. For higher energy cascades, energy dependent average light yield parameterizations were used to predict the number of photons produced by each cascade and to draw the generation angles of each individual photon produced. These parametrizations were derived using dedicated GEANT4 simulations and are presented in [104]. Propagation of atmospheric muons was performed using a dedicated software package called PROPOSAL [105]. This software simulated the stochastic and continuous energy losses of the 135 muons, propagating them through the ice until they eventually stop and/or decay. Stochastic losses produce cascades that are treated as described previously, whereas parameterizations were used to predict the Cerenkov photons produced by the โ€™bare muonโ€™ and continuous energy loss secondaries, and to draw the generation angles of each individual photon. These parametrizations were also derived using dedicated GEANT4 simulations and are presented in [106]. This logic was also originally used in the ๐œˆ ๐œ‡ simulation. However, it was later modified for the purpose of neutrino/antineutrino classification in the NMO analysis presented here. Since PROPOSAL simulations do not include muon capture processes, muons were instead propagated directly using GEANT4. It was also decided to simulate cascades using GEANT4 at all energies in these simulations, since muon pairs produced in cascades could constitute a possible background to decay signals from the outgoing muon in CC interactions. Only about a third of all ๐œˆ ๐œ‡ events originally simulated were processed using this new scheme, given time constrains. As a result, the simulated ๐œˆ ๐œ‡ sample used in the NMO analysis had a lower effective livetime than the ๐œˆ๐‘’ and ๐œˆ๐œ samples. D.4 Photon Propagation and Detection The next step in the simulations was the propagation of Cherenkov photons through the ice, and their detection by the IceCube and DeepCore DOMs. It was run in the same fashion on atmospheric muon and atmospheric neutrino events. The modeling of the optical properties of the natural Antarctic ice, the refrozen hole ice, and the overall detection efficiency of the DOMs, were all factors that could modify the collection rates of photons in the detector. The baseline ice model used for the simulations was introduced in section 3.4. The baseline hole ice parametrization used values ๐‘ 0 = 0.101569 and ๐‘ 1 = โˆ’0.049344. These parameters are expected within a range โˆ’2.0 < ๐‘ 0 < 1.0 and โˆ’0.2 < ๐‘ 1 < 0.2 [1]. The baseline in-situ DOM efficiency arises from several effects including the PMT wavelength acceptance, the DOM angular acceptance and the local ice scattering [1]. A correction factor was applied in simulation to scale this efficiency. Uncertainties on the order of 10% with respect to the calculated DOM efficiencies are typically considered when performing oscillation analyses in IceCube. 136 Simulation set ID ๐‘0 ๐‘1 DOM efficiency factor Baseline set 0.101569 -0.049344 1.00 (1xx000) 1xx001 0.101569 -0.049344 0.90 1xx002 0.101569 -0.049344 0.95 1xx003 0.101569 -0.049344 1.05 1xx004 0.101569 -0.049344 1.10 1xx100 -0.0648 -0.1088 1.0 1xx101 -0.4839 -0.0171 1.0 1xx102 0.2803 -0.0754 1.0 1xx103 0.1122 0.0035 1.0 1xx104 -0.0498 -0.0543 0.88 1xx105 -0.3729 0.0349 0.93 1xx106 0.2965 -0.0363 0.97 1xx107 0.1244 -0.1132 1.03 1xx109 -0.3142 -0.0773 1.12 Atm. muon set 0.101569 -0.049344 1.00 (139011) Table D.1: The systematic MC sets generated for the purpose of neutrino oscillations analyses in IceCube with the neutrino sample used in the NMO analysis. The DOM efficiency factor is a correction factor multiplied to the nominal value of the overall in-situ DOM efficiency in the simulations. The values xx are to be replaced by 20,47,60 for ๐œˆ๐‘’ ,๐œˆ ๐œ‡ and ๐œˆ๐œ events. The ๐œˆ ๐œ‡ set has a different name convention since it was reprocessed to include muon capture physics. There were no systematic sets produced for atmospheric muon events for the purpose of the NMO analysis. Since the event rates in the detector change according to modeling of the detector ice and overall DOM efficiency, several systematic MC sets were generated to account for the effect of model variations on calculated neutrino oscillation sensitivities (see table D.1). No systematic sets were produced for atmospheric muon simulations for the purpose of the NMO analysis as their contribution is insignificant at the final analysis sample level (see Appendix E). Each MC set was then further processed to include simulations of PMT-noise, DOM readout, online processing at the SPS and the standard offline processing applied to data after it is sent to the Northern Hemisphere. A MC set of dedicated noise-only events (events triggered by noise signals) 137 was also generated to develop noise rejection strategies when developing the neutrino sample used in the NMO analysis. The set had an effective livetime of โˆผ 1 month [2]. 138 APPENDIX E DEEPCORE NMO ANALYSIS - ANALYSIS SAMPLE To achieve a sample of neutrino events in DeepCore, events due to noise and atmospheric muons, which dominate the sample after online processing at the SPS, have to be removed, until their rate is decreased by a factor of โˆผ 105 . The following provides a detailed description of the criteria used at each of the processing levels. Unless specified otherwise, the information provided in this chapter is based on the detailed review presented in [1]. The selection criteria was tested on the baseline MC simulation sets for atmospheric muons, atmospheric neutrinos, and noise triggered events, discussed in Appendix D prior to the reprocessing of the ๐œˆ ๐œ‡ sample. E.1 Level 2 (L2) Online processing at the SPS is commonly referred as level 1 (L1) processing within the collaboration and is common to all data. Similarly, standard offline filtering, commonly referred as L2 processing, is applied to all data after its transferred to the northern hemisphere. The filter at L2 relevant for the neutrino sample selected for this work is the DeepCore filter. It searches for events satisfying the SMT3 trigger [71] and produces a cleaned pulse series within each event, where pulses that were likely to have originated from noise are removed. The filter then loops over all pulses outside of the cleaned set, and adds them to the cleaned set if they are found within a designated radius R and time window T around a cleaned pulse (defined as R=150 m and T=1 ๐œ‡s). The procedure is repeated for the newly added cleaned pulses, until no new pulses are added. The resultant pulse series is called a "SRT-cleaned" series (Seeded RT) [107] and is further restricted to a time window around the SMT3 trigger. Only pulses within a static time window [-5000, 4000] ns around the SMT3 trigger are kept. The purpose of pulse cleaning is to remove isolated pulses, which are more likely to originate from dark noise. Pulses in the cleaned series are further divided into two sets: signal pulses located inside the DeepCore fiducial volume and veto pulses found in the baseline IceCube modules and the DeepCore veto cap. These two pulse sets are then used in an algorithm that identifies causal correlations between veto pulses and signal pulses, which are a clear indicator of an atmospheric muon reaching 139 the DeepCore volume. The algorithm first identifies the center-of-gravity (CoG) of the signal pulses, defined by the following equations [70]: ๐‘ 1 โˆ‘๏ธ ๐‘ฅยฎCoG = ๐‘ฅยฎ๐‘– , (E.1) ๐‘ ๐‘– ๐‘ 1 โˆ‘๏ธ  |ยฎ ๐‘ฅ๐‘– โˆ’ ๐‘ฅยฎCoG |  ๐‘กCoG = ๐‘ก๐‘– โˆ’ , (E.2) ๐‘ ๐‘ ice ๐‘– where ๐‘ฅยฎ๐‘– is the position of the DOM that detected the i-th pulse (in detector coordinates, with an origin roughly at the detector center), ๐‘ก๐‘– the time of the i-th pulse, N is the total number of pulses and ๐‘ ice the speed of light in ice. Veto pulses are then checked for causal correlations with the CoG, consistent with a particle traveling at a speed โˆผ ๐‘. Each veto pulse found to satisfy these correlations is marked and events with several of these will be ultimately discarded in the next processing level, since they are likely to contain atmospheric muon events. The neutrino sample used in this work generally only utilizes starting events, where secondaries are produced within the DeepCore volume as a result of a neutrino interaction [1]. E.2 Level 3 (L3) Several cuts are applied at this processing level to remove background events, including those that are usually improperly simulated and contribute to poor data/MC agreement: coincident muon events, muon bundles, and those satisfying the SLOP trigger [68]. Coincident muon events typically deposit light over a larger time window than neutrino events, since they contain two uncorrelated physical signals. Cuts are applied over two variables related to event duration to remove these events (see table E.1). For the selection criteria at L3, the extended DeepCore (see figure 3.3) is used to define the fiducial volume, while optical modules outside this region are used for vetoing. A large fraction of the atmospheric muon events (individual and bundles) depositing light inside DeepCore are expected to also deposit light in their path towards DeepCore, that is, in the veto region. Several cuts are therefore performed to remove events with a large fraction of veto pulses, particularly when pulses inside the DeepCore fiducial volume can be correlated to pulses in the veto region. For 140 Cut Variable Cut Value Background Description Cleaned time < 5000 ns Coincident events Difference between first and last pulse length time in cleaned series Uncleaned time < 13000 ns Coincident events Difference between first and last pulse length time NAbove200 < 10 Muons Number of pulses at Z > -200 m be- fore the DeepCore trigger time VertexGuess Z < -120 m Muons Z position of earliest pulse in cleaned series CausalVetoHits <7 Muons Number of veto pulses identified in the DeepCore filter at L2 Veto/Fiducial < 1.5 Muons Ratio of cleaned pulses outside/inside ratio the DeepCore fiducial volume C2HR6 > 0.37 Muons Fraction of all hit DOMs in the first 600 ns after the DeepCore trigger time RT Veto ==True Muons See description in main text NchCleaned โ‰ฅ6 Noise Number of DOMs with a cleaned pulse MicroCount >2 Noise The 300 ns time window in each event hits with most pulses is identified (in win- dow [-4,5] ๐œ‡s around DeepCore trig- ger). MicroCount hits corresponds to number of pulses in that window. Fiducial hits >2 Noise Number of DOMS that detected pulses in the cleaned series within the fiducial volume. NoiseEngine == True Noise See description in main text hits Table E.1: Cut variables at L3 processing, including the targeted background for each cut and a brief description of each variable [1]. Positions are displayed in detector coordinates. 141 example, there is a special cut that removes events based on the presence of clusters of light in the veto region by first exploiting the SRT-cleaning introduced in Appendix E.1 using each veto pulse as a seed to a SRT-cleaned series. The largest series generated this way is then marked as a cluster of light in the veto an used for a simple selection criteria, called the RT veto. The RT veto is passed when an event has more than 100 hit DOMs in the fiducial region and less than four pulses in its veto cluster, or when it has 4 pulses in its veto cluster but more than 75 hit DOMs in the fiducial volume [1]. Other cuts to remove atmospheric muons try to verify that the events started within DeepCore by checking the fraction of early pulses located inside its fiducial volume or the position of the earliest pulse. All cuts related to atmospheric muon removal are listed in table E.1. Noise triggered events are removed primarily by applying the assumption that events with a higher number of pulses in the cleaned series are more likely to contain a physics event. Most of the associated cuts involve simple variables, as summarized in table E.1, but one involves a more complex algorithm, the so called NoiseEngine. This particular algorithm looks for directionality in the cleaned pulse series light deposition pattern, which is expected for physics events. This is achieved by identifying the 300 ns time window within the event with a higher number of cleaned pulses and pairing each cleaned pulse within this window to all other cleaned pulses without pairs within the same DOM. A direction vector is calculated for each pair using time and direction information and a zenith and azimuth value retrieved. These are then binned and the distribution is evaluated for values above a certain threshold. If such bin can be found, the NoiseEngine criteria for directionality is satisfied [70]. Lastly, SLOP triggered events are removed to avoid events that are known to return poor data/MC agreement, and in the case of experimental data, an additional cut is included to remove events with so-called flaring DOMs, that is, instances of single optical modules with very high charge readings as a result of internal effects [1]. Rejection levels at L3 for atmospheric muons and noise were calculated using MC events and identified to be โˆผ 95% and above 99%, respectively. At this level, about 60% of the atmospheric 142 (a) (b) (c) Figure E.1: Distributions for three different cut variables used in the L3 selection. Distributions show noise, atmospheric muon background (MuonGun) and neutrino rates in MC compared to total rates in data [1]. Pass1/2 refers to data prior to and after a re-calibration of the DAQ in IceCube (in 2015) to achieve better data/MC agreement in the single photoelectron curve for DOMs [108]. Data ratios to noise, muon and total event rates in MC are also presented. a) Distribution for number of pulses at Z > -200 m before the DeepCore trigger time, introduced to remove muon backgrounds. Notice the improved data/MC agreement as well. b) Distribution of number of DOMs with a cleaned pulse, introduced to remove noise triggered events. This particular cut doesnโ€™t improve data/MC agreement. c) Distribution for time difference between first and last pulse in cleaned series, introduced to remove coincident muon events. Notice the high data/MC disagreement in the removed region, a result of poor simulation of coincident events. 143 Feature Description NchCleaned See table E.1 MicroCount hits Similar to the variable listed in table E.1, but we look for a 200 ns time window with most pulses (inside window [-3.5,4] ๐œ‡s around DeepCore trigger) in this case. iLineFit speed Reconstructed particle speed, as calculated by the recon- struction algorithm iLineFit (described in Appendix F) Fill ratio from mean Starting from a reconstructed interaction vertex (see Ap- pendix F), the mean of the distances between pulses and vertex is calculated, defining a radius R. This variable is the ratio of hit DOMs to not hit DOMs inside a sphere of radius R, centered at the reco vertex (โˆผ 1 for spherical symmetry). FullTimeLength Ratio Ratio of cleaned time length to uncleaned time length (see definitions in table E.1). Table E.2: Table indicating all training features for the noise classifier at L4 (adapted from [1]). neutrino events remain. Overall event rates before and after L3 are 13903 ๐œ‡Hz and 547 ๐œ‡Hz, respectively (in MC), meaning only about 3.9% of all events pass the selection criteria. Further processing levels are nevertheless required, since neutrinos constitute only 1% of the remaining sample at L3, while noise and atmospheric muon events constitute โˆผ 6.7% and โˆผ 92.3% of the sample, respectively. Example distributions showing data/MC agreement in the accepted and rejected regions for several cut variables used in the L3 selection can be found in figure E.1. E.3 Level 4 (L4) After L3, overall data/MC agreement is improved significantly, allowing the use of machine learning-based algorithms to increase the purity of neutrino events in the sample. Two independent BDTs were trained for this level: a noise discrimination classifier and an atmospheric muon discrimination classifier. The noise classifier used the training features listed in table E.2. Simulations of noise triggered events (โˆผ 1 month effective live-time) and the neutrino MC sample were used for the training (33.3% of events) and testing (66.7% of events) of the BDT. Events weights were pre-processed 144 (a) (b) Figure E.2: BDT scores for the noise rejection classifier at L4 (a) and event rates relative to the cut value on the BDT scores for the same classifier (b) [1]. before training to allow balanced weighted samples (total weighted events for noise and neutrinos are equal). Events with a BDT score below 0.7 are cut at this level. About 99.2% of the noise events at L3 are removed after this cut, while about 96% of all neutrino events survive. The BDT score distributions are shown in figure E.2. The muon classifier used the training features listed in table E.3. Neutrino events surviving the noise classifier cut are used for training. However, instead of using simulations for muon events, the classifier was trained for discrimination against data events, since atmospheric muon levels are still about 99% of the sample after the noise classifier. There is an added benefit of obtaining classification against poorly simulated physics events that might be missing in the MC. Data from several years of processing is used for training, but only 3 months of data per year. The same fractions for training and testing neutrino events as in the noise classifier were used. Events with a BDT score below 0.65 are cut; resulting in about 94% of the atmospheric muon events being removed (from the MC sample), while 87 % of all neutrino events survive. BDT score distributions can be found in figure E.3 145 Feature Description NchCleaned See table E.1 NAbove200 See table E.1 CausalVetoHits See table E.1 RT Veto See table E.1 Accumulated time Amount of time required for 75% of an eventโ€™s charge to be accumulated (only uses pulses from the cleaned series). VICH An algorithm identifies a "causal veto region" in a two dimensional space based on the relative position and time of each pulse to the pulse closest to the trigger time. Sev- eral pulses within this region are expected for atmospheric muon events. The training feature is the number of DOMs with pulses inside this phase space veto region. CoG Z Z coordinate of the CoG calculated with the cleaned pulse series (see equation E.1). RMS Z RMS of the Z coordinate of all pulses in the cleaned series. Z travel Maximum vertical distance between two DOMs that de- tected pulses in the cleaned series First HLC rho Radial distance (in cylindrical coordinates) between the center-most string in IceCube and the DOM that detected the first HLC pulse in an event. Table E.3: Training features for the muon classifier at L4. Table shown in detector coordinates (adapted from [1]) After L4 cuts, neutrinos constitute about 13.7 % of the remaining sample, while the percentages of noise and atmospheric muon events are โˆผ 0.9 % and โˆผ 85.4 % respectively. E.4 Level 5 (L5) At this level the main focus is to finally obtain a neutrino dominated sample by removing a large fraction of the remaining atmospheric muon events. These events have survived veto related cuts and, thus, are mostly due to muon events that arrive to DeepCore without depositing significant light in the veto. This is possible given the geometrical arrangement of the baseline IceCube detector. There are corridors along the baseline array through which muons can travel largely undetected, as shown in figure 3.3, until they arrive to the DeepCore fiducial volume. L5 cuts mainly aim to identify and remove events where the is a high chance of a muon traveling through a corridor. 146 (a) (b) Figure E.3: BDT scores for the atmospheric muon rejection classifier at level 4 (a) and event rates relative to cut value on BDT scores for the same classifier (b) [1]. Corridor cuts depend on an algorithm that identifies all possible azimuthal directions through which a muon could arrive undetected to each particular DeepCore string. The algorithm finds the CoG of all pulses in the cleaned series within DeepCore, and identifies the closest string. IceCube strings near to the corridors related to the identified DeepCore string are checked, and if more than two DOMs in a single corridor have isolated pulses, the event is discarded. Another method to remove events suspected to be corridor atmospheric muons is to identify events with a tentative starting position (starting events within DeepCore are an indicative of neutrino events) well within the fiducial DeepCore volume. Events with estimated vertex beyond a radial distance of 150m from the DeepCore center, or outside the range Z = [-490,-220] m (in detector coordinates), are removed. Finally, more stringent cuts on the BDT classifier scores for both noise (> 0.85) and atmospheric muons (> 0.9) are applied at this level. After L5 cuts, neutrinos constitute about 68.4 % of all events in the sample, whereas atmospheric muons and noise are now 29.4 % and 2.2 % of the sample, respectively. Only โˆผ 9.6 % of all events at L4 survive. 147 E.5 Level 6 (L6) Information regarding L6 and level 7 is based on the discussion in [2], unless stated otherwise within the text. At L6, the SANTA reconstruction (see Appendix F.2) is performed and cuts applied based on the results. SANTA is the final directional reconstruction applied to the neutrino sample. An additional step within L5 processing, but directly related to L6, is the generation of a pulse series to be used for the SANTA reconstruction. This is important because the hypothesis used by this reconstruction is based on clusters of direct light pulses, that is non-scattered photon signals, and requires a minimum of 5 DOMs with identified direct pulses. This stringent requirement results in a selection cut that removes about a third of all events surviving L5. To identify a pulse series comprised of direct photons, the following pulse cleaning steps are implemented in succession, as explained in [109]: โ€ข The cleaning algorithm first marks strings with 3 or more hit DOMs. Only those strings are considered when constructing the direct light pulse series. โ€ข A cleaning routine is run independently over each string. It starts by calculating the median pulse time over each string and selecting only pulses within [-250,2000] ns. Very late or early pulses are likely to be due to scattered light or noise. โ€ข The selected pulse with highest charge, or hot-spot, is then identified. The position of the DOM that detected it is used as an estimate for the point of closest distance between the string and the particle traveling through the detector (or cascade). The time of this pulse is also saved as a reference "time of closest approach". โ€ข Pulses in the DOMs above and below the hot-spot are then tested. For each pulse, its difference in time and vertical distance to the hot-spot is used to calculate an expected speed of propagation between them. If this speed is not consistent with the speed of light (within some tolerance), the pulse is discarded. 148 โ€ข All accepted pulses are considered reference pulses and are compared to all later surviving pulses by checking that the expected speed of propagation between reference and other pulses is consistent with c. โ€ข Finally, all surviving pulses are tested for their distance to the closest DOM (within their string) that also has a surviving pulse. If this "closest DOM" is more than 6 DOMs without clean pulses away, the tested pulse is removed. The result is a pulse series of direct light clusters along a couple of select strings (typically those that were close enough to the propagating particle). Events are then reconstructed using SANTA and the "direct pulse" cleaned series. After removing events without the SANTA reconstruction, neutrinos constitute about 75.3% of all events in the sample, whereas atmospheric muons are now about 24.7% of the sample. E.6 Level 7 (L7) The final reconstructions are run at this level. First, the reconstructed endpoint under a finite track hypothesis is reconstructed for each event using a dedicated algorithm called FiniteReco, described in Appendix F.3. SANTA provides an interaction vertex that is used by the FiniteReco algorithm as a first guess for its interaction vertex. At this level there is both directional and length information related to the neutrino interaction secondaries. This information is then used to seed a final reconstruction called LEERA, that assumes a mixed hypothesis consisting of a cascade at the particle interaction vertex together with an outgoing track. The algorithm is used to infer the cascade and track energy associated with the event. At this point, the energy and direction of the neutrino interaction secondaries is known and can be used to probe the oscillation probabilities of the incoming atmospheric neutrinos. However, a final piece of information is still missing from this picture, the neutrino flavor. Indeed, oscillation probabilities are different for ๐œˆ๐‘’ and ๐œˆ ๐œ‡ atmospheric neutrinos and, as a result, expectations are different for cascade and track-like events. Particle Identification (PID) is therefore important to separate these two populations, and is achieved by running a previously trained PID BDT. This 149 Feature Description Leera PID Leeraโ€™s hypothesis consists of both a track and cascade, with a common vertex. Track length (L) and cascade energy are fitted simultaneously. A simple PID estima- tor derived from this reconstruction is the likelihood ratio between the best fit cascade+track hypothesis and the cas- cade only hypothesis (L=0). The larger the ratio, the more likely an event is to contains an outgoing muon (track sig- nature). SANTA PID SANTA can be run under a track or cascade hypothesis. The ๐œ’2 ratio between the best fit (modified) ๐œ’2 under each hypothesis (see Appendix G) can be used as a PID estimator, following a similar logic to the Leera PID. Reco Length Reconstructed length from LEERA. Evidently, events with larger reconstructed lengths are more likely to contain a track topology. ๐‘ฃ๐‘’๐‘Ÿ๐‘ก๐‘’๐‘ฅ ๐œŒ36 Radial distance between reconstructed interaction vertex and the central string in DeepCore. Reco Vertex Z Depth of the reconstructed vertex. ๐‘’๐‘›๐‘‘๐‘๐‘œ๐‘–๐‘›๐‘ก ๐œŒ36 Radial distance between reconstructed muon endpoint and the central string in DeepCore. Reco Endpoint Z Depth of the reconstructed muon endpoint. Table E.4: Table describing all training features used in the PID classifier that is applied at level 7 processing (adapted from [2]). classifier is trained on half of the events of the MC simulation used for the analysis and was tested to ensure the performance on the rest of the sample was comparable to that of the training sample. The input variables used to train this BDT are provided in table E.4 and were selected after dedicated studies than searched for variables that could help in the separation of cascade and track topologies. PID BDT score distributions can be found in figure E.4. The final step at L7 is a series of cuts aimed to generate a data sample that is useful for oscillation analyses. The criteria is aimed to keep only track-like events. It starts by the definition of the analysis energy (6.31 to 158.49 GeV), zenith and PID score ranges: 150 Figure E.4: BDT scores for the PID classifier used at L7, where ๐œˆ ๐œ‡ CC events are considered tracks and all other event types cascades [2]. โ€ข ๐ธ โˆˆ [6.31, 158.49]GeV โ€ข cos(๐œƒ) โˆˆ [โˆ’1., 0.1] โ€ข PID > 0.55 (removed most cascade-like events, as seen in figure E.4) These ranges are chosen to allow sensitivity to the oscillation probabilities related to ๐œˆ ๐œ‡ disap- pearance [2] and events outside them are removed. Afterwards, a cut on the quality of the SANTA reconstruction was applied, and finally two cuts are applied to further reduce atmospheric muon contamination. First, a more stringent cut on the atmospheric muon rejection algorithm at L4 (BDT > 0.97), and then a cut aimed at removing neutrino events where an atmospheric muon reaches the detector simultaneously. This second cut removed events where there are pulses in the topmost IceCube DOMs consistent with a downgoing signal or with more than 7 pulses in the outermost strings in the detector. 151 Figure E.5: Showing Data/MC agreement in cosine zenith (left) and energy (right) after scaling distributions to match total data and MC event rates [2]. The discrepancy at the higher analysis range energies was considered acceptable since some systematics accounted for in oscillation analyses are energy dependent. After L7, the noise contamination is negligible and the muon contamination is only 2.1%. Sim- ulated event rates were found to be approximately 103 ๐œ‡Hz using the baseline MC set, atmospheric flux model, neutrino interaction model parameters and best fit oscillation parameters (under NO) from table 2.1 [2]. In comparison, event rates with experimental data were found to be about 93 ๐œ‡Hz [2]. Data/MC agreement is observed once event rates are scaled to match (see figure E.5) and thus the uncertainty in the total event rate is accounted for as a systematic (see section 5.1.1). 152 APPENDIX F RECONSTRUCTION ALGORITHMS The reconstruction algorithms used to retrieve the NMO analysis sample, and consequently in the NMO analysis presented in this work, are described in this section. Reconstructions are used to retrieve the physical properties of neutrino interaction secondaries produced in an event (primarily cascades or muon tracks), including energy, direction (zenith and azimuth), vertex, and time, among other relevant quantities. This is important for oscillation analyses, where the observables of interest (the energy and direction of the atmospheric neutrinos) are correlated to the energy and direction of the neutrino interaction secondaries. The reconstructions utilized to retrieve the analysis sample reviewed in Appendix E used a sequence of increasingly complex algorithms, where simpler ones were used as first guesses, or seeds, for more complex hypotheses. F.1 LineFit The following description of the LineFit algorithm is based on the review presented in [110]. LineFit is a simple reconstruction algorithm that ignores the cone geometry typical of Cherenkov emission by a moving particle in the ice. Instead, a simpler approximation is explored: a plane wave produced at the interaction vertex (ยฎ ๐‘ฅ0 and ๐‘ก0 ) and traveling with the particle, at its speed ๐‘ฃยฎ. Under the LineFit hypothesis, at some time, t, the light wavefront always includes the point: ๐‘ฅยฎ(๐‘ก) = ๐‘ฅยฎ0 + (๐‘ก โˆ’ ๐‘ก0 )ยฎ๐‘ฃ . (F.1) When an optical module at a given position ๐‘ฅยฎ๐‘› detects a pulse at some time ๐‘ก ๐‘› , LineFit approximates the location of the receiving optical module at ๐‘ฅยฎ(๐‘ก ๐‘› ). The algorithm then performs a least-squares fit, minimizing the function: ๐‘ โˆ‘๏ธ ๐œ’2 = ๐‘ฅ ๐‘› โˆ’ ๐‘ฅยฎ0 โˆ’ (๐‘ก๐‘– โˆ’ ๐‘ก0 )ยฎ๐‘ฃ ) 2 , (ยฎ (F.2) ๐‘› where all SRT-cleaned pulses (see appendix E.1) were included in this term when running the algorithm with the analysis sample, to determine the event vertex and direction. This minimization 153 problem has a closed analytic solution: ๐‘ฅยฎ0 = โŸจยฎ ๐‘ฅ ๐‘› โŸฉ โˆ’ โŸจ๐‘ก ๐‘› โŸฉยฎ๐‘ฃ , (F.3) โŸจ๐‘ก ๐‘› ๐‘ฅยฎ๐‘› โŸฉ โˆ’ โŸจ๐‘ก ๐‘› โŸฉโŸจยฎ ๐‘ฅ๐‘› โŸฉ ๐‘ฃยฎ = , (F.4) โŸจ๐‘ก ๐‘›2 โŸฉ โˆ’ โŸจ๐‘ก ๐‘› โŸฉ 2 (where the terms within brackets represent averages over all SRT-pulses), making LineFit an exceptionally fast algorithm that can be used as a seed for more complex reconstructions. Note that directionality can be inferred from the speed vector components. The expected results for tracks and cascades are โˆฅยฎ๐‘ฃ โˆฅ โ‰ˆ 0 and โˆฅยฎ๐‘ฃ โˆฅ = ๐‘ respectively . The LineFit algorithm used within the analysis sample processing is run at L2. F.2 SANTA The SANTA reconstruction consists of two distinct hypotheses for the light deposition patterns of tracks and cascades (reviewed in detail in Appendix G) along a single detector string. In simple terms, the track hypothesis comes from the intersection of a Cherenkov cone with a string, whereas the cascade hypothesis assumes a point source of light. These hypotheses are based entirely on geometric considerations and are valid only for unscattered Cherenkov photons [111]. The algorithm identifies strings with several direct pulses (a pulse series with only direct pulses is found in advance using the algorithm reviewed in Appendix E.5) and compares the observed light deposition pattern along each of those strings to the expected light deposition pattern from a given track or cascade hypothesis (defined by the vertex and direction of the particle traveling through the detector). A loss function quantifying the level of agreement between hypotheses and observations (equation G.21) is then minimized to identify the reconstructed vertex and direction (for tracks only) of the event. The SANTA algorithm can be run in two modes: single string (SS) fit or multiple string (MS) fit. For SS, the hypotheses are symmetric with respect to rotation around the string axis (z-axis in IceCube coordinates) [111]. The variables fit with the track hypothesis are interaction vertex (x,y,z and time) and direction (zenith and azimuth). Azimuth and x-y vertex positions cannot be resolved with a SS fit. In this case, the direction reconstructed with LineFit is used to identify the azimuth 154 SANTA error at analysis level SANTA error at analysis level 60 60 e+ e CC median + CC median 50 ยฑ25% 50 ยฑ25% 40 40 ] 30 ] 30 true [ 20 true [ 20 reco 10 reco 10 0 0 10 10 20 20 101 102 101 102 E [GeV] E [GeV] SANTA error at analysis level SANTA error at analysis level 60 60 + CC median All + NC median 50 ยฑ25% 50 ยฑ25% 40 40 ] 30 ] 30 true [ 20 true [ 20 reco 10 reco 10 0 0 10 10 20 20 101 102 101 102 E [GeV] E [GeV] Figure F.1: Showing SANTAโ€™s directional resolution for different event types in the analysis sample used for the NMO analysis discussed in this work. The ๐œˆ ๐œ‡ distribution corresponds to the reprocessed MC sample with muon capture processes included. direction and from there the x-y position of the vertex. For the cascade hypothesis, only the vertex is reconstructed [2]. A SS fit requires at least 5 direct pulses within a single string, whereas a MS fit requires two strings with at least 3 direct pulses per string. The criteria used in the analysis sample was to always prefer MS fits over SS fits. Nevertheless, SS fits are used in events where the conditions for a MS fit are not satisfied [2]. Figure F.1 illustrates SANTAโ€™s directional resolution in the analysis sample at L7. F.3 FiniteReco FiniteReco (FR) is an algorithm dedicated exclusively to reconstructing start and stop points along a track hypothesis. It works in two independent stages: a first guess algorithm, referred 155 Figure F.2: The procedure used by I3StartStopPoint to identify the points along an infinite track seed hypothesis corresponding to the start and stop point of the propagating particle [113]. to as I3StartStopPoint, that relies purely on geometrical considerations and a likelihood based reconstruction, that uses the first guess as seed. The I3StartStopPoint algorithm requires a seed infinite track reconstruction and a cleaned pulse series. For the analysis sample, these are the track reconstruction with SANTA and the SRT-cleaned pulse series. All DOMs in a radius R (=300 m) around the seed track are checked for pulses within the cleaned pulse series. Identified hit DOMs are then projected onto the track using the Cherenkov angle and the two projected points further away from each other are identified as the start and stop point (direction is used to label them properly) [112]. The procedure is illustrated in figure F.2. The FR likelihood-based reconstruction compares the light depositions in an event to expected light depositions at each optical module given a finite track hypothesis. Each hypothesis is uniquely defined by a set of parameters { ๐‘Ž},ยฎ that includes the vertex time and position, the track direction and its length [112]. In an iteration of the algorithm, only length is fit while all other parameters are kept constant at their seed values (seeded with I3StartStopPointโ€™s vertex and length and SANTAโ€™s direction within the oscNext processing). The likelihood function used in the FR algorithm can be understood by first considering a single DOM at a position ๐‘ฅยฎ๐‘– . Lets assume that given a set of parameters { ๐‘Ž}, ยฎ we can estimate the mean number of photons expected to arrive at ๐‘ฅยฎ๐‘– , called ๐œ†๐‘– ( ๐‘Ž). ยฎ The detected number of counts, ๐‘›๐‘– , in the 156 DOM will then follow a Poisson distribution of the form: ๐‘› ๐œ†๐‘– ๐‘– ๐‘’ โˆ’๐œ†๐‘– ๐‘ƒ(๐‘›๐‘– |๐œ†๐‘– ( ๐‘Ž)) ยฎ = . (F.5) ๐‘›๐‘– ! The next step is to then consider the probability that the DOM doesnโ€™t detect a pulse, which is simply given by: ๐‘ƒ๐‘– (no hit| ๐‘Ž) ยฎ = ๐‘ƒ(๐‘›๐‘– = 0|๐œ†๐‘– ( ๐‘Ž))ยฎ = ๐‘’ โˆ’๐œ†๐‘– ( ๐‘Ž)ยฎ , (F.6) and the probability the DOM detects one or more pulses, given by 1 โˆ’ ๐‘ƒ๐‘– (no hit| ๐‘Ž). ยฎ The (approx- imately constant) noise rate of each individual DOM can be used to infer the expected number of noise pulses ๐œŒ๐‘– in the DOM, and from there the corrected hit/no hit probabilities, where ๐œ†๐‘– ( ๐‘Ž) ยฎ is replaced by ๐œ†๐‘– ( ๐‘Ž) ยฎ + ๐œŒ๐‘– . As explained in [112], FR would use a hit/no hit likelihood of the form: ร– ร– ยฎ = ๐ฟ( ๐‘Ž) [1 โˆ’ ๐‘๐‘– (no hit| ๐‘Ž, ยฎ ๐œŒ๐‘– )] [ ๐‘๐‘– (no hit| ๐‘Ž, ยฎ ๐œŒ๐‘– )], (F.7) ๐‘–โˆˆ๐ป ๐‘–โˆˆ๐‘€ where H and M are the sets of all hit and non hit DOMs respectively. The likelihood function is minimized over each event to derive the best fit value for the particle length. This procedure is, in practice, performed twice once keeping the vertex at a constant point (endpoint reconstruction), and a subsequent time keeping the stop point constant (vertex reconstruction). The remaining detail about FR is how the light expectations ๐œ†๐‘– ( ๐‘Ž) ยฎ can be calculated. There are several options, as reviewed in [112]. The one chosen for the analysis sample processing and relevant to the work presented in this thesis is the usage of so-called "photonrec" tables [114]. These tables are generated by simulation of point-like Cherenkov emitters with unit energy, each described by a propagation direction (zenith) and a position (depth). Dedicated simulations were run where photons are propagated directly though a detector medium for a set of different position and orientation combinations, and for each the expected photon flux at several positions relative to the emitter are calculated and the results tabulated. The tabulated results can then be interpolated to derive expected photon fluxes at an arbitrary location within the detector, given a point-like Cherenkov emitter with some arbitrary position and direction. To obtain the expected light fluxes 157 LEERA error at analysis level LEERA error at analysis level 1.00 1.00 e+ e CC median + CC median 0.75 ยฑ25% 0.75 ยฑ25% 0.50 0.50 0.25 0.25 (E , reco E )/E (E , reco E )/E 0.00 0.00 0.25 0.25 0.50 0.50 0.75 0.75 1.00 1.00 101 102 101 102 E [GeV] E [GeV] LEERA error at analysis level LEERA error at analysis level 1.00 1.00 + CC median All + NC median 0.75 ยฑ25% 0.75 ยฑ25% 0.50 0.50 0.25 0.25 (E , reco E )/E (E , reco E )/E 0.00 0.00 0.25 0.25 0.50 0.50 0.75 0.75 1.00 1.00 101 102 101 102 E [GeV] E [GeV] Figure F.3: Showing LEERAโ€™s energy resolution for different event types in the analysis sample used for the NMO analysis discussed in this work. The ๐œˆ ๐œ‡ distribution corresponds to the reprocessed MC sample with muon capture processes included. from a finite muon track, the expected light yields of the infinitesimal point like emitters are integrated over the muon length [114]. The advantage of the photon tables is that they include ice properties (scattering and absorption) in the calculation of expected light yields, although the cost is they are computationally intensive to generate. Energy can additionally be included in the light yield calculations to scale the intensity of the photon flux. For minimum ionizing muon hypotheses, however, it can calculated directly from the muon length. The finite muon hypothesis within FR, as can be seen from the description above, assumes minimum ionizing muon tracks. No description of stochastic energy losses is included. 158 F.4 LEERA Reconstructions discussed so far can be used to estimate the vertex, direction, time and length of a track (or cascade) within the detector. However, none were designed to directly infer the event energy. This final step in the reconstruction pipeline within oscNext is performed by the LEERA likelihood-based algorithm. LEERA uses a mixed hypothesis, consisting of a cascade at the particle interaction vertex together with an outgoing track (starting at the same vertex). The reconstruction is seeded with both the start point and length from FR and the direction found using SANTA. Both cascade energy and track length (stop point fixed, interaction vertex allowed to move) are fit by minimizing a hit/no hit likelihood (as in equation F.7). Notice this algorithm, once again, uses a minimum ionizing track hypothesis; thus, track energy can be inferred from the track length. A fundamental difference in the likelihood term minimized with LEERA, with respect to FR, is the calculation method for the light yields. Both algorithms use table-based light yields, but where FR uses photonrec tables, LEERA uses the more recently generated spline tables [115]. These are produced by approximation of the multidimensional photon tables using basis splines and, as explained in [115], help reduce numerical artifacts that arise from binning during tabulation and the memory requirements for table storage. Figure F.3 illustrates LEERAโ€™s energy resolution in the analysis sample at L7. 159 APPENDIX G SANTAโ€™S TRACK AND CASCADE FITS The following discussion is fully based on the description for the SANTA track and cascade hypotheses that is presented in [111] and [2]. G.1 Track Hypothesis To understand the track reconstruction performed by SANTA, consider an infinite track hypoth- esis, defined by a single straight line, representing a charged particleโ€™s trajectory. The position of the particle at some time t is given by: ยฎ = ๐‘ยฎ0 + ๐‘(๐‘ก โˆ’ ๐‘ก0 ) ๐‘ข, ๐‘(๐‘ก) ห† (G.1) where ๐‘ยฎ0 = ๐‘(๐‘กยฎ 0 ) is a point along the trajectory, ๐‘ก 0 is some reference time and ๐‘ขห† a unit vector in the direction of motion. Notice the hypothesis assumes a particle traveling at the speed of light c. On the other hand, a point along a string, at height z (in detector coordinates), can be described by a position vector: (๐ฟ ๐‘ฅ , ๐ฟ ๐‘ฆ , ๐‘ง), (G.2) where ๐ฟ ๐‘ฅ , ๐ฟ ๐‘ฆ are constant (x and y positions of the string in detector coordinates). The line segment of closest approach between two lines is perpendicular to both of them. A unit vector connecting the points of closest approach (pointing from the string towards the track) will be defined by: ๐‘›ห† = ๐‘ขห† ร— ๐‘งห†, (G.3) as long as the track is not parallel to the string (in which case any position along the track corresponds to that of closest approach). Given the (horizontal) line segment of closest approach connects the two line segments, we find the following equation relating the point of closest approach along the track ๐‘(๐‘ก ยฎ ๐‘ ) to the position in the string at the height of closest approach ๐‘ง ๐‘ : (๐ฟ ๐‘ฅ , ๐ฟ ๐‘ฆ , ๐‘ง ๐‘ ) + ๐‘›๐‘‘ ยฎ ๐‘ ), ห† ๐‘ = ๐‘(๐‘ก (G.4) 160 where ๐‘‘ ๐‘ , ๐‘ก ๐‘ are the distance and time of closest approach respectively. This system of equations can be solved to derive the following equations, relating the variables we want to reconstruct (direction and vertex) to the distance, time and depth of closest approach between the track and the string: 1h i ๐‘ก ๐‘ = ๐‘ก0 + ๐ฟ ๐‘ฅ sin(๐œƒ) cos ๐œ™ + ๐ฟ ๐‘ฆ sin(๐œƒ) sin ๐œ™ + ๐‘ง ๐‘ cos(๐œƒ) โˆ’ ๐‘ยฎ0 ยท ๐‘ขยฎ , (G.5) ๐‘ ๐‘ 0,๐‘ง โˆ’ cos(๐œƒ)( ๐‘ยฎ0 ยท ๐‘ข) ยฎ + (๐ฟ ๐‘ฅ sin(๐œƒ) cos ๐œ™ + ๐ฟ ๐‘ฆ sin(๐œƒ) sin ๐œ™) cos(๐œƒ) ๐‘ง๐‘ = , (G.6) sin2 (๐œƒ) โˆš๏ธƒ ๐‘‘๐‘ = ( ๐‘ ๐‘ฅ (๐‘ก ๐‘ ) โˆ’ ๐ฟ ๐‘ฅ ) 2 + ( ๐‘ ๐‘ฆ (๐‘ก ๐‘ ) โˆ’ ๐ฟ ๐‘ฆ ) 2 , (G.7) where we have expressed the direction vector of the track in spherical coordinates, with ๐œƒ and ๐œ™ the zenith and azimuth angles respectively (defined as the direction the particle is traveling towards to). The equations hold for non-vertical tracks. For the case in which a track travels vertically, there is freedom to choose any position along the track and time as those of closest approach. The choice in SANTA is: ๐‘ง ๐‘ = ๐‘ 0,๐‘ง and ๐‘ก ๐‘ = ๐‘ก0 . Dependence on azimuth angle disappears from these equations when fitted over a single string. In that case we are free to define the coordinate system centered at the string position (๐ฟ ๐‘ฅ ,๐ฟ ๐‘ฆ ) = (0,0), which results in the following simplified equations: 1h i ๐‘ก ๐‘ = ๐‘ก0 + ๐‘ง ๐‘ cos(๐œƒ) โˆ’ ๐‘ยฎ0 ยท ๐‘ขยฎ , (G.8) ๐‘ ๐‘ 0,๐‘ง โˆ’ cos(๐œƒ) ( ๐‘ยฎ0 ยท ๐‘ข) ยฎ ๐‘ง๐‘ = , (G.9) sin2 (๐œƒ) โˆš๏ธƒ ๐‘‘๐‘ = ( ๐‘ ๐‘ฅ (๐‘ก ๐‘ )) 2 + ( ๐‘ ๐‘ฆ (๐‘ก ๐‘ )) 2 , (G.10) where: ๐‘ 0,๐‘ง โˆ’ sin2 (๐œƒ)๐‘ง ๐‘ ๐‘ยฎ0 ยท ๐‘ขยฎ = , (G.11) cos ๐œƒ is independent on the azimuth angle as well. The cylindrical symmetry doesnโ€™t allow us to uniquely identify the x,y coordinates of the particle vertex. The next step involves deriving the expected arrival times of photons into the string, at a height z, given some particle direction and vertex. The main geometrical feature that SANTA exploits 161 Figure G.1: Schematic showing the geometry considered by SANTA when fitting over a single string [111]. is the fact that these times can be fully parameterized in terms of the distance, time and depth of closest approach, together with the zenith angle of the particle. This allows to fit for these variables using photon arrival times. Notice the azimuth angle cannot be determined with a single string as a result. To obtain the detection time of photons, start by considering figure G.1. For simplicity,when fitting over a single string, lets set (๐ฟ ๐‘ฅ ,๐ฟ ๐‘ฆ ) = (0,0). Imagine a photon is produced at a cherenkov angle ๐œƒ ๐‘ at some time ๐‘ก โˆ— . Notice the vector along the path of the photon starts at a position ๐‘(๐‘ก ยฎ โˆ—) and ends at the detection point along the string, at a height z. The dot product of this vector and 162 the direction vector ๐‘ขยฎ is given by: ๐‘‘ ๐›พ (๐‘ง) cos(๐œƒ ๐‘ ) = [(0, 0, ๐‘ง) โˆ’ ๐‘(๐‘ก ยฎ โˆ— )] ยท ๐‘ขยฎ = โˆ’ ๐‘(๐‘กยฎ โˆ— ) ยท ๐‘ขยฎ + ๐‘ง cos(๐œƒ), (G.12) where ๐‘‘ ๐›พ (๐‘ง) is the distance travelled by the photon. It can be show using equation G.1 at ๐‘ก โˆ— and equation G.8, that: ยฎ โˆ— ) ยท ๐‘ขยฎ = ๐‘ง ๐‘ cos(๐œƒ) + (๐‘ก โˆ— โˆ’ ๐‘ก ๐‘ )๐‘ ๐‘(๐‘ก (G.13) Using equations 3.1, G.12 and G.13 we can then express the generation time of the photon in terms of the travelled photon distance, the zenith angle and the closest approach variables: 1 ๐‘ก โˆ— = ๐‘ก ๐‘ + [(๐‘ง โˆ’ ๐‘ง ๐‘ ) cos(๐œƒ) โˆ’ ๐‘‘ ๐›พ (๐‘ง)/๐‘›], (G.14) ๐‘ where n is the phase index of refraction of the medium. Since the difference between the detection time of the photon ๐‘ก ๐›พ and the generation time of the photons ๐‘ก โˆ— is the propagation time, and photons travel at a speed c/n is ice, we finally have: 1  1 ๐‘ก ๐›พ (๐‘ง) = ๐‘ก ๐‘ + [(๐‘ง โˆ’ ๐‘ง ๐‘ ) cos(๐œƒ) + ๐‘› โˆ’ ๐‘‘ ๐›พ (๐‘ง)], (G.15) ๐‘ ๐‘› From similar geometric considerations it can also be shown that (see [111]): โˆš๏ธƒ ๐‘› ๐‘‘ ๐›พ (๐‘ง) = โˆš ๐‘‘ ๐‘2 + (๐‘ง โˆ’ ๐‘ง ๐‘ ) 2 sin2 (๐œƒ). (G.16) ๐‘›2 โˆ’ 1 Equations G.15 and G.16 return the expected detection times (๐‘ก ๐›พ ) of direct photons along the string, given a track hypothesis (detection time curves obtained for an example track hypothesis can be found in figure G.2). Notice the track is parameterized in these equations as a function of ๐‘ก ๐‘ , ๐‘ง ๐‘ , ๐‘‘ ๐‘ and the zenith angle, as mentioned before. Single string (SS) fits in SANTA use photon detection times across each individual string to infer these four parameters. Fits that rely on multiple strings (MS), on the other hand, will be able to disentangle the vertex and direction fully from these four parameters, from the relations G.5 to G.7 over several strings. 163 Figure G.2: The figure shows the times of detection of photons for a track hypothesis with ๐‘‘ ๐‘ = 10. m, ๐‘ง ๐‘ = 0 and ๐‘ก ๐‘ = 0, for several zenith angles, as defined by equations G.15 and G.16. In the figure, ๐œƒ = 0โ—ฆ for muons traveling upwards, that is, towards the detector surface. Horizontal and downgoing muons defined by 90โ—ฆ and 180โ—ฆ respectively. In addition, SANTA fits over the expected charge in each DOM, given the expected arrival angles of photons into the DOMs and their angular acceptance profile. The arrival angle ๐œƒ ๐›พ (as shown in figure G.1), can once again be derived from geometrical considerations and is found in [111] to be: ๐‘ง โˆ’ ๐‘ง ๐‘ ๐‘๐‘œ๐‘ (๐œƒ) cos ๐œƒ ๐›พ = sin2 (๐œƒ)  + , (G.17) ๐‘‘ ๐›พ (๐‘ง) ๐‘› where again, this quantity can be parameterized in terms of the closest approach variables and the zenith angle. G.2 Cascade Hypothesis The cascade hypothesis for SANTA assumes a point source of light, at the interaction point. It parametrizes the interaction time and vertex as those of closest approach. In this case, it can be easily shown that the following relations hold for the detection times and expected arrival angles of photons into the string: ๐‘› ๐‘ก ๐›พ (๐‘ง) = ๐‘ก ๐‘ + ๐‘‘ ๐›พ (๐‘ง), (G.18) ๐‘ 164 Figure G.3: The figure shows the times of detection of photons for cascade and horizontal track hypothesis with ๐‘‘ ๐‘ = 10. m, ๐‘ง ๐‘ = 0 and ๐‘ก ๐‘ = 0. โˆš๏ธƒ ๐‘‘ ๐›พ (๐‘ง) = ๐‘‘ ๐‘2 + (๐‘ง โˆ’ ๐‘ง ๐‘ ) 2 , (G.19)  (๐‘ง โˆ’ ๐‘ง ๐‘ ) cos ๐œƒ ๐›พ = (G.20) ๐‘‘ ๐›พ (๐‘ง) Notice no directionality is fitted when using this hypothesis. For multiple strings, the x and y positions of the interaction point can be disentangled, since equation G.7 will hold for each individual string. The profile of arrival times for photons over a string under a cascade hypothesis is exemplified in figure G.3, where there is also a comparison to the light profile for a horizontal track with the same closest approach parameters. Usually both track and cascade hypotheses are tested when using SANTA. G.3 Fitting To fit for the the closest approach variables (and the particle direction for track hypotheses), the SANTA algorithm minimizes over the following loss function: ๐‘‘โ€ฒ ๐‘žโ€ฒ ๐‘ โˆ‘๏ธ h ๐‘ก ๐›พ (๐‘ง ๐‘› |ฮ˜) โˆ’ ๐‘ก ๐‘›  2 i ๐‘„(ฮ˜) = ฮฆ ๐ถ2 + ๐‘› ร ๐‘› , (G.21) ๐ถ๐œŽ๐‘ก ๐‘‘0 1 ๐‘ ๐‘žโ€ฒ ๐‘›=1 ๐‘ ๐‘˜=1 ๐‘˜ โˆš๏ธƒ ๐‘‘๐‘›โ€ฒ = ๐‘‘12 + ๐‘‘ ๐›พ2 , (G.22) 165 2๐‘ž ๐‘› ๐‘žโ€ฒ๐‘› =  (G.23) cos ๐œƒ ๐›พ + 1 where ฮ˜ represents the set of parameters that define the track or cascade hypothesis. The left term in this function is a comparison between the expected (๐‘ก ๐›พ ) and detected (๐‘ก ๐‘› ) photon times, for each optical module with a cleaned pulse (cleaning for direct light pulses described in Appendix E.5),   where ฮฆ(๐‘… 2 ) = log 1 + ๐‘… 2 is a loss function applied to the detection time residuals, ๐œŽ๐‘ก is the time uncertainty of pulses and C is a scaling factor. The right term in equation G.21 is a penalty term to the product of the detected charge on the DOM (๐‘žโ€ฒ๐‘› ; corrected by the DOM angular acceptance profile) to the expected travel distance of the photons (๐‘‘๐‘›โ€ฒ ; corrected to include a minimum distance ๐‘‘1 , roughly equivalent to the DOM radius). If the product is too large, the DOM observes charge that would be expected from a photon travel path much shorter than the one being tested. In those cases, the penalty term is also large and the hypothesis disfavoured. The product need to be large relative to a baseline expected product between these two quantities. This is given by the average detected charge from all pulses considered by the algorithm, multiplied by a normalization factor ๐‘‘0 = 7m. The fit can be performed over multiple strings or a single string. In the second case, the azimuth angle (and ultimately the x,y position of the vertex) are derived from a LineFit reconstruction. The single string fits require at least 5 pulses on a string, while the multiple string fits need two strings minimum, each with at least 3 pulses. 166