Articles | Volume 18, issue 1
https://doi.org/10.5194/tc-18-75-2024
https://doi.org/10.5194/tc-18-75-2024
Research article
 | 
04 Jan 2024
Research article |  | 04 Jan 2024

In situ estimation of ice crystal properties at the South Pole using LED calibration data from the IceCube Neutrino Observatory

Rasha Abbasi, Markus Ackermann, Jenni Adams, Nakul Aggarwal, Juanan Aguilar, Markus Ahlers, Maryon Ahrens, Jean-Marco Alameddine, Antonio Augusto Alves Junior, Najia Moureen Binte Amin, Karen Andeen, Tyler Anderson, Gisela Anton, Carlos Argüelles, Yosuke Ashida, Sofia Athanasiadou, Spencer Axani, Xinhua Bai, Aswathi Balagopal V, Moreno Baricevic, Steve Barwick, Vedant Basu, Ryan Bay, James Beatty, Karl Heinz Becker, Julia Becker Tjus, Jakob Beise, Chiara Bellenghi, Samuel Benda, Segev BenZvi, David Berley, Elisa Bernardini, Dave Besson, Gary Binder, Daniel Bindig, Erik Blaufuss, Summer Blot, Federico Bontempo, Julia Book, Jürgen Borowka, Caterina Boscolo Meneguolo, Sebastian Böser, Olga Botner, Jakob Böttcher, Etienne Bourbeau, Jim Braun, Bennett Brinson, Jannes Brostean-Kaiser, Ryan Burley, Raffaela Busse, Michael Campana, Erin Carnie-Bronca, Chujie Chen, Zheyang Chen, Dmitry Chirkin, Koun Choi, Brian Clark, Lew Classen, Alan Coleman, Gabriel Collin, Amy Connolly, Janet Conrad, Paul Coppin, Pablo Correa, Stefan Countryman, Doug Cowen, Robert Cross, Christian Dappen, Pranav Dave, Catherine De Clercq, James DeLaunay, Diyaselis Delgado López, Hans Dembinski, Kunal Deoskar, Abhishek Desai, Paolo Desiati, Krijn de Vries, Gwenhael de Wasseige, Tyce DeYoung, Alejandro Diaz, Juan Carlos Díaz-Vélez, Markus Dittmer, Hrvoje Dujmovic, Michael DuVernois, Thomas Ehrhardt, Philipp Eller, Ralph Engel, Hannah Erpenbeck, John Evans, Paul Evenson, Kwok Lung Fan, Ali Fazely, Anatoli Fedynitch, Nora Feigl, Sebastian Fiedlschuster, Aaron Fienberg, Chad Finley, Leander Fischer, Derek Fox, Anna Franckowiak, Elizabeth Friedman, Alexander Fritz, Philipp Fürst, Tom Gaisser, Jay Gallagher, Erik Ganster, Alfonso Garcia, Simone Garrappa, Lisa Gerhardt, Ava Ghadimi, Christian Glaser, Thorsten Glüsenkamp, Theo Glauch, Noah Goehlke, Javier Gonzalez, Sreetama Goswami, Darren Grant, Shannon Gray, Timothée Grégoire, Spencer Griswold, Christoph Günther, Pascal Gutjahr, Christian Haack, Allan Hallgren, Robert Halliday, Lasse Halve, Francis Halzen, Hassane Hamdaoui, Martin Ha Minh, Kael Hanson, John Hardin, Alexander Harnisch, Patrick Hatch, Andreas Haungs, Klaus Helbing, Jonas Hellrung, Felix Henningsen, Lars Heuermann, Stephanie Hickford, Colton Hill, Gary Hill, Kara Hoffman, Kotoyo Hoshina, Wenjie Hou, Thomas Huber, Klas Hultqvist, Mirco Hünnefeld, Raamis Hussain, Karolin Hymon, Seongjin In, Nadege Iovine, Aya Ishihara, Matti Jansson, George Japaridze, Minjin Jeong, Miaochen Jin, Ben Jones, Donghwa Kang, Woosik Kang, Xinyue Kang, Alexander Kappes, David Kappesser, Leonora Kardum, Timo Karg, Martina Karl, Albrecht Karle, Uli Katz, Matt Kauer, John Kelley, Ali Kheirandish, Ken'ichi Kin, Joanna Kiryluk, Spencer Klein, Alina Kochocki, Ramesh Koirala, Hermann Kolanoski, Tomas Kontrimas, Lutz Köpke, Claudio Kopper, Jason Koskinen, Paras Koundal, Michael Kovacevich, Marek Kowalski, Tetiana Kozynets, Emmett Krupczak, Emma Kun, Naoko Kurahashi, Neha Lad, Cristina Lagunas Gualda, Michael Larson, Frederik Lauber, Jeffrey Lazar, Jiwoong Lee, Kayla Leonard, Agnieszka Leszczyńska, Massimiliano Lincetto, Qinrui Liu, Maria Liubarska, Elisa Lohfink, Christina Love, Cristian Jesus Lozano Mariscal, Lu Lu, Francesco Lucarelli, Andrew Ludwig, William Luszczak, Yang Lyu, Wing Yan Ma, Jim Madsen, Kendall Mahn, Yuya Makino, Sarah Mancina, Wenceslas Marie Sainte, Ioana Mariş, Szabolcs Marka, Zsuzsa Marka, Matthew Marsee, Ivan Martinez-Soler, Reina Maruyama, Thomas McElroy, Frank McNally, James Vincent Mead, Kevin Meagher, Sarah Mechbal, Andres Medina, Maximilian Meier, Stephan Meighen-Berger, Yarno Merckx, Jessie Micallef, Daniela Mockler, Teresa Montaruli, Roger Moore, Bob Morse, Marjon Moulai, Tista Mukherjee, Richard Naab, Ryo Nagai, Uwe Naumann, Amid Nayerhoda, Jannis Necker, Miriam Neumann, Hans Niederhausen, Mehr Nisa, Sarah Nowicki, Anna Obertacke Pollmann, Marie Oehler, Bob Oeyen, Alex Olivas, Rasmus Orsoe, Jesse Osborn, Erin O'Sullivan, Hershal Pandya, Daria Pankova, Nahee Park, Grant Parker, Ek Narayan Paudel, Larissa Paul, Carlos Pérez de los Heros, Lilly Peters, Josh Peterson, Saskia Philippen, Sarah Pieper, Alex Pizzuto, Matthias Plum, Yuiry Popovych, Alessio Porcelli, Maria Prado Rodriguez, Brandon Pries, Rachel Procter-Murphy, Gerald Przybylski, Christoph Raab, John Rack-Helleis, Mohamed Rameez, Katherine Rawlins, Zoe Rechav, Abdul Rehman, Patrick Reichherzer, Giovanni Renzi, Elisa Resconi, Simeon Reusch, Wolfgang Rhode, Mike Richman, Benedikt Riedel, Ella Roberts, Sally Robertson, Steven Rodan, Gerrit Roellinghoff, Martin Rongen, Carsten Rott, Tim Ruhe, Li Ruohan, Dirk Ryckbosch, Devyn Rysewyk Cantu, Ibrahim Safa, Julian Saffer, Daniel Salazar-Gallegos, Pranav Sampathkumar, Sebastian Sanchez Herrera, Alexander Sandrock, Marcos Santander, Sourav Sarkar, Subir Sarkar, Merlin Schaufel, Harald Schieler, Sebastian Schindler, Berit Schlüter, Torsten Schmidt, Judith Schneider, Frank Schröder, Lisa Schumacher, Georg Schwefer, Steve Sclafani, Dave Seckel, Surujhdeo Seunarine, Ankur Sharma, Shefali Shefali, Nobuhiro Shimizu, Manuel Silva, Barbara Skrzypek, Ben Smithers, Robert Snihur, Jan Soedingrekso, Andreas Søgaard, Dennis Soldin, Christian Spannfellner, Glenn Spiczak, Christian Spiering, Michael Stamatikos, Todor Stanev, Robert Stein, Thorsten Stezelberger, Timo Stürwald, Thomas Stuttard, Greg Sullivan, Ignacio Taboada, Samvel Ter-Antonyan, Will Thompson, Jessie Thwaites, Serap Tilav, Kirsten Tollefson, Christoph Tönnis, Simona Toscano, Delia Tosi, Alexander Trettin, Chun Fai Tung, Roxanne Turcotte, Jean Pierre Twagirayezu, Bunheng Ty, Martin Unland Elorrieta, Karriem Upshaw, Nora Valtonen-Mattila, Justin Vandenbroucke, Nick van Eijndhoven, David Vannerom, Jakob van Santen, Javi Vara, Joshua Veitch-Michaelis, Stef Verpoest, Doga Veske, Christian Walck, Winnie Wang, Timothy Blake Watson, Chris Weaver, Philip Weigel, Andreas Weindl, Jan Weldert, Chris Wendt, Johannes Werthebach, Mark Weyrauch, Nathan Whitehorn, Christopher Wiebusch, Nathan Willey, Dawn Williams, Martin Wolf, Gerrit Wrede, Johan Wulff, Xianwu Xu, Juan Pablo Yanez, Emre Yildizci, Shigeru Yoshida, Shiqi Yu, Tianlu Yuan, Zelong Zhang, and Pavel Zhelnin
Abstract

The IceCube Neutrino Observatory instruments about 1 km3 of deep, glacial ice at the geographic South Pole. It uses 5160 photomultipliers to detect Cherenkov light emitted by charged relativistic particles. An unexpected light propagation effect observed by the experiment is an anisotropic attenuation, which is aligned with the local flow direction of the ice. We examine birefringent light propagation through the polycrystalline ice microstructure as a possible explanation for this effect. The predictions of a first-principles model developed for this purpose, in particular curved light trajectories resulting from asymmetric diffusion, provide a qualitatively good match to the main features of the data. This in turn allows us to deduce ice crystal properties. Since the wavelength of the detected light is short compared to the crystal size, these crystal properties include not only the crystal orientation fabric, but also the average crystal size and shape, as a function of depth. By adding small empirical corrections to this first-principles model, a quantitatively accurate description of the optical properties of the IceCube glacial ice is obtained. In this paper, we present the experimental signature of ice optical anisotropy observed in IceCube light-emitting diode (LED) calibration data, the theory and parameterization of the birefringence effect, the fitting procedures of these parameterizations to experimental data, and the inferred crystal properties.

Dates
1 Introduction

The 2021 IPCC report (Intergovernmental Panel on Climate Change2021) highlights the need to understand the dynamics of ice sheets in order to predict their contribution to sea level rise in a changing climate. Ice flows under its own weight, either through basal sliding or through plastic deformation, which is mediated by the deformations of individual grains as well as interactions between grains (e.g., Cuffey2010). The viscosity of an individual ice crystal strongly depends on the direction of the applied strain, and it will most readily deform as shear is applied orthogonal to the c axis (crystal symmetry axis, normal to the hexagonal basal planes), leading to slip of the individual basal planes (e.g., McConnel1891; Hobbs2010; Petrenko and Whitworth2002). In polycrystalline ice subjected to strain, the crystals may undergo lattice rotation or recrystallization, both of which result in non-isotropic c-axis distributions and a bulk anisotropic viscosity (e.g., Faria et al.2014b).

In this work we only consider scenarios in which the c axes are distributed isotropically (uniform fabric), are aligned in a single direction (single-pole fabric), or lie in a plane (girdle fabric). The latter is of primary importance for the studied ice.

The crystal orientation fabric is experimentally most commonly observed through the use of polarized-light microscopy on thin sections of ice core samples (e.g., Alley1988; Wilson et al.2003; Langway1958; Wilen et al.2003).

The average crystal size and elongation can also be quantified directly through microscopy (e.g., Fitzpatrick et al.2014).

While ice core analysis uniquely delivers ground-truth information, it is limited by its small sampling volume and often unable to resolve the absolute direction of fabric orientation as the core orientation is not preserved in the drilling process (e.g., Westhoff et al.2021). Volumetric quantities such as grain volumes and shapes are generally not directly accessible through the commonly employed techniques. Grain sizes and elongations evaluated through the microscopy of thin slices cut from ice cores in turn often depend on the sample plane.

Ice fabric can be imaged not only in ice cores. It also leads to a directionality in the propagation of sound and electromagnetic radiation. The mechanical anisotropy of ice results in a fabric-dependent speed of sound, as has, for example, been measured using a sonic logger in boreholes (Kluskiewicz et al.2017). Ice crystals also are a birefringent material such that any incoming electromagnetic radiation is separated into an ordinary and extraordinary ray of perpendicular polarizations with respect to the c axis and which propagate with different refractive indices. Today this is primarily employed by polarimetric radar systems to infer fabric properties (e.g., Fujita et al.2006; Matsuoka et al.2003; Jordan et al.2019; Young et al.2021) through periodic power anomalies detected as a result of the direction and polarization-dependent delay in the propagation of radio waves.

Recently, as part of ice calibration measurements for the IceCube Neutrino Observatory (Aartsen et al.2017), Chirkin (2013d) described the observation of an “ice optical anisotropy”. At receivers 125 m away from isotropic 400 nm emitters, about twice as much light is observed for emitter–receiver pairs oriented along the glacial flow axis versus orthogonal to the flow axis (see Fig. 4). The effect was originally modeled as a direction-dependent modification to impurity-induced Mie scattering quantities, either through a modification of the scattering function as proposed by Chirkin (2013d) or through the introduction of a direction-dependent absorption as introduced by Rongen (2019). As also shown by Rongen (2019), both parameterizations lack a thorough theoretical justification and resulted in an incomplete description of the IceCube data (see Fig. 6).

First attempts to attribute the observed effect not to Mie scattering but to the ice-intrinsic birefringence have been made by Chirkin and Rongen (2020). Here the optical anisotropy results from the cumulative diffusion that a beam of light experiences as it is refracted or reflected on many grain boundaries in a birefringent polycrystal with a preferential c-axis distribution. The wavelength of ∼400 nm employed in the IceCube calibration studies is significantly smaller than the average grain size, which is expected to be on the millimeter scale. Thus, the spacing of grain boundaries and the distribution of encountered grain boundary orientations, both of which are a function of the average grain shape, must be accounted for in addition to the fabric.

In this scenario the diffusion is found to be strongest when photons initially propagate along the ice flow axis and smallest when initially propagating orthogonal to the flow axis. In addition photons are, on average, deflected towards the flow axis. The deflection per unit distance increases for stronger girdle fabrics, a larger average crystal elongation, or a smaller average crystal size. For crystal configurations and/or realizations where the deflection outweighs the additional diffusion along the flow axis compared to the diffusion along the orthogonal direction, the photon flux along the flow axis will increase with distance compared to the photon flux along the orthogonal axis. This interplay between diffusion and deflection leaves a unique imprint in the spatial and temporal light signatures recorded by IceCube. Due to computational limitations, a grain-resolving anisotropic optical model has been parameterized by Abbasi et al. (2021) using diffusion functions. These functions in turn have been applied as an extension to the existing homogeneous ice optical simulation. These simulations, assuming different ice crystal realizations, have then been compared to LED flasher data, which partially constrain the crystal fabric, size, and elongation.

Work on this model has so far been performed and proceedings published (Chirkin2013d; Chirkin and Rongen2020; Abbasi et al.2021) in the context of detector calibration for the measurements performed by IceCube. With this paper we will summarize the full extent of past and ongoing modeling of the ice optical anisotropy to a geophysical audience for the first time. The described measurements may be unique to IceCube and thus not easily adopted as a tool in glaciology. Nevertheless we believe that they yield an interesting complementary view on ice physical properties and through comparison to ice core data, in particular from SPC14 (Casey et al.2014) drilled ∼1 km from the IceCube array, will be informative for the modeling of ice dynamics.

This paper has the following structure: Sect. 2 introduces the IceCube Neutrino Observatory (Sects. 2.1 and 2.2) and how it employs ice as a detection medium (Sect. 2.3). Section 3 describes the properties of the LED calibration data used in this study (Sect. 3.1), explains the photon propagation software used to generate simulated data (Sect. 3.2), and details the likelihood analysis comparing simulated to experimental data in order to infer ice properties (Sect. 3.3). The state of the isotropic, layered model used to describe the ice optical properties prior to the discovery of the ice optical anisotropy is briefly reviewed in Sect. 3.4. The experimental signature of the ice optical anisotropy (Sect. 4.1) and early modeling attempts (Sect. 4.3) are summarized in Sect. 4. The newly developed model to account for the ice optical anisotropy based on the ice-intrinsic birefringence is described starting with Sect. 5. Section 5.2 explains the electromagnetic theory governing the birefringence in polycrystals, while Sect. 5.3 introduces a software package to simulate the resulting diffusion patterns. Section 5.4 compares the experimental signatures and conceptual understanding of the underlying optics to birefringence observations in radar sounding, a field most readers are probably more familiar with. Section 6 explains how the diffusion patterns are applied in the IceCube photon propagation simulation (Sects. 6.1 and 6.2) and how crystal properties have been inferred (Sect. 6.3). The resulting ice optical model is described in Sect. 7. Section 8 discusses shortcomings of the model as well as future measurements in upcoming IceCube extensions and through drill-hole logging.

2 The IceCube Neutrino Observatory

2.1 Scientific context: neutrino astronomy

The IceCube Neutrino Observatory operates within the context of astroparticle physics and multi-messenger astronomy. While astronomy is most commonly associated with the observation of the universe in visible light, today the entire electromagnetic spectrum ranging from radio waves to hard X-rays and ultrahigh-energy gamma rays is exploited, with each spectral range giving a complementary insight. Infrared radiation, for example, is only weakly attenuated by interstellar dust (e.g., Li and Draine2001), allowing for the imaging of objects otherwise obscured by dust clouds.

In addition to photons, the quanta of light, other stable messenger particles are also observed. Most prominently cosmic rays, primarily protons, have now been found at energies exceeding 5×1019 eV, the equivalent of roughly 8 J per particle (Aab et al.2020). While these ultrahigh-energy cosmic rays offer the promise to probe the highest-energy processes in the universe, they are deflected by magnetic fields along their journey from source to detection (e.g., Aartsen et al.2015). Thus their arrival directions at Earth cannot be easily traced to their origins, making the identification of the sources of high-energy cosmic rays one of the biggest challenges in astroparticle physics.

Associated with the production of high-energy protons, one expects the production of high-energy neutrinos, or astrophysical neutrinos (e.g., Margolis et al.1978). These are electrically neutral, elementary particles belonging to the family of leptons (as counterparts to electrons, muons, and taus). As they are electrically neutral, they are not deflected in magnetic fields and thus point back to their point of origin. Additionally they only interact through the weak force and as a result can traverse vast astronomical distances without their flux being significantly attenuated. While these properties ensure that neutrinos carry unbiased information about the highest-energy regions of the universe, these same properties also make them exceptionally hard to detect, requiring cubic-kilometer-scale detectors to intercept a few dozen astrophysical neutrinos per year (e.g., Markov1960). Detectors of this scale can only be built into natural media such as ocean water or glacial ice, which need to be characterized in situ as presented here.

https://tc.copernicus.org/articles/18/75/2024/tc-18-75-2024-f01

Figure 1Overview of the IceCube detector. The 86 cables of the deep in-ice array, called strings, are indicated as gray lines, with black dots for the 60 DOMs per string. The lateral spacing of the strings is ∼125 m with a vertical spacing between DOMs of 17 m. A central part of the array, called DeepCore, is more densely instrumented. The detector is capped by a surface detector, aimed at cosmic ray physics, called IceTop. Figure credit: IceCube.

2.2 The detector

The IceCube Neutrino Observatory (Aartsen et al.2017) has, among other science goals, been built to explore the cosmos using high-energy astrophysical neutrinos. Located about 1 km from the geographic South Pole, it is logistically supported by the Amundsen–Scott South Pole Station. IceCube features a surface detector, called IceTop, as well as the deep in-ice array of interest here, which consists of 5160 optical sensors instrumenting a 1 km3 volume of ice at depths of 1450 m to 2450 m. The instrumentation layout is shown in Fig. 1. Each sensor, called a digital optical module (DOM) (Stokstad2005; Abbasi et al.2010, 2009), is equipped with a 10 in. photomultiplier tube sensitive to light between approximately 300–600 nm and all required readout electronics to be able to time-stamp the arrival time of individual photons to within 2 ns. It addition each DOM features 12 LEDs which can emit light pulses of known intensity and duration into the ice and which are used to calibrate the optical properties of the instrumented ice, as detailed in this paper. Construction took 6 years, with 86 holes of 60 cm diameter being drilled using hot-water drilling (Benson et al.2014). Cables called “strings” were instrumented with 60 DOMs each and deployed in the boreholes.

The top 1450 m was left without instrumentation because of the strongly scattering ice that exists in this region. The depth where most bubbles have converted to air hydrates was determined to be roughly 1350 m by the predecessor experiment AMANDA (Ackermann et al.2006).

Upon a neutrino interaction in the ice, charged particles with relativistic velocities are created, which emit blue light along their path through a process called Cherenkov radiation (Cherenkov1937). A small fraction of this light, after propagating through the ice, reaches some of the sensors and is detected. Reconstruction of the particle properties, namely energy and direction, relies on a precise understanding of the optical properties of the instrumented ice. Generally the particle energy is proportional to the amount of detected light, while the arrival direction is inferred from the geometric deposition of the light as well as its timing information (Aartsen et al.2013d).

Since its completion in 2010, the IceCube detector has been in continuous operation with an up-time exceeding 99 %. On average around 2000 particle events are detected and reconstructed per second, with the vast majority of these being particle showers induced by cosmic rays striking Earth's atmosphere and only a vanishing fraction (approximately hundreds per year) being astrophysical neutrinos. Using IceCube data, a wide range of results have been obtained. Those include, among others, the discovery of a high-energy astrophysical neutrino flux (Aartsen et al.2014) and first associations of high-energy neutrinos with astrophysical objects (Aartsen et al.2018a), competitive measurements of neutrino oscillation parameters (Aartsen et al.2018b), and world-leading limits on possible dark-matter properties (Albert et al.2020).

2.3 Glacial ice as an optical medium

IceCube detects individual photons that are produced through Cherenkov radiation or as emitted by the calibration LEDs. On their way from their source to a potential detection at a DOM these photons are subject to absorption and scattering in the ice, shaping both the intensity pattern in the detector and the arrival time distributions on every module.

Absorption is characterized by a wavelength-dependent λ absorption length λa(λ), the propagated distance at which the survival probability of a photon drops to 1/e. In contrast, scattering does not reduce the photon count but results in discrete direction changes at an average distance of λb(λ), the geometric scattering length. Scattering is further described by the scattering function, a probability density distribution describing the probability of deflection angles in each scattering process. Neglecting its functional form, the scattering function is described through the average deflection angle or asymmetry parameter g=〈cos θ. The effective scattering length λeff, denoting the distance at which an initially directional beam becomes diffuse independent of the scattering function, is given as (Aartsen et al.2013c)

(1) λ eff ( λ ) = λ b ( λ ) / 1 - g ( λ ) .

As pure ice itself is only very weakly absorbing (Warren and Brandt2008) (and as we will see later also effectively weakly scattering), the light propagation is dominated by Mie scattering on impurities. In this scenario absorption and scattering strengths are commonly denoted by coefficients (a=1/λa and be=1/λeff), which are proportional to the impurity concentration (Ackermann et al.2006). The primary impurity constituents contributing to absorption and scattering were identified by He and Price (1998) to be mineral dust, marine salt, and acid droplets as well as soot. These constituents range from nanometer to micrometer in size, with their combined size distribution resulting in a very strong forward scattering with g≈0.95 at the relevant wavelengths around 400 nm (He and Price1998). The impurities have been deposited with the snow precipitation over the past 100 kyr, which was compressed into the ice that is present today at the relevant depths. The impurity composition and concentration, and thus also the optical properties, accordingly trace the global climatological conditions such as dust and aerosols in the atmosphere in the past. This stratigraphy was traced at millimeter resolution using a laser dust logger deployed down seven IceCube drill holes as described by Aartsen et al. (2013a).

While not contributing to absorption, air hydrates also contribute to scattering. Their number density is large, and their large size (Uchida et al.2011), compared to the typical wavelengths considered, results in isotropic scattering. Yet due to the small difference in refractive index (Uchida et al.1995) compared to ice they contribute at most a few percent to the overall scattering coefficient (He and Price1998). Thus, scattering on air hydrates was previously not modeled explicitly and its impact was effectively incorporated into the overall scattering coefficients. Diffusion through scattering on grain boundaries was also already quantitatively estimated by He and Price (1998) to contribute about as much as air hydrates to the overall scattering coefficient. At the time the average deflection process described in this work was not known and thus its large importance not realized. The quantitative contribution of diffusion in the polycrystal to the overall scattering coefficient as derived in this work is given in Sect. 7.

Describing the depth dependence

The detailed stratigraphy associated with the yearly layering cannot be constrained through IceCube data, nor is it needed in order to accurately describe the photon propagation over large distances exceeding tens of meters. Instead, average properties in 10 m depth increments, here called “ice layers”, are being considered. The absolute depths of these layers, such as shown in Fig. 3, are referenced to a location in the center of the surface area of the detector. At any other location in the detector the same layers are found at slightly different depths following the layer undulations as will be described in Sect. 3.4.1. Each layer is described by its dust-induced absorption and scattering coefficients at a wavelength of 400 nm. These are scaled to other wavelengths as described by Aartsen et al. (2013c).

While all parameters are in principle depth-dependent, e.g., the asymmetry factor g due to changes in the impurity composition, some are deemed constant enough to be described by a single global value or functional parameterization. These are the coefficients describing the wavelength dependence and the parameterization of the scattering function, achieved through a mixture of the Henyey–Greenstein (Henyey and Greenstein1941) and simplified Liu (Liu1994) approximations of Mie scattering, as well as its asymmetry parameters g. Thus six global parameters (three for the wavelength dependencies, one ice intrinsic absorption in the infrared, two for the scattering function – g and the mixing ratio) and about 100 layers within the instrumented volume, with individual dust-induced absorption and scattering coefficients each, are required to describe the layered ice properties.

3 Deriving ice optical properties from LED calibration data

3.1 LED calibration data

As will be described in Sect. 3.4 the absorption and effective scattering lengths encountered at IceCube depths range up to 400 and 100 m, respectively. The limited volume of the ice cores thus does not allow for a direct measurement of optical properties, even though they are able to provide information on the impurity constituents and their size distributions. To enable in situ calibration of the ice optical properties, each of the 5160 DOMs deployed in ice is equipped with 12 light-emitting diodes (LEDs) that are positioned on a “flasher” board and can emit light one at a time or in simultaneous combinations. The LEDs are placed in pairs at 60 increments in azimuth, with one LED at a 48 elevation angle and the other pointing horizontally into the ice. Most of the LEDs emit light centered at the 405 nm wavelength in a cone of about 9.7 width (root mean square – rms). The duration and intensity of the light flashes can be configured and range between 6 and 70 ns (full width at half maximum – FWHM) and up to 1.2×1010 photons per flash.

For this study data with all available LEDs flashing individually and at the highest possible intensity have been used. Upon an LED flash the arrival times of photons received in all other DOMs are recorded. An example light curve, histogramming the measured arrival times, for one emitter–receiver pair is shown in Fig. 2.

3.2 Photon propagation simulation

From the recorded LED data, ice properties are inferred by comparing the data to an expectation given different hypothesized optical properties and ice crystal orientations. For a point-like emitter in the far field (dλeff) and given a weak absorption coefficient compared to the scattering coefficient, the arrival time distribution u(t), which is the density function belonging to the light curves at a distance d from an isotropic source, is described by a Green's function (Ackermann et al.2006) as

(2) u ( d , t ) = 1 ( 4 π D t ) 3 / 2 exp - d 2 4 D t exp - t c ice λ a ,

where D=ciceλeff/3 is the diffusion constant. As evident from this equation, the time of the rising edge is generally sensitive to the scattering coefficient, while the slope of the tail is determined by the absorption coefficient. While this behavior is generally also observed outside the far field, the Green's function is inaccurate in the semi-diffuse regime given by the clean, layered ice and at the sensor spacings used in IceCube. Thus, the photon propagation needs to be fully modeled in simulation. This is achieved through the use of photon propagation software, namely the photon propagation code (PPC) (Chirkin2013a).

PPC aims to be a full first-principles simulation, tracking each photon individually and as accurately as possible. For every created photon the total lifetime, or absorption weight in multiples of absorption lengths, is sampled from an exponential distribution with unity scale. Next the distance to the next scattering process is determined in the same fashion and the photon is moved through a depth-layered ice model along its current propagation direction towards the next scattering center. For each layer traversed, the length multiplied by the local absorption and scattering coefficient is subtracted from the current absorption and scattering weight. When the scattering weight reaches zero, the scattering site has been reached and the photon is deflected according to the modeled scattering function. The scattering transport process is repeated until the photon is either absorbed, as the absorption weight reaches an epsilon cut-off value, or the photon is incident on a DOM and stored for later processing.

PPC has been in active development and use since 2009. As photons propagate independently of each other, their simulation is an ideal use case for parallelization using Graphics Processing Units (GPUs). Using a single GPU, the full paths of ∼108 photons can be simulated per second, corresponding to simulating one full LED flash in 100 s. Computational resources are still the limiting factor in these studies, in particular when it comes to evaluating systematic uncertainties through repeated analysis under slightly perturbed assumptions. For this study simulations amounting to roughly 400 000 GPU hours have been performed on the IceCube computing cluster.

3.3 Likelihood analysis

The photon propagation described in the previous section enables reproducing (LED) events in simulation, given a set of model parameters including a realization of the ice properties. Most ice calibration studies perform an optimization of the ice assumptions by minimizing the discrepancies between simulated and measured LED events. In practice, the best estimators for the ice properties are obtained through a log-likelihood minimization, where a single likelihood value is computed for every pair of emitter and receiver DOMs. For this purpose, the experimental and simulated events are averaged over the number of repetitions in this LED configuration (usually around 200 in data and 10 in simulation). The light curve of each receiving DOM is then binned in time using a Bayesian blocking (Scargle1998) algorithm, where each bin is multiples of 25 ns long and balances maximizing photon statistics per bin with accurately describing the rate of change in photon counts at the rising and trailing edges.

Example flasher light curve
https://tc.copernicus.org/articles/18/75/2024/tc-18-75-2024-f02

Figure 2Example flasher light curve in 25 ns binning. DOM 50 on string 1 emits light which is detected by DOM 55 on string 8 about 150 m away. The data are averaged over 240 repetitions. The simulation is averaged over 10 repetitions.

Download

The per-event average expectation in each bin is a function of the sampled ice properties and nuisance parameters, such as a per-LED light yield, a timing offset of the light emission with regard to the LED trigger, and the absolute LED orientations. The likelihood function used for comparing this expectation to data is given by

(3) - ln L = i s i ln s i / n s μ s i + d i ln d i / n d μ d i + 1 2 σ 2 ln μ d i μ s i 2 ,

where i denotes a receiver DOM and time bin of its light curve, si and di the photon count in simulation and data for this bin, respectively, ns and nd the simulation repetitions and number of data events, σ the model error, and μs and μd the simulation and data expectation values; −ln ℒ is abbreviated as LLH in the following.

The model error takes into account potential discrepancies in reproducing data with a simulation that may be incomplete or may use nonideal parameterizations. Using the model error, it is assumed that a difference between the expectation values of simulations and data can exist even at the best-fit point, μsμd(si+di)/(ns+nd). This is modeled through the penalty term in the likelihood (Chirkin2013b). This extension also requires an optimization of the now in principle independent expectation values within the likelihood calculation and is performed as described by Chirkin (2013b).

This likelihood (Chirkin2013b) improves on a common Poisson likelihood by taking into account the uncertainty of the expectation caused by the small statistics of the simulated data compared to the experimental data. Therefore, the expectation is optimized including the knowledge of the limited statistics of both the simulated and experimental data. In the limit of infinite statistics of simulated data this likelihood converges to a saturated Poisson likelihood.

The parameters of the ice model are generally obtained through likelihood scans, where each scan point is one realization of the ice model parameters tested against flasher data. The timing offset and LED intensity nuisance parameters are optimized for each realization analytically and through a number of low statistics iterations.

The likelihood method described above does not, in general, fulfill Wilks' theorem (Wilks1938), which would, under certain conditions, allow one to approximate the distribution of the likelihood ratio between the best-fit and null hypotheses with a chi-squared distribution. As such, the log-likelihood contour of a one-dimensional likelihood scan enclosing the minimum by ΔLLH of 1 does not represent a 1σ statistical uncertainty. Instead, the spread in LLH values equivalent to the 1σ uncertainty is obtained by re-simulating a realization close to the optimum a number of times and computing the standard deviation of the resulting LLH values.

Fitting the flasher data, the statistical errors in the ice properties, in particular the layered absorption and scattering coefficients, are entirely due to the limited simulation statistics but generally remain below 1 %. Thus the statistical error is subdominant compared to systematic biases introduced through incomplete modeling. This bias is hard to quantify, in particular due to the enormous computational cost. Taking into account the limited knowledge of the relative detection efficiencies of the DOMs, the discrepancy between fitted values using only horizontal or only tilted LEDs and different realizations of the modeled scattering function, the systematic uncertainty on the scale of absorption and scattering coefficients is estimated to be around 5 %.

https://tc.copernicus.org/articles/18/75/2024/tc-18-75-2024-f03

Figure 3Stratigraphy of fitted absorption and scattering strength. Properties above the detector (<1450 m) are taken from AMANDA measurements (Ackermann et al.2006) or are extrapolated from dust logger data. Properties below (>2450 m) are extrapolated using the stratigraphy as obtained from the EDML ice core (Bay et al.2010) and ice age vs. depth curve from Price et al. (2000).

Download

3.4 The South Pole Ice Model (SPICE)

Employing the experimental and analysis methods described above, absolute absorption and scattering coefficients and their wavelength scaling have been measured for all instrumented depths as described in detail by Ackermann et al. (2006) and Aartsen et al. (2013c). The resulting model, called the South Pole Ice Model (SPICE), continues to be updated and refined as new aspects of the instrumentation such as the properties of the refrozen drill columns (Chirkin et al.2021) as well as previously unconsidered features in the ice begin to be modeled. The stratigraphy, used as the starting point for this study, is shown in Fig. 3. At the instrumented depths, absorption lengths mostly exceed 100 m, with the most significant exception being a region at around 2000 m, in IceCube commonly referred to as “the dust layer”. This has been associated with a period of continuously elevated dust concentrations during a stadial around 65 000 years ago (Ackermann et al.2006).

While primarily developed and employed for the simulation of particle interactions, the deduced model parameters are also informative of ice properties in general. Most prominently the lowest measured absorption coefficients now serve as a reference for an upper limit on ice-intrinsic absorption as compiled by Warren and Brandt (2008). The technique of time-resolved photon counting has recently also been adopted by Allgaier et al. (2022) to deduce impurity concentrations in firn.

3.4.1 Layer undulation

One relevant complication is the undulations of layers of equal optical properties within the instrumented volume. As established from ground-penetrating radar sounding (e.g., Fujita et al.1999) ice isochrons can be traced over thousands of kilometers. While the ice surface is generally flat, deeper layers tend to gradually follow the topography of the underlying bedrock, with additional features such as upwarping and folds in basal ice (e.g., Cooper et al.2019; Dow et al.2018; MacGregor et al.2015).

Available radar data generally do not have the spatial resolution required to map features within the instrumented volume of IceCube. Instead the depth offset of characteristic features as observed in the dust logger data from seven different IceCube holes has been used to interpolate the depth-dependent layer undulations assuming an undisturbed chronological layering as described by Aartsen et al. (2013a). Layers with roughly constant scattering and absorption change in depth by as much as 60 m as one moves across the ∼1 km detector. This gradient is mainly found along the SW direction, orthogonal to the flow direction. At the location of IceCube, the ice flows in the direction grid NW at a rate of about 10 m yr−1 (Lilien et al.2018), slowly draining into the Weddell Sea after flowing through the Pensacola–Pole Basin (Paxman et al.2019).

Within the context of the ice model, the depth offset at which a given ice layer is encountered relative to the stratigraphy as defined in the center of the detector is generally referred to as “tilt”. The orientation of the main gradient is termed “tilt direction”. Within the context of this work the tilt model as described by Aartsen et al. (2013c) is employed.

https://tc.copernicus.org/articles/18/75/2024/tc-18-75-2024-f04

Figure 4Ice optical anisotropy seen as azimuth-dependent intensity excess in flasher data. Each dot is the observed intensity ratio for one pair of light-emitting and light-detecting DOMs comparing data to a simulation with no anisotropy modeling enabled. The tilt and flow directions are shown for reference.

Download

4 The ice optical anisotropy

4.1 Experimental signature

Given the optical modeling discussed so far, the amount of light received from an isotropic source should not depend on the direction of the receiver with respect to the emitter. However, if we consider many DOMs, each with their 12 calibration LEDs and at a random azimuthal orientations in the refrozen drill holes, as isotropic emitters and we average observations along different directions of emitter–receiver pairs of DOMs, we find a significant directional dependence. About twice as much light is observed along the direction of the ice flow compared to the orthogonal ice tilt direction when measured at distances of ∼125 m, as se