Valency
alpha
PDF

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

Rasha AbbasiMarkus AckermannJenni AdamsNakul AggarwalJuanan Aguilar
Fresh — published within the last 30 daysv1 publishedCC BY 4.063 min read
AI contributions — not recorded for this paperView details

Abstract

The IceCube Neutrino Observatory instruments about 1 km 3 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.

Published by Copernicus Publications on behalf of the European Geosciences Union. Cristian Jesus Lozano Mariscal42^{42}, Lu Lu39^{39}★, Francesco Lucarelli28^{28}★, Andrew Ludwig24,35^{24,35}★, William Luszczak39^{39}★, Yang Lyu8,9^{8,9}★, Wing Yan Ma63^{63}★, Jim Madsen39^{39}★, Kendall Mahn24^{24}★, Yuya Makino39^{39}★, Sarah Mancina39^{39}★, Wenceslas Marie Sainte39^{39}★, Ioana Mariş12^{12}★, Szabolcs Marka45^{45}★, Zsuzsa Marka45^{45}★, Matthew Marsee58^{58}★, Ivan Martinez-Soler14^{14}★, Reina Maruyama44^{44}★, Thomas McElroy25^{25}★, Frank McNally37^{37}★, James Vincent Mead22^{22}★, Kevin Meagher39^{39}★, Sarah Mechbal63^{63}★, Andres Medina21^{21}★, Maximilian Meier16^{16}★, Stephan Meighen-Berger27^{27}★, Yarno Merckx13^{13}★, Jessie Micallef24^{24}★, Daniela Mockler12^{12}★, Teresa Montaruli28^{28}★, Roger Moore25^{25}★, Bob Morse39^{39}★, Marjon Moulai39^{39}★, Tista Mukherjee31^{31}★, Richard Naab63^{63}★, Ryo Nagai16^{16}★, Uwe Naumann62^{62}★, Amid Nayerhoda47^{47}★, Jannis Necker63^{63}★, Miriam Neumann42^{42}★, Hans Niederhausen24^{24}★, Mehr Nisa24^{24}★, Sarah Nowicki24^{24}★, Anna Obertacke Pollmann62^{62}★, Marie Oehler31^{31}★, Bob Oeyen29^{29}★, Alex Olivas19^{19}★, Rasmus Orsoe27^{27}★, Jesse Osborn39^{39}★, Erin O’Sullivan61^{61}★, Hershal Pandya43^{43}★, Daria Pankova60^{60}★, Nahee Park33^{33}★, Grant Parker4^{4}★, Ek Narayan Paudel43^{43}★, Larissa Paul41^{41}★, Carlos Pérez de los Heros61^{61}★, Lilly Peters1^{1}★, Josh Peterson39^{39}★, Saskia Philippen1^{1}★, Sarah Pieper62^{62}★, Alex Pizzuto39^{39}★, Matthias Plum49^{49}★, Yuiry Popovych40^{40}★, Alessio Porcelli29^{29}★, Maria Prado Rodriguez39^{39}★, Brandon Pries24^{24}★, Rachel Procter-Murphy19^{19}★, Gerald Przybylski9^{9}★, Christoph Raab12^{12}★, John Rack-Helleis40^{40}★, Mohamed Rameez22^{22}★, Katherine Rawlins3^{3}★, Zoe Rechav39^{39}★, Abdul Rehman43^{43}★, Patrick Reichherzer11^{11}★, Giovanni Renzi12^{12}★, Elisa Resconi27^{27}★, Simeon Reusch63^{63}★, Wolfgang Rhode23^{23}★, Mike Richman48^{48}★, Benedikt Riedel39^{39}★, Ella Roberts2^{2}★, Sally Robertson8,9^{8,9}★, Steven Rodan55^{55}★, Gerrit Roellinghoff55^{55}★, Martin Rongen26,40^{26,40}★, Carsten Rott52,55^{52,55}★, Tim Ruhe23^{23}★, Li Ruohan27^{27}★, Dirk Ryckbosch29^{29}★, Devyn Rysewyk Cantu24^{24}★, Ibrahim Safa14,39^{14,39}★, Julian Saffer32^{32}★, Daniel Salazar-Gallegos24^{24}★, Pranav Sampathkumar31^{31}★, Sebastian Sanchez Herrera24^{24}★, Alexander Sandrock23^{23}★, Marcos Santander58^{58}★, Sourav Sarkar25^{25}★, Subir Sarkar46^{46}, Merlin Schaufel1^{1}★, Harald Schieler31^{31}★, Sebastian Schindler26^{26}★, Berit Schlüter42^{42}★, Torsten Schmidt19^{19}★, Judith Schneider26^{26}★, Frank Schröder31,43^{31,43}★, Lisa Schumacher27^{27}★, Georg Schwefer1^{1}★, Steve Sclafani48^{48}★, Dave Seckel43^{43}★, Surujhdeo Seunarine50^{50}★, Ankur Sharma61^{61}★, Shefali Shefali32^{32}★, Nobuhiro Shimizu16^{16}★, Manuel Silva39^{39}★, Barbara Skrzypek14^{14}★, Ben Smithers4^{4}★, Robert Snihur39^{39}★, Jan Soedingrekso23^{23}★, Andreas Søgaard22^{22}★, Dennis Soldin32^{32}★, Christian Spannfellner27^{27}★, Glenn Spiczak50^{50}★, Christian Spiering63^{63}★, Michael Stamatikos21^{21}★, Todor Stanev43^{43}★, Robert Stein63^{63}★, Thorsten Stelzberger9^{9}★, Timo Stürwald62^{62}★, Thomas Stuttard22^{22}★, Greg Sullivan19^{19}★, Ignacio Taboada6^{6}★, Samvel Ter-Antonyan7^{7}★, Will Thompson14^{14}★, Jessie Thwaites39^{39}★, Serap Tilav43^{43}★, Kirsten Tollefson24^{24}★, Christoph Tönnis56^{56}★, Simona Toscano12^{12}★, Delia Tosi39^{39}★, Alexander Trettin63^{63}★, Chun Fai Tung6^{6}★, Roxanne Turcotte31^{31}★, Jean Pierre Twagirayezu24^{24}★, Bunheng Ty39^{39}★, Martin Unland Elorrieta42^{42}★, Karriem Upshaw7^{7}★, Nora Valtonen-Mattila61^{61}★, Justin Vandenbroucke39^{39}★, Nick van Eindhoven13^{13}★, David Vannerom15^{15}★, Jakob van Santen63^{63}★, Javi Vara42^{42}★, Joshua Veitch-Michaelis9^{9}★, Stef Verpoest29^{29}★, Doga Veske45^{45}★, Christian Walck53^{53}★, Winnie Wang39^{39}★, Timothy Blake Watson4^{4}★, Chris Weaver24^{24}★, Philip Weigel15^{15}★, Andreas Weindl31^{31}★, Jan Weldert40^{40}★, Chris Wendt39^{39}★, Johannes Werthebach23^{23}★, Mark Weyrauch31^{31}★, Nathan Whitehorn24,35^{24,35}★, Christopher Wiebusch1^{1}★, Nathan Willey24^{24}★, Dawn Williams58^{58}★, Martin Wolf39^{39}★, Gerrit Wrede26^{26}★, Johan Wulff1,1^{1,1}★, Xianwu Xu7^{7}★, Juan Pablo Yanez25^{25}★, Emre Yildizci39^{39}★, Shigeru Yoshida16^{16}★, Shiqi Yu24^{24}★, Tianlu Yuan39^{39}★, Zelong Zhang54^{54}★, and Pavel Zhelnin14^{14}★

1^{1}III. Physikalisches Institut, RWTH Aachen University, 52056 Aachen, Germany

2^{2}Department of Physics, University of Adelaide, Adelaide, 5005, Australia

3^{3}Dept. of Physics and Astronomy, University of Alaska Anchorage, 3211 Providence Dr., Anchorage, AK 99508, USA

4^{4}Dept. of Physics, University of Texas at Arlington, 502 Yates St., Science Hall Rm 108, Box 19059, Arlington, TX 76019, USA

5^{5}CTPS, Clark-Atlanta University, Atlanta, GA 30314, USA

6^{6}School of Physics and Center for Relativistic Astrophysics, Georgia Institute of Technology, Atlanta, GA 30332, USA

7^{7}Dept. of Physics, Southern University, Baton Rouge, LA 70813, USA

8^{8}Dept. of Physics, University of California, Berkeley, CA 94720, USA

9^{9}Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

10^{10}Institut für Physik, Humboldt-Universität zu Berlin, 12489 Berlin, Germany

11^{11}Fakultät für Physik & Astronomie, Ruhr-Universität Bochum, 44780 Bochum, Germany

12^{12}Science Faculty CP230, Université Libre de Bruxelles, 1050 Brussels, Belgium

13^{13}Dienst ELEM, Vrije Universiteit Brussel (VUB), 1050 Brussels, Belgium {sup}14`Department of Physics and Laboratory for Particle Physics and Cosmology, Harvard University, Cambridge, MA 02138, USA

{sup}15`Dept. of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

{sup}16`Dept. of Physics and The International Center for Hadron Astrophysics, Chiba University, Chiba 263-8522, Japan

{sup}17`Department of Physics, Loyola University Chicago, Chicago, IL 60660, USA

{sup}18`Dept. of Physics and Astronomy, University of Canterbury, Private Bag 4800, Christchurch, New Zealand

{sup}19`Dept. of Physics, University of Maryland, College Park, MD 20742, USA

{sup}20`Dept. of Astronomy, Ohio State University, Columbus, OH 43210, USA

{sup}21`Dept. of Physics and Center for Cosmology and Astro-Particle Physics, Ohio State University, Columbus, OH 43210, USA

{sup}22`Niels Bohr Institute, University of Copenhagen, 2100 Copenhagen, Denmark

{sup}23`Dept. of Physics, TU Dortmund University, 44221 Dortmund, Germany

{sup}24`Dept. of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA

{sup}25`Dept. of Physics, University of Alberta, Edmonton, Alberta, Canada T6G 2E1

{sup}26`Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität Erlangen–Nürnberg, 91058 Erlangen, Germany

{sup}27`Physik-department, Technische Universität München, 85748 Garching, Germany

{sup}28`Département de physique nucléaire et corpusculaire, Université de Genève, 1211 Genève, Switzerland

{sup}29`Dept. of Physics and Astronomy, University of Gent, 9000 Gent, Belgium

{sup}30`Dept. of Physics and Astronomy, University of California, Irvine, CA 92697, USA

{sup}31`Karlsruhe Institute of Technology, Institute for Astroparticle Physics, 76021 Karlsruhe, Germany

{sup}32`Karlsruhe Institute of Technology, Institute of Experimental Particle Physics, 76021 Karlsruhe, Germany

{sup}33`Dept. of Physics, Engineering Physics, and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada

{sup}34`Dept. of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA

{sup}35`Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095, USA

{sup}36`Centre for Cosmology, Particle Physics and Phenomenology – CP3, Université catholique de Louvain, Louvain-la-Neuve, Belgium

{sup}37`Department of Physics, Mercer University, Macon, GA 31207-0001, USA

{sup}38`Dept. of Astronomy, University of Wisconsin–Madison, Madison, WI 53706, USA

{sup}39`Dept. of Physics and Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin–Madison, Madison, WI 53706, USA

{sup}40`Institute of Physics, University of Mainz, Staudinger Weg 7, 55099 Mainz, Germany

{sup}41`Department of Physics, Marquette University, Milwaukee, WI 53201, USA

{sup}42`Institut für Kernphysik, Westfälische Wilhelms-Universität Münster, 48149 Münster, Germany

{sup}43`Bartol Research Institute and Dept. of Physics and Astronomy, University of Delaware, Newark, DE 19716, USA

{sup}44`Dept. of Physics, Yale University, New Haven, CT 06520, USA

{sup}45`Columbia Astrophysics and Nevis Laboratories, Columbia University, New York, NY 10027, USA

{sup}46`Dept. of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, UK

{sup}47`Dipartimento di Fisica e Astronomia Galileo Galilei, Università Degli Studi di Padova, 35122 Padova PD, Italy

{sup}48`Dept. of Physics, Drexel University, 3141 Chestnut Street, Philadelphia, PA 19104, USA

{sup}49`Physics Department, South Dakota School of Mines and Technology, Rapid City, SD 57701, USA

{sup}50`Dept. of Physics, University of Wisconsin, River Falls, WI 54022, USA

{sup}51`Dept. of Physics and Astronomy, University of Rochester, Rochester, NY 14627, USA

{sup}52`Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA

{sup}53`Oskar Klein Centre and Dept. of Physics, Stockholm University, 10691 Stockholm, Sweden

{sup}54`Dept. of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA

{sup}55`Dept. of Physics, Sungkyunkwan University, Suwon 16419, Korea

{sup}56`Institute of Basic Science, Sungkyunkwan University, Suwon 16419, Korea

{sup}57`Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan

{sup}58`Dept. of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA

{sup}59`Dept. of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802, USA

{sup}60`Dept. of Physics, Pennsylvania State University, University Park, PA 16802, USA

{sup}61`Dept. of Physics and Astronomy, Uppsala University, Box 516, 75120 Uppsala, Sweden

{sup}62`Dept. of Physics, University of Wuppertal, 42119 Wuppertal, Germany

{sup}63`DESY, 15738 Zeuthen, Germany

{sup}a`also at: Earthquake Research Institute, University of Tokyo, Bunkyo, Tokyo 113-0032, Japan

†deceased, 20 February 2022

Introduction

The 2021 IPCC report (Intergovernmental Panel on Climate Change, 2021) 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., [32]). 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 cc axis (crystal symmetry axis, normal to the hexagonal basal planes), leading to slip of the individual basal planes (e.g., [67]; [49]; [69]). In polycrystalline ice subjected to strain, the crystals may undergo lattice rotation or recrystallization, both of which result in non-isotropic cc-axis distributions and a bulk anisotropic viscosity (e.g., [35]).

In this work we only consider scenarios in which the cc 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., [16]; [93]; [56]; [91]).

The average crystal size and elongation can also be quantified directly through microscopy (e.g., [38]).

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., [90]). 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 ([54]). 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 cc axis and which propagate with different refractive indices. Today this is primarily employed by polarimetric radar systems to infer fabric properties (e.g., [41]; [66]; [53]; [96]) 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 ([8]), [25] described the observation of an “ice optical anisotropy”. At receivers 125 m away from isotropic 400 mm 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 Figure 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 [27] or through the introduction of a direction-dependent absorption as introduced by [75]. As also shown by [75], both parameterizations lack a thorough theoretical justification and resulted in an incomplete description of the IceCube data (see Figure 6).

Ice optical anisotropy seen as azimuth-dependent intensity excess in flasher data

Figure 4. Ice 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.

Comparison of fit quality achieved with different models for the ice optical anisotropy

Figure 6. Comparison of fit quality achieved with different models for the ice optical anisotropy. Shown are photon arrival time distributions (summed counts in 25 ns time bins) for all nearest-pair emitters and receivers, roughly aligned along and perpendicular to the ice flow. As more emitter–receiver pairs are included in the perpendicular case compared to the case along the ice flow, the total photon counts are not directly comparable between the two plots and should instead be compared to the curve titled “flasher data” within each plot. While the array geometry is well aligned with the flow axis, the nearest inter-module propagation direction perpendicular to the flow is roughly 30° off. The “absorption” and “scattering” models represent ad hoc, directional modifications to Mie scattering and absorption but are unable to describe timing and intensity simultaneously. “Birefringence” refers to the microstructure-based effect introduced in this paper. A combination of the absorption and birefringence model yields the closest match to data to date.

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 cc-axis distribution. The wavelength of ∼400\sim400 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 [30] 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 particularly constrain the crystal fabric, size, and elongation.

Work on this model has so far been performed and proceedings published (Chirkin, 2013d; Chirkin and Rongen, 2020; [30]) 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 SPCI4 ([22]) drilled ∼1\sim1 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.

The IceCube Neutrino Observatory

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 Draine, 2001), 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×10195 \times10^{19} eV, the equivalent of roughly 88 J per particle ([1]). 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., [7]). 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., [64]). 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., [65]). 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.

The detector

The IceCube Neutrino Observatory ([8]) 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^3 volume of ice at depths of 1450 m to 2450 m. The instrumentation layout is shown in Figure 1. Each sensor, called a digital optical module (DOM) ([81]; [12], [11]), 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 ([20]). Cables called “strings” were instrumented with 60 DOMs each and deployed in the boreholes.

Overview of the IceCube detector

Figure 1. Overview 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\sim125 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.

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 ([14]).

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 ([23]). 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 ([5]).

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 ([6]) and first associations of high-energy neutrinos with astrophysical objects ([9]), competitive measurements of neutrino oscillation parameters ([10]), and world-leading limits on possible dark-matter properties ([15]).

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 λ\lambda absorption length λa(λ)\lambda_a(\lambda), the propagated distance at which the survival probability of a photon drops to 1/e1/e. In contrast, scattering does not reduce the photon count but results in discrete direction changes at an average distance of λb(λ)\lambda_b(\lambda), 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⁡θ⟩g = \langle\cos\theta\rangle. The effective scattering length λeff\lambda_{\mathrm{eff}}, denoting the distance at which an initially directional beam becomes diffuse independent of the scattering function, is given as [5]

λeff(λ)=λb(λ)1−g(λ).(1)\lambda_{\mathrm{eff}}(\lambda) = \frac{\lambda_b(\lambda)}{1-g(\lambda)}. \tag*{(1)}

As pure ice itself is only very weakly absorbing [87] (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/λaa = 1/\lambda_a and be=1/λeffb_e = 1/\lambda_{\mathrm{eff}}), which are proportional to the impurity concentration [14]. The primary impurity constituents contributing to absorption and scattering were identified by [46] 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.95g \approx0.95 at the relevant wavelengths around 400 nm [46]. 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 [2].

While not contributing to absorption, air hydrates also contribute to scattering. Their number density is large, and their large size [85], compared to the typical wavelengths considered, results in isotropic scattering. Yet due to the small difference in refractive index [84] compared to ice they contribute at most a few percent to the overall scattering coefficient [46]. 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 [46] 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 Figure 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 [2].

Stratigraphy of fitted absorption and scattering strength

Figure 3. Stratigraphy of fitted absorption and scattering strength. Properties above the detector (<1450 m< 1450\ \mathrm{m}) are taken from AMANDA measurements [14] or are extrapolated from dust logger data. Properties below (>2450 m> 2450\ \mathrm{m}) are extrapolated using the stratigraphy as obtained from the EDML ice core [19] and ice age vs. depth curve from Price et al. (2000).

While all parameters are in principle depth-dependent, e.g., the asymmetry factor gg 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 [48] and simplified Liu [62] approximations of Mie scattering, as well as its asymmetry parameters gg. Thus six global parameters (three for the wavelength dependencies, one ice intrinsic absorption in the infrared, two for the scattering function – gg 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.

Deriving ice optical properties from LED calibration data

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 distri-ેશન butions. 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×10101.2 \times10^{10} 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 Figure 2.

Example flasher light curve

Figure 2. Example 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.

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≫λeffd \gg\lambda_{\mathrm{eff}}) and given a weak absorption coefficient compared to the scattering coefficient, the arrival time distribution u(t)u(t), which is the density function belonging to the light curves at a distance dd from an isotropic source, is described by a Green’s function [14] as

u(d,t)=1(4π⋅Dt)3/2⋅exp⁡(−d24Dt)⋅exp⁡(−tciceλa)(2)u(d,t) = \frac{1}{(4\pi\cdot D t)^{3/2}} \cdot\exp\left(-\frac{d^2}{4Dt}\right) \cdot\exp\left(-\frac{t c_{\mathrm{ice}}}{\lambda_{\mathrm{a}}}\right) \tag*{(2)}

where D=ciceλeff/3D = c_{\mathrm{ice}}\lambda_{\mathrm{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) [24].

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\sim10^8 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.

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 [78] 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.

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

−ln⁡L=∑i[siln⁡si/nsμis+diln⁡di/ndμid+12σ2(ln⁡μidμis)2](3)-\ln\mathcal{L} = \sum_i \left[s_i \ln\frac{s_i/n_s}{\mu_i^s} + d_i \ln\frac{d_i/n_d}{\mu_i^d} + \frac{1}{2\sigma^2}\left(\ln\frac{\mu_i^d}{\mu_i^s}\right)^2\right] \tag*{(3)}

where ii denotes a receiver DOM and time bin of its light curve, sis_i and did_i the photon count in simulation and data for this bin, respectively, nsn_s and ndn_d the simulation repetitions and number of data events, σ\sigma the model error, and μs\mu_s and μd\mu_d the simulation and data expectation values; −ln⁡L-\ln\mathcal{L} 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)\mu_s \ne\mu_d \ne(s_i+d_i)/(n_s+n_d). This is modeled through the penalty term in the likelihood [25]. This extension also requires an optimization of the now in principle independent expectation values within the likelihood calculation and is performed as described by [25].

This likelihood [25] 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 [92], 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\Delta\mathrm{LLH} of 1 does not represent a 1σ1\sigma statistical uncertainty. Instead, the spread in LLH values equivalent to the 1σ1\sigma 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 %.

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 [30] 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 Figure 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 [14].

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.

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., [40]) 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., [31]; [33]; [63]).

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 [2]. Layers with roughly constant scattering and absorption change in depth by as much as 60 m60\ \mathrm{m} as one moves across the ∼1 km\sim1\ \mathrm{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−110\ \mathrm{m}\ \mathrm{yr}^{-1} [60], slowly draining into the Weddell Sea after flowing through the Pensacola–Pole Basin [68].

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 [5] is employed.

The ice optical anisotropy

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\sim125\ \mathrm{m}, as seen in Figure 4. This ice optical anisotropy was first discussed in 2013 (Chirkin, 2013d). The experimental arrival time distributions are nearly unchanged compared to a simulation expectation without anisotropy (as will be evident in Figure 6).

The anisotropy axis

A determination of the axis of the ice optical anisotropy can be achieved independent of any model assumption by fitting the phase of the sinusoidal intensity modulation as shown in Figure 4. To obtain spatial resolution, the data are binned in emitting DOMs, either within a tilt-corrected depth range or by string number. Thus the data are dominated by propagation in a given depth range or in the vicinity of a given string. Figure 5 shows the resulting anisotropy axes.

Measured anisotropy axes as a function of lateral position and depth

Figure 5. Measured anisotropy axes as a function of lateral position averaging over all depth (a) and as a function of depth averaging over all strings (b). Strings on the perimeter of the detector have been excluded, as the lack of symmetric neighbors leads to potentially biased results. The sensitivity is greatly reduced in the region of strongest scattering around 2000 m. The dashed line indicates the anisotropy angle averaged over all strings.

The anisotropy axis is seen to have constant direction throughout the entire detector and is considered constant for all following investigations. The resolution is around 1∘1^\circ everywhere, except in the strongly scattering and absorbing dust layer. Edge strings are also disregarded as the lack of symmetric neighbors potentially leads to biased results.

The absolute direction is 130∘130^\circ in the IceCube coordinate system (an azimuth of 0∘0^\circ is defined with respect to the positive xx axis in Figure 5 and runs counterclockwise), equivalent to the 40∘40^\circ W meridian in the universal polar stereographic coordinate system, and is in excellent agreement with present-day flow direction as measured using a GPS stake field by [60] et al. (2018).

As a part of the models described in Sects. 4 and 5 a possible elevation angle to the anisotropy axis has been considered. In both cases a near-constant elevation angle of 5∘5^\circ on average has been fitted. However, this fit is difficult to completely disentangle from effects that may arise as a result of mis-modeling of the layer undulations or the optical properties of the refrozen drill holes ([76]). As the resulting improvement in data–simulation agreement was seen to be small, this additional complication is not further considered here. As will be explained later, this elevation angle would directly relate to an elevation angle of the crystal orientation fabric.

Early empirical modeling

Following the paradigm that ice optical properties are driven by Mie scattering on impurities, early attempts tried to model the anisotropy through directional modifications of absorption and scattering. In the original parameterization presented by Chirkin (2013d), it was argued that due to time and space reversal symmetries the absorption length and geometric scattering length cannot be direction-dependent. Therefore the anisotropy was implemented as a modification to the scattering function, the only remaining Mie scattering parameter. This effectively results in a change in the effective scattering coefficient as a function of the propagation direction. Photons propagating along the flow axis experience less scattering than photons propagating along the tilt axis or inclined from the horizontal.

While not derived from first-principles Mie calculations, the parameterization was justified to be a plausible result of elongated impurities becoming preferentially aligned by the flow and thus introducing a direction dependence to the scattering function. While several glaciological studies ([70]; [80]; [42]) explore the shapes of impurities, elongations for different impurities are not well established, nor is there to our knowledge any evidence of elongated impurities becoming oriented with the flow. Alternatively a directionality of Mie scattering may be believed to be the result of inhomogeneous impurity distributions, with some impurity types known to preferentially aggregate on the grain boundaries [82], [34] . Yet the derivation of Mie scattering properties only depends on the volumetric particle densities and is independent of homogeneity. In the context of studying the ice optical anisotropy, [75] explicitly tested this in a number of simulated toy experiments and verified that inhomogeneous impurity distributions cannot lead to large-scale anisotropy.

An evaluation of the data–simulation agreement is shown in Figure 6. It shows summed photon arrival time distributions for all nearest emitter–receiver pairs, roughly aligned along and perpendicular to the ice flow for a variety of anisotropy models and the employed flasher data. The scattering-based anisotropy model results in more intensity being observed along the flow axis. However, substantial disagreement remains between the model and the observed data. As scattering is reduced in the flow direction light arrives earlier on average. The resulting change in the rising edge position is strongly penalized in the fit and limits the amount of intensity that can be recovered.

To reduce the shift of the rising edge, a directional modification to Mie absorption was considered as an alternative by [75]. A factor of 11 modulation of the absorption coefficient was required to fit the data, which seems unphysical. As evident from Figure 6, this model results in a delayed rising edge for propagation along the flow direction as desired and did result in an improved data description compared to the scattering-based model described earlier, but it is also unable to fully match the intensity difference to data.

To conclude, while resulting in partially successful effective descriptions, directional modifications to Mie scattering or absorption cannot reproduce observations nor are such modifications well motivated on first principles.

Light diffusion in birefringent polycrystals

The electromagnetics of uniaxial, birefringent crystals

Departing from the paradigm that optical properties are purely driven by impurities, let us consider the impact of the microstructure of the ice itself on light propagation.

Light diffusion in birefringent, polycrystalline materials was discussed as early as 1955 by [73] and Viswanathan (1955). While the literature agrees that the combined effect of ray splitting on many crystal interfaces will lead to a continuous beam diffusion, the resulting diffusion patterns remained largely unexplored. [71] already considered this average overall diffusion in the context of Cherenkov neutrino telescopes but disregarded it as sub-dominant compared to scattering on impurities.

In a homogeneous, transparent, and non-magnetic medium the relation between the electric field and the displacement field as well as the magnetic fields is given as (Landau and Lifshitz, 1960) [55].

B=H,D=εE.(4)B = H,\quad D = \varepsilon E. \tag*{(4)}

As the dielectric tensor ε\varepsilon is symmetric, one can always find a coordinate system where it is diagonal ε=diag⁡(nx2,ny2,nz2)\varepsilon= \operatorname{diag}(n_x^2,n_y^2,n_z^2), with nin_i being the refractive index along the given axis. Uniaxial crystals, such as ice in glacial environments, have two distinct refractive indices: nx=ny≡no≠nz≡nen_x = n_y \equiv n_o \ne n_z \equiv n_e. The axis of the refractive index nen_e defines the optical axis and coincides with the cc axis.

A light ray entering a uniaxial crystal is split into an ordinary wave and an extraordinary wave of orthogonal polarizations. Figure 7 visualizes the orientations of all electromagnetic vectors, the plane spanned by the optical axis cc, and the wave vector kk is highlighted in gray. The electric field vector EE and the displacement vector DD for the ordinary wave are always co-linear with each other and perpendicular to both the optical axis of the crystal and the parallel propagation vectors kk and SS. However, the electric field EE for the extraordinary wave is not, in general, perpendicular to the propagation vector kk. It lies in the plane formed by the propagation vector and the displacement vector. The electric field vectors of these waves are mutually orthogonal [97]. The energy flow is given by the Poynting vector S=c4πE×HS = \frac{c}{4\pi} E \times H. For the extraordinary wave, the Poynting vector SS is not parallel to kk.

Orientation of all electromagnetic vectors for the ordinary and extraordinary ray with respect to the crystal axis ($c$ axis)

Figure 7. Orientation of all electromagnetic vectors for the ordinary and extraordinary ray with respect to the crystal axis (cc axis). See text for a detailed explanation of this figure.

While the ordinary ray always propagates with the ordinary refractive index non_o, the refractive index of the extraordinary ray depends on the opening angle θ\theta between the optical axis and the wave vector kk (as described in a later section with (7)). The difference to non_o is largest when the optical axis and the wave vector are perpendicular. In this case the extraordinary ray propagates with the refractive index nen_e.

The birefringence strength can be expressed as

β=(neno)2−1.(5)\beta= \left(\frac{n_e}{n_o}\right)^2 - 1. \tag*{(5)}

For ice β≈2×10−3\beta\approx2 \times10^{-3} across the entire visible wavelength spectrum. Refractive indices at specific wavelengths can be found in Table 1 [69].

Wavelength λ\lambda (nm)non_onen_eβ\beta
4051.31851.32002.3×10−32.3 \times10^{-3}
4361.31611.31762.3×10−32.3 \times10^{-3}
4921.31281.31432.3×10−32.3 \times10^{-3}
5461.31051.31192.1×10−32.1 \times10^{-3}
6241.30911.31052.1×10−32.1 \times10^{-3}
6911.30671.30812.1×10−32.1 \times10^{-3}

Table 1. Refractive indices of ice taken from Petrenko and Whitworth (2002) [69].

Analytic calculation of a single grain boundary transition

Assuming an arbitrary ray incident on a plane interface, we first calculate the four possible wave vectors, the ordinary and extraordinary refracted rays, and the ordinary and extraordinary reflected rays. Given the wave vectors, the four associated Poynting vectors are calculated from the boundary conditions, yielding the energy flow and as such probable photon directions.

Wave vectors

Figure 8 shows the situation at hand: an incoming wave vector kk intersects the interface and is split into four outgoing wave vectors rr. The coordinate system can always be chosen such that the surface normal nn is along the yy axis and that the surface components of kk, and as such rr, are along the xx axis. Here we implicitly assume, as an approximation, that the boundary surface is a perfect plane infinite in its extension and, without a loss of generality, that the incoming and outgoing waves are all plane waves.

Sketch of wave vectors for the incident, reflected, and refracted rays

Figure 8. Sketch of wave vectors for the incident, reflected, and refracted rays. The surface component is identical for all rays.

Because of translational symmetry of the interface surface, the surface components of all wave vectors are identical (Landau and Lifshitz, 1960): kx=rxk_x = r_x. As the wave number is given by k=2πλk = \frac{2\pi}{\lambda}, we can define a vector n\mathbf{n} such that k=ωn/c\mathbf{k} = \omega\mathbf{n}/c, whose magnitude nn is the direction-dependent refractive index n=ϵ(θ)n = \sqrt{\epsilon(\theta)}. As such the magnitude of the wave vector is proportional to the refractive index and we shall simplify ∣k∣=n|\mathbf{k}| = n in the following.

Outgoing ordinary rays

Given the magnitude non_o and surface component kxk_x of the wave vector the yy component is

ry=±no2−kx2.(6)r_y = \pm\sqrt{n_o^2-k_x^2}. \tag*{(6)}

The outgoing ordinary ray of an inbound ordinary ray is not deflected, as it does not see a change in refractive index. In the case of no birefringence, one obtains Snell’s law for refraction and the usual law for reflection (ry=−kyr_y = -k_y).

Outgoing extraordinary ray

Determining ryr_y for the extraordinary rays follows the same logic, only with a refractive index which depends on the opening angle θ\theta between the outgoing wave vector r=(rx,ry)\mathbf{r} = (r_x,r_y) and the optical axis a=(ax,ay,az)\mathbf{a} = (a_x,a_y,a_z):

1n2=1ne2+(1no2−1ne2)⋅cos⁡2θ.(7)\frac{1}{n^2} = \frac{1}{n_e^2} + \left(\frac{1}{n_o^2} - \frac{1}{n_e^2}\right)\cdot\cos^2\theta. \tag*{(7)}

The optical axis is given by the optical axis of medium 1 for the reflected ray and of medium 2 for the refracted ray. Rewriting cos⁡(θ)\cos(\theta) as scalar product between the wave vector and the optical axis gives

1ne2+(1no2−1ne2)⋅(axrx+ayry)2n2−1n2=0.(8)\frac{1}{n_e^2}+\left(\frac{1}{n_o^2}-\frac{1}{n_e^2}\right)\cdot\frac{(a_xr_x+a_yr_y)^2}{n^2}-\frac{1}{n^2}=0. \tag*{(8)}

Here n2=r2=rx2+ry2n^2=r^2=r_x^2+r_y^2. The solution is

ry=−βaxayrx±D1+βay2,(9)r_y=\frac{-\beta a_xa_yr_x\pm\sqrt{D}}{1+\beta a_y^2}, \tag*{(9)}

with

D=(βaxayrx)2−(1+βay2)(rx2⋅(1+βax2)−ne2)(10)D=(\beta a_xa_yr_x)^2-(1+\beta a_y^2)(r_x^2\cdot(1+\beta a_x^2)-n_e^2) \tag*{(10)}
=ne2⋅(1+βay2)−rx2⋅(1+β⋅(ax2+ay2)).(11)=n_e^2\cdot(1+\beta a_y^2)-r_x^2\cdot(1+\beta\cdot(a_x^2+a_y^2)). \tag*{(11)}

Of the two solutions the direction appropriate for the reflected or refracted ray is chosen and the other discarded. In the case of no birefringence (β=0\beta=0) we again obtain the solution for the ordinary ray.

Poynting vectors

Once the wave vector directions are determined, the boundary continuity conditions can be written for normal components of D\mathbf{D} and B\mathbf{B} and for tangential components of E\mathbf{E} and H\mathbf{H}. If n\mathbf{n} is a normal vector perpendicular to the interface surface, we have

n⋅D1=n⋅D2,n⋅B1=n⋅B2,(12)\mathbf{n}\cdot\mathbf{D}_1=\mathbf{n}\cdot\mathbf{D}_2,\quad\mathbf{n}\cdot\mathbf{B}_1=\mathbf{n}\cdot\mathbf{B}_2, \tag*{(12)}
n×E1=n×E2,n×H1=n×H2.(13)\mathbf{n}\times\mathbf{E}_1=\mathbf{n}\times\mathbf{E}_2,\quad\mathbf{n}\times\mathbf{H}_1=\mathbf{n}\times\mathbf{H}_2. \tag*{(13)}

Here the subscript 1 indicates the total sum of fields for incident and reflected waves, and the subscript 2 indicates the fields of the refracted waves propagating away from the boundary surface in the second medium. Since B=H\mathbf{B}=\mathbf{H} two of the equations above simply imply that B1=B2\mathbf{B}_1=\mathbf{B}_2 and H1=H2\mathbf{H}_1=\mathbf{H}_2. Together with the boundary conditions for D\mathbf{D} and E\mathbf{E}, this is a system of six linear equations. These equations are sufficient to determine the amplitudes of four outgoing waves: two reflected (ordinary and extraordinary) and two refracted (also ordinary and extraordinary). Since we only have four unknowns, two of these equations are necessarily co-linear with the rest if the wave vectors were determined correctly.

From the solution to the linear equation system the Poynting vectors and as such the photon directions of the (up to) four outgoing rays are calculated. The relative intensity of these rays, as usually denoted in Fresnel coefficients, is derived from the Poynting theorem, which for our case (no moving charges, no temporal change in total energy) is given as

∯∂VS⋅dA=0,(14)\oiint_{\partial V}\mathbf{S}\cdot\mathrm{d}\mathbf{A}=0, \tag*{(14)}

where ∂V\partial V is the boundary of a volume VV surrounding the interface. The choice of volume is arbitrary. A simple choice is a box around the interface. In the limit of an infinitely thin but wide box, it is evident that the sum of Poynting vector components normal to the interface plane is conserved.

Evanescent waves, i.e., waves with a complex wave vector, which decay away from the boundary surface and arise when the discriminant in the wave vector equation (Eq. 10) is negative, will necessarily yield vanishing contributions to such a sum. As the photon interacts with a boundary there is a brief flow of energy along the surface boundary within evanescent solutions (if any), but no energy flows away from the boundary within such solutions. The evanescent waves need to be considered when solving the boundary conditions as given in Eq. (13).

After deriving the solution presented here, we learned of the paper by Zhang and Caulfield (1996) and found that our approach is similar to the one they described.

Simulating diffusion patterns

Based on the calculations above, a photon propagation simulation for birefringent polycrystals was implemented in C++. At each grain transition, the outgoing photon is then chosen randomly, with probabilities proportional to the (up to) four normal components of the non-evanescent Poynting vectors to account for the relative intensities.

The resulting diffusion patterns, defined as the distribution of photon directions after crossing a given number of grains, depend on two factors related to the polycrystal configuration.

The assumed probability density distribution of c-axis orientations, which is the crystal orientation fabric, determines the refractive indices a photon will encounter. As measured c-axis distributions offer limited statistics and are restricted to the encountered fabric states, it is necessary here to statistically sample generic c-axis distributions. Appendix A briefly summarizes the different kinds of fabric and describes the approach developed to sample an arbitrary number of c axes based on Woodcock parameters log⁡(S1/S2)\log(S_1/S_2) and log⁡(S2/S3)\log(S_2/S_3), the usually published statistical moments associated with the fabric orientation tensor. Woodcock parameters for the ice at the South Pole are available from the South Pole Ice Core SPC14 [22], drilled by the SPICEcore project in 2014–2016 at a location ∼1\sim1 km from the IceCube array using the Intermediate Depth Drill designed and deployed by the U.S. Ice Drilling Program (IDP) [52].

It reached a final depth of 1751 m [94], which corresponds to a depth of ∼1820\sim1820 m in the IceCube ice model (see Figure 3), accounting for the layer undulation between the two reference points. The c-axis distributions have been measured by [86] at all depths and show an exceptionally clean girdle fabric at the overlapping depth as summarized in Figure A1.

Fabric versus depth trajectory with example c-axis distributions at selected depths

Figure A1. Fabric versus depth trajectory as measured in the SPC14 ice core. Individual cc-axis distributions at example depths are shown superimposed in the Schmidt equal-area projection. The ice features a prominent girdle fabric at the overlapping depths instrumented by IceCube starting at ∼1450\sim1450 m. Adapted from [86].

As evident from Snell’s law, in addition to the change in refractive index, the slope of the interface surface also dictates the refraction angle at a grain boundary transition. Thus the distribution of grain boundary plane orientations, resulting from a given grain shape, needs to be modeled in addition to the crystal orientation fabric. Appendix B shows that the surface orientation density of an ensemble of ice crystals, simulated as a polyhedral tessellation of a volume, can be approximated using a triaxial ellipsoid that represents the average shape. For a generalized ellipsoid the diffusion patterns are thus not only a function of the opening angle between the initial photon direction and the flow (as expected from the crystal orientation fabric), but also depend on the absolute zenith and azimuth orientation of the propagation direction with respect to the flow. Employing an alternative parameterization, developed prior to the one introduced in Sect. 6.1, it was determined early on that fully triaxial ellipsoids offer no advantage to describing the flasher data compared to prolate spheroids, where the major axis is aligned with the flow and the horizontal and vertical minor axes are identical. These spheroids, described by the size of the major axis and an elongation, are what we restrict ourselves to here. Grain size and shape distributions have not yet been fully published by the SPICEcore collaboration but are expected from preliminary material shown at conferences [18] as well as other cores (e.g., [89]; [61]; [82]; [35]) to be on the millimeter scale with elongations of at most a factor of 2. Both fabric and grain shape are not directly taken from ice core data but are determined from the flasher data (see Sect. 6.3).

Simulated diffusion patterns after crossing 1000 grain boundaries for four initial propagation directions relative to the flow axis and assuming on average spherical grains as well as a perfect girdle fabric are shown in Figure 9. The overall diffusion is largest when propagating along the flow direction and becomes continuously smaller towards the tilt direction. For intermediate angles the distribution is slightly asymmetric, resulting in a mean deflection towards the ice flow axis. The diffusion being largest along the flow axis results in a reduction of intensity in this direction, which is contrary to observations. The deflection, however, slowly diverts intensity from the tilt direction and overpopulates the flow (see Figure 10 direction). Thus a good fit to the data should be obtainable by finding the right combination of crystal orientation fabric, shape, and crystal size as it changes the number of crystals per distance (see Sect. 6.1).

Example diffusion patterns after photon propagation through 1000 crystals

Figure 9. Example diffusion patterns after photon propagation through 1000 crystals (roughly equivalent to 1 m) with a perfect girdle distribution of cc-axis orientations. The initial directions of the emitted photon point perpendicularly out of the picture, with an opening angle to the flow as indicated. The figures histogram the final direction vectors of 10810^{8} photons each. The change in diffusion (width of the distributions) as well as the subtle effect of photon scattering towards the ice flow (towards the right) can be seen.

Artist illustration visualizing the deflection concept

Figure 10. Artist illustration visualizing the deflection concept. Without birefringence light streams out radially from an isotropic light source. With birefringence rays get slowly deflected towards the flow axis. The effects of scattering and diffusion are not shown. The hexagonal pattern of the IceCube array around the light source is shown.

To validate our calculations and implementation, a polycrystal was realized in Zemax, a commercial optics simulation program, using a polycrystal tessellation simulated using Neper: Polycrystal Generation and Meshing [72] and exporting each interlocking monocrystal as a CAD object. The same quantitative behavior as described above is reproduced. This approach, however, does not allow for a flexible configuration and is slow to simulate reasonable photon statistics.

Comparison to fabric-induced anisotropies in radar measurements

Before incorporating the diffusion patterns into the overall IceCube simulation and fitting new ice parameters, we will discuss some conceptual differences of the birefringence-induced optical anisotropy in comparison to birefringence effects in radar measurements, which many readers may be more familiar with.

When probing the ice with radio waves the employed wavelength is orders of magnitude larger than the crystal size. Thus the waves do not interact with individual grains, and propagation is only influenced by the bulk dielectric tensor, weighting the per-crystal dielectric tensor with their relative occurrence. Since the birefringence strength β=(ne/no)2−1\beta=(n_e/n_o)^2-1 is an order of magnitude larger in radio (β∼1%\beta\sim1\%) compared to optical waves (β∼0.2%\beta\sim0.2\%), the available observables are primarily direction-dependent timing delays (either of the entire pulse or measured as a phase difference) and – for polarimetric systems – changes in the received polarization with respect to the emitted polarization.

Given the timing precision of IceCube and given the low birefringence strength in the optical regime, the effect of birefringence on timing will not be relevant here. Even assuming the unrealistic case in which one ray propagates purely with the ordinary and another with the extraordinary refractive index, the propagation delay over 250 m would only amount to ∼1\sim1 ns, which is undetectable with IceCube. Polarization is also not an available observable using IceCube data. Since each crystal effectively acts as a polarization analyzer and a large number of these are randomly sequenced, the diffusion patterns also do not depend on the initial polarization.

Instead, since the wavelength is small compared to the crystal size, light rays experience the individual grains as distinct objects and slowly diffuse through the continued refractions and reflections at the grain boundaries. Given a mean elongation or equivalently a preferential cc-axis distribution in addition to the diffusion, rays get on average slowly deflected towards the elongation axis. To our knowledge this is a newly discovered optical effect not described in the literature before.

The birefringence ice model

Parameterizing diffusion patterns

While the simulation described above in principle scales to arbitrary crystal counts, it is computationally unfeasible to explicitly simulate every single grain boundary with every simulated photon traveling dozens to hundreds of meters. For this reason, an analytic parameterization was developed, which allows describing the cumulative effect at large scales. Diffusion patterns have been simulated for a wide range of spheroid elongations (1–3) and fabric parameters (spanning the plane of Woodcock parameters between 0.1 and 4 in both dimensions). As evident from the example in Figure 9, these diffusion patterns have a strong central core with a broad large-angle tail. The tail is dominated by single large-angle reflections and as such scales linearly with the number of crystals traversed. We found that the precise simulation of the tail is unimportant, in particular as shape uncertainties of the Mie scattering function far outweigh the errors introduced by a simple parameterization. Therefore the distribution is modeled as a 2D Gaussian on a sphere, lending itself to usual scaling (with distance) relationships for mean displacement and width. The distributions are very slightly skewed towards the flow axis and are slightly better described by a skewed Gaussian. A number of more complicated functions were also fit with good success in precisely describing the underlying distribution. Figure 9 in fact uses a function with 10 parameters to illustrate all features of the distribution without statistical fluctuations. These were, however, abandoned, as no simple distance scaling could be established.

The three parameters of the diffusion pattern modeled with the 2D Gaussian on a sphere are the two widths (in the directions towards the flow, σx\sigma_x, and perpendicular to it, σy\sigma_y) and a single mean deflection towards the flow, mxm_x. The mean deflection in the perpendicular direction was zero for all cases that we chose to include in the final model (i.e., single-axis ellipsoids for particle shape and selected crystal fabric configurations). Because we mainly simulate small deflections (ignoring the long tails), we simulated the 2D Gaussian in Cartesian coordinates and then projected that to the sphere with an inverse stereographic projection. The three quantities were fitted to the following functions of angle η\eta of the initial photon direction with respect to the ice flow for simulations with a fixed number of 1000 crystal crossings.

mx=α⋅arctan⁡(δ⋅sin⁡ηcos⁡η)⋅exp⁡(−βsin⁡η+γcos⁡η)(15)m_x = \alpha\cdot\arctan(\delta\cdot\sin\eta\cos\eta) \cdot\exp(-\beta\sin\eta+ \gamma\cos\eta) \tag*{(15)}
σx=Ax⋅exp⁡(−Bx⋅[arctan⁡(Dxsin⁡η)]Cx)(16)\sigma_x = A_x \cdot\exp\left(-B_x \cdot[\arctan(D_x \sin\eta)]^{C_x}\right) \tag*{(16)}
σy=Ay⋅exp⁡(−By⋅[arctan⁡(Dysin⁡η)]Cy)(17)\sigma_y = A_y \cdot\exp\left(-B_y \cdot[\arctan(D_y \sin\eta)]^{C_y}\right) \tag*{(17)}

These functions were found to describe all considered crystal realizations with only 12 free parameters (Ax…DxA_x \ldots D_x, Ay…DyA_y \ldots D_y and α…δ\alpha\ldots\delta). Figure 11 shows the mean deflection for nine crystal configurations. Note that increasing elongation has a stronger effect compared to a strengthening fabric, i.e., increasing the value of the Woodcock parameter ln⁡(S2/S3)\ln(S_2/S_3).

Applying diffusion patterns in photon propagation

During photon propagation simulation, directions are only updated upon scattering. To minimize the additional computational burden, the new birefringence anisotropy is discretized and also evaluated only at the scattering sites. This requires scaling the diffusion, deflection, and displacement derived from simulation through 1000 grains to the number of traversed grains between two scattering sites. This introduces a new model parameter, the average grain size, and also requires taking into account the different average crystal chord lengths as a function of propagation direction (as described in [75]), further increasing the importance of elongation over fabric.

The grain size distribution, which is the size distribution of ice mono-crystals, defines the distance between interface crossings. As would be expected from a diffusion process and was confirmed in simulation, the deflection scales linearly and the diffusion scales with the square root of the number of traversed grains nn (σx,y∝n\sigma_{x,y} \propto\sqrt{n} and mx∝nm_x \propto n). The overall ice diffusion strength, including both Mie scattering and the birefringence-induced diffusion, has previously been measured to great accuracy. To decouple the fitting of anisotropy properties from this overall ice, the effective scattering Mie coefficient was reduced by the amount resulting from the birefringence-induced light diffusion assuming on average isotropic photon directions.

Updating not only a photon’s direction with deflection due to birefringence, but also the photon coordinates (as it shifts transversely with respect to straight-path expectation) at the next Mie scattering site, improves the agreement with data in the final fit. Due to the simple physics of cumulative photon deflections, the effect can be simulated at a small additional computational cost and with no additional parameters. Assuming without loss of generality that all birefringence deflections happen at constant distance interval Δl\Delta l and that these can be sampled from the same distribution (which depends on the initial photon direction), as the individual and even final calculated deflections are very small, we can express the new photon direction n\boldsymbol{n} and coordinates r\boldsymbol{r} after NN deflections as

n=n0+∑i=1NΔni,(18)\boldsymbol{n} = \boldsymbol{n}_0 + \sum_{i=1}^{N} \Delta\boldsymbol{n}_i, \tag*{(18)}
r=∑i=1Nni⋅Δl=Δl⋅N⋅n0+Δl⋅∑i=1N∑j=1iΔnj.(19)\boldsymbol{r} = \sum_{i=1}^{N} \boldsymbol{n}_i \cdot\Delta l = \Delta l \cdot N \cdot\boldsymbol{n}_0 + \Delta l \cdot\sum_{i=1}^{N} \sum_{j=1}^{i} \Delta\boldsymbol{n}_j. \tag*{(19)}

The second term in each of the two expressions above describes a cumulative direction change δn\delta\boldsymbol{n} and relative coordinate update δr\delta\boldsymbol{r}, respectively (we note that the total distance traveled is L=Δl⋅NL = \Delta l \cdot N). We can now calculate that in the limit of large NN we get

⟨δr⟩=⟨δn⟩L2,(20)\langle\delta\boldsymbol{r} \rangle= \langle\delta\boldsymbol{n} \rangle\frac{L}{2}, \tag*{(20)}
⟨Δ(δr−δnL2)2⟩=⟨Δ(δn)2⟩L212,(21)\left\langle\Delta\left(\delta\boldsymbol{r} - \delta\boldsymbol{n}\frac{L}{2}\right)^2 \right\rangle= \left\langle\Delta(\delta\boldsymbol{n})^2 \right\rangle\frac{L^2}{12}, \tag*{(21)}
⟨Δ(δr−δnL2)⋅Δ(δn)⟩=0.(22)\left\langle\Delta\left(\delta\boldsymbol{r} - \delta\boldsymbol{n}\frac{L}{2}\right) \cdot\Delta(\delta\boldsymbol{n}) \right\rangle= 0. \tag*{(22)}

Δ\Delta in the equations above is the variation (difference) from the mean of the quantity immediately following in brackets. These equations indicate that the coordinate update δr\delta\boldsymbol{r} can be sampled from a distribution with a mean given by the first equation (which could be approximated by propagating the photon half the distance with initial direction vector and the other half with the final direction vector) and variance given by the second equation. Because there is no correlation between the residual in the variance and the deflection vector, as shown by the third equation, the variance can be sampled using the already tabulated birefringence parameters independently from sampling the variance of the deflection vector.

Deflection $m_x$ as a function of opening angle to the flow for crystal configurations

Figure 11. Deflection mxm_x as described in the text as a function of opening angle to the flow for a number of crystal configurations. The black curves were fitted through the blue simulated points using the functional form introduced in (17). Note the different ordinate scales per row.

Fitting to flasher data

Besides the anisotropy direction already discussed in Sect. 4.2, the model described above requires four parameters to specify a birefringence anisotropy realization: crystal size and elongation and the two Woodcock parameters ln⁡(S1/S2)\ln(S_1/S_2) and ln⁡(S2/S3)\ln(S_2/S_3). Additionally allowing for a correction to the previously established total absorption and scattering coefficients adds two more parameters. As minimizing all six parameters for all 100 depth layers in the ice model is not computationally feasible, we need to simplify the model by identifying some parameters which are either depth-independent or have a small effect on the data–simulation agreement.

This is done through pre-fits, which either vary all parameters for a single exemplary layer or fit the depth dependence of a single parameter while keeping all other parameters fixed. The required pre-fits, as well as the final depth evaluation, were performed following the method described in [5] and summarized in Sect. 3.3. This primarily entails minimizing the summed LLH comparing the single-LED data set (where all 12 LEDs were flashed one at a time on all in-ice DOMs) with the full photon propagation simulation of these events taking into account precisely known DOM orientations as measured in [26]. Fits for individual layers were carried out by only including LEDs situated within the considered (tilt-corrected) ice layer in the LLH summation. This method offers a reduced depth resolution compared to [5] but reduces computation time while making use of the full data. An example LLH space at a depth of ∼1500\sim1500 m is shown in Figure 12. During the pre-fits the following behavior was noted: given a girdle fabric (ln⁡(S1/S2)≪ln⁡(S2/S3)\ln(S_1/S_2) \ll\ln(S_2/S_3)), the actual fabric strength has a small effect and cannot be distinguished by the data. Accordingly the fabric has been fixed to values as measured in the deepest sections of the South Pole Ice Core SPC14 ([86]) (ln⁡(S1/S2)=0.1\ln(S_1/S_2) = 0.1 and ln⁡(S2/S3)=4\ln(S_2/S_3) = 4). The fit is largely degenerate in crystal elongation and size, with small, near-spherical crystals yielding similar results to larger, more elongated realizations. Thus, the elongation was fixed to 1.4, which is a good fit at all layers and is a reasonable value given the largest value measured in the deepest parts of SPC14 (∼1.24\sim1.24) and the observed trend of increasing elongations up to that depth [17].

Log-likelihood space for layer 41 and a subset of parameters

Figure 12. Log-likelihood (LLH) space, where each point quantifies the agreement of a simulation with a set of assumed parameter values against data, for one ice layer and a subset of parameters. Each panel shows a marginalized 2D space, each point being a simulated ice realization, color-coded by its LLH distance from the best fit. In this example, the absorption anisotropy (κ1\kappa_1) is a free parameter (corresponding to the final model). This example is particularly detailed and is used to understand the behavior of the pre-fits. In particular, note the strong degeneracy in crystal elongation and size (parameterized as the scale of the major axis). Near-spherical crystals yield similar results to larger, more elongated realizations. The final fit for size, scattering, and absorption correction as performed for all layers generally contains around 100 tested realizations per layer.

Fitting the remaining parameters (absorption and scattering corrections and crystal size) for all layers yields a significant improvement as seen, for example, in the average light curves in Figure 6 (birefringence-only line). The best fit still features clearly visible discrepancies, such as an elevated intensity in the peak region in the case of propagation along the flow direction and too little intensity in the peak region in the case of propagation perpendicular to the flow direction. Problematically, the crystal sizes required to obtain this result are on the order of 0.1 mm and as such far smaller than expected from the overlapping SPC14 depths [17].

After thoroughly checking both the assumptions and implementation of the birefringence model, it was decided to reintroduce scattering as well as absorption anisotropy,

both following the formalism of [25], into the fit. As would be expected from the timing behavior, the fit does not make use of the scattering anisotropy, but surprisingly the absorption anisotropy is mixed into the birefringence model with a significant nonzero contribution. The fitted strength of the absorption anisotropy is nearly depth-independent with a directional modulation of the absorption coefficients by a factor of 2.45. This means a departure from a first-principles model but was adopted for its improvement in data–simulation agreement. After including the absorption anisotropy, absorption and scattering corrections and the crystal size were again fitted for all layers.

Resulting ice model

Figure 13 depicts the best-fit stratigraphy of grain sizes. The overall grain size of ∼\sim1 mm and the increase in size at larger depths, where ice crystals are generally larger, are as generally expected and measured in glaciology (e.g., [?]; [17]). In addition an anticorrelation between crystal size and impurity concentrations, as mapped by optical properties, can be observed. This follows the expectation that impurity-related processes such as impurity drag hinder grain growth (e.g., [34]). As noted previously, the fit is largely degenerate in elongation and size. As a result the overall size scale is somewhat unconstrained. Repeating the fit under the assumption of an elongation of 1.7 instead of 1.4, for example, results in 26 % larger circle-equivalent diameters on average.

Best-fit crystal sizes as deduced in this analysis

Figure 13. Best-fit crystal sizes as deduced in this analysis. The sphere-equivalent diameter denotes the diameter of a sphere with volume equivalent to the fitted spheroid describing the average crystal size and elongation at each depth. Error bars denote the statistical uncertainty only.

Averaged over all instrumented depths, light diffusion in the birefringent ice polycrystal amounts to an effective scattering coefficient of 2.47×10−2 m−12.47 \times10^{-2}\ \mathrm{m}^{-1}, accounting for ∼8.5\sim8.5 % of the total scattering present in the ice on average. The comparatively strong isotropizing effect of Mie scattering also explains why the intensity on the tilt axis is never fully depleted.

As shown in Figure 6, the new model significantly improves in matching the flasher data light curves in terms of both timing and total intensity with regards to older models and overall achieves excellent data–simulation agreement. While these light curves only represent some of the full data, the average relative deviation of each model from the data in the plots as shown is 8.5 %, 3.3 %, and 2.4 % for the scattering-function-based model, the birefringence-only model, and after including the ad hoc absorption anisotropy, respectively.

Widespread application in physics analyses requires large-scale simulations and is still in preparation. Nevertheless, first tests employing the ice model in direct-fit reconstructions ([26]) of high-energy events ([3]) indicate that the improved data–simulation agreement seen in flasher data also translates to more accurate descriptions of neutrino events.

Outlook

The model presented here is the first time that the ice microstructure has been included in the modeling of ice optical properties at macroscopic scales. Due to the need to include absorption anisotropy in order to arrive at reasonable grain sizes, for which no first-principles explanation is known, there appear to be remaining additional physical effects not fully accounted for by the first-principles model. At this point it remains unclear whether the anisotropic Mie absorption is real or if it is an artifact from incomplete modeling of birefringence effects. It is currently assumed, for example, that the deduced ice crystal properties follow the same layer undulations as the other ice optical properties and are not simply a function of absolute depth. This assumption may be reevaluated in future works.

Inclusion of ice-intrinsic attenuation in the electromagnetic calculations in Sect. 5 may already change the overall diffusion patterns. In addition, birefringent materials also exhibit di-attenuation where the imaginary index of refraction is polarization- and direction-dependent (Grechushnikov and Konstantinova, 1988). The overall imaginary refractive index of ice is largely unknown (with upper limits derived from IceCube/AMANDA measurements as mentioned earlier; Ackermann et al., 2006) and di-attenuation of ice in the optical has, to our knowledge, not been studied at all. A first step in exploring these options will be to include per crystal (di-)attenuation in the electromagnetic modeling, with the complex refractive indices as free parameters and fitting required values given different assumptions on the crystal orientation fabric. Yet, since the ice-intrinsic absorption accounts for at most 10 % of the overall absorption and the fitted absorption anisotropy is stronger than that, di-attenuation is unlikely to fully explain the observed effect.

The only other known and currently neglected birefringence effect is photoelasticity. Photoelasticity describes the change in refractive index due to applied stresses and is a property of all dielectric media, including ice. Ice is anecdotally known (e.g., [49]) to exhibit strong photoelasticity compared to its intrinsic birefringence strength. Yet, the stress optical parameters have so far only been measured by Ravi-Chandar et al. (1994) for light of an unspecified wavelength at unspecified temperature and only for light propagating along the cc axis. Ravi-Chandar et al. (1994) arrived at a material fringe value of ∼67 kN m−1\sim67\ \mathrm{kN}\,\mathrm{m}^{-1}. Taking this measurement at face value, unrealistically large internal stress of roughly 200 MPa would be required to match the unstressed difference in refractive index.

To investigate the potential relevance of photoelasticity for light diffusion in deep glacial ice, the first step will be to repeat the Ravi-Chandar et al. (1994) measurement and extend it to light propagating orthogonal to the cc axis. If photoelasticity adds a significant contribution, it would allow the presented measurement to also probe the stress state of the sampled ice in addition to the already studied microstructure.

IceCube Upgrade

The IceCube Upgrade [51], planned for deployment in 2025/2026, marks the first extension of the IceCube detector. Over 700 additional modules, including a number of stand-alone calibration devices [47, 76], will be deployed on seven additional strings. Of particular interest for the anisotropy are 11 so-called pencil beam devices. They allow a laser-like beam to be directed in arbitrary directions, enabling sweeps over receiver directions. The birefringence-induced deflection yields a unique signature, where the emission direction of maximum received intensity is offset from the geometric direction to the receiver. Measuring sweeping profiles for several emitter–receiver pairs at different orientations will allow us to disentangle absorption and birefringence contributions to the anisotropy with high precision.

Borehole logging

The described measurement is particularly tailored to the IceCube experiment. Nevertheless, the optical anisotropy effect may still prove to be a useful tool for glaciology. As described by Rongen et al. (2020), most likely fabric-induced azimuthal anisotropy was also observed in the back-scattered intensity recorded by an optical dust logger deployed down the SPC14 drill hole.

To date, the measurement has only been described qualitatively. An accurate simulation of the back-scattering scenario would need to include a good model for the large-angle tail of the Mie scattering function, which is currently poorly constrained from IceCube data. Given a better understanding of back-scattering processes, for example derived using the pencil beam described above, optical logging of drill holes could become a complementary tool for fabric, crystal size, and elongation studies and find wider application in glaciology.

Conclusions

Measurements of ice optical properties in the context of the calibration of the IceCube Neutrino Observatory and its predecessor AMANDA offer unique insights into the properties of glacial ice. In the past, modeling and measurements focused on the impact of airborne impurities as deposited with the original snow accumulation on absorption and scattering and their stratigraphy. This, in particular, yielded the most stringent upper limit as compiled by Warren and Brandt (2008) on the absorption coefficient of pure ice, as measured in the deepest parts of the detector.

Here we have described the observation of an ice optical anisotropy, a direction-dependent intensity modulation aligned with the local ice flow axis. The effect has been identified to largely result from diffusion within the polycrystalline ice microstructure, resulting in a previously unknown optical effect: a slow but continuous deflection towards the normal vector of the girdle plane of the crystal orientation fabric. Combining prior knowledge about the crystal orientation fabric and average grain elongation as obtained from SPC14, the depth-dependent average crystal size has been fitted to IceCube LED calibration data. The resulting depth evolution conforms to the expectation of larger crystals at greater depth and an inverse correlation with impurity concentrations.

The first-principles birefringence explanation was not able to fully describe the experimental data. This has been improved upon by including ad hoc Mie absorption anisotropy, for which no first-principles explanation is known. The origin of this remaining discrepancy will hopefully be resolved using upcoming instrumentation in the IceCube Upgrade, modeling of ice-intrinsic di-attenuation, and future lab measurements regarding the photoelasticity of ice.

Overall the large variety of measurements performed in close vicinity to the Amundsen–Scott South Pole Station (optical data from IceCube and its upcoming detector upgrade, the SPC14 ice core, ground-penetrating radar data from PolarGap, Forsberg et al. (2015); GPS stake field data, Lilien et al. (2018) make the geographic South Pole a unique laboratory for comparative measurements. Yet to date, the overlap in sampled depth between SPC14 and IceCube is unfortunately too small to allow for quantitative comparison. This may be resolved by future drilling projects such as a potential deployment of the Rapid Access Ice Drill (RAID) (Goodge et al., 2021).

Appendix A: Sampling c-axis distributions from the eigenvalues of ice fabric orientation tensors

One can describe the crystal orientation fabric of NN cc axes, measured in an ice sample, by NN unit vectors ni\mathbf{n}_i, with components nixn_{ix}, niyn_{iy}, and nizn_{iz}. Note that ni\mathbf{n}_i is equivalent to −ni-\mathbf{n}_i as the vector can be chosen to point along either direction of the axis. By convention ni\mathbf{n}_i is chosen to point upward. This ensemble of vectors can be represented via the following matrix (Scheidegger, 1965).

a=[∑nix2∑nix⋅niy∑nix⋅niz∑niy⋅nix∑niy2∑niy⋅niz∑niz⋅nix∑niz⋅niy∑niz2](23)a = \begin{bmatrix} \sum n_{ix}^{2} & \sum n_{ix}\cdot n_{iy} & \sum n_{ix}\cdot n_{iz} \\ \sum n_{iy}\cdot n_{ix} & \sum n_{iy}^{2} & \sum n_{iy}\cdot n_{iz} \\ \sum n_{iz}\cdot n_{ix} & \sum n_{iz}\cdot n_{iy} & \sum n_{iz}^{2} \end{bmatrix} \tag*{(23)}

The normalized form A=a/NA = a/N is called the second-order orientation tensor. It was introduced in glaciology through

TypeAbsorptionScatteringModeling
ImpuritiesTotal
SootstrongRayleigh (isotropic)combined absorption and scattering coefficients in 10 m tilt corrected layers (see Sects. 2 and 3)
Mineral duststrongMie (forward)
SaltsweakMie (forward)
AcidsweakMie (forward)
Polycrystalline microstructurenoneAsymmetric diffusionscattering and deflection (see Sects. 5 and 6)

Table A1. Conceptual overview of different constituents considered as part of the ice optical modeling. For details on the behavior of different impurities see [14] and [46]. The polycrystalline microstructure leading to asymmetric diffusion is newly considered in this work.

[43]. AA has three eigenvectors and three corresponding eigenvalues S1S_1, S2S_2, and S3S_3, with S1+S2+S3=1S_1+S_2+S_3=1.

The axes of the coordinate system in which the cc axes are evaluated can be chosen such that the xx axis points along the mean cc-axis direction ∑ini/N\sum_i \mathbf{n}_i/N, that the zz axis points along the pole to the best-fit girdle to the distribution (see [95]), and that the yy axis is orthogonal to the other two. In this case the coordinate axes are the eigenvectors, Sj=∑inij2S_j=\sum_i n_{ij}^{2}, and the eigenvalues follow a strict ordering such that S1≥S2≥S3S_1 \ge S_2 \ge S_3.

A perfectly uniform girdle or single-pole fabric features the following relations between the eigenvalues.

  • uniform: S1≈S2≈S3≈1/3S_1 \approx S_2 \approx S_3 \approx1/3

  • single pole: S1≈1S_1 \approx1; S2≈S3≈0S_2 \approx S_3 \approx0

  • girdle: S1≈S2≈0.5S_1 \approx S_2 \approx0.5; S3≈0S_3 \approx0

[95] realized that many commonly encountered fabric states can be visualized in a 2D plot, as only two of the three eigenvalues are independent. He suggested the representation where the abscissa is given as ln⁡(S2/S3)\ln(S_2/S_3) and the ordinate is given as ln⁡(S1/S2)\ln(S_1/S_2). In this representation uniform cc-axis distributions are found at the origin of the plot. The distance from the origin C=ln⁡(S1/S3)C=\ln(S_1/S_3) is called the strength parameter. Girdle fabrics are found to the lower right, while single-pole fabrics reside to the upper left. The type of fabric can also be quantified by the so-called Woodcock shape parameter K=ln⁡(S1/S2)/ln⁡(S2/S3)K=\ln(S_1/S_2)/\ln(S_2/S_3). Large KK values denote a single-pole fabric. KK values smaller than 1 denote a girdle fabric. Figure A1 presents the fabric versus depth evolution as measured at the geographic South Pole in this representation.

Neither the orientation tensor nor its eigenvalues retain the full information on the ensemble of underlying cc axes. Thus, an assumption on the functional form of the fabric has to be made when trying to sample a distribution. Here we focus on describing random, girdle, and single-pole distributions, as well as combinations of these, as those are the types most commonly encountered in ice fabric measurements, with SPC14 in particular featuring a very strong girdle fabric.

The book Statistical analysis of spherical data by [37] gives a good overview of commonly used probability density functions (PDFs) for directional data. Of the presented PDFs the Watson (1965) distribution seems most applicable for our case due to the following.

  1. It can represent both unimodal and single-pole as well as rotational symmetric girdle data.

  2. An (approximate) parameter estimation exists based on eigenvalues alone.

In its standardized form the PDF, evaluated on a spherical coordinate system with the polar angle θ\theta and the azimuth angle ϕ\phi, has only one free parameter κ\kappa and is given as

f(θ,ϕ)=Cwexp⁡(κ⋅cos⁡2θ)sin⁡θ,(24)f(\theta,\phi)=C_{\mathrm{w}}\exp(\kappa\cdot\cos^{2}\theta)\sin\theta, \tag*{(24)}

with the normalization constant

Cw=14π∫01exp⁡(κ⋅u2)du.(25)C_{\mathrm{w}}=\frac{1}{4\pi\int_{0}^{1}\exp(\kappa\cdot u^{2})\mathrm{d}u}. \tag*{(25)}

In the following, the Python package available at https://github.com/duncandc/watson_distribution (last access: 20 December 2023) is used to sample from the Watson distribution. Alternatively, the sampling approach described in [37] may be used.

At κ=0\kappa=0 the direction distribution is perfectly uniform. For positive κ\kappa the distribution is bimodal in vector space, which is equivalent to a single-pole distribution in axis space, and has the highest probability at the poles. For negative κ\kappa values the distribution is girdle with the directions equally distributed around the Equator.

[21] showed that for a purely single-pole distribution the κ\kappa parameter can be estimated from the eigenvalues as

κ={3.75⋅(3⋅S1−1),13≤S1≤0.34−5.95+14.9S1+1.481−S1−11.05S12,0.34<S1≤0.64−7.96+21.5⋅S1+1.481−S1−13.25⋅S12,S1>0.64,(26)\kappa= \begin{cases} 3.75\cdot(3\cdot S_{1}-1), & \frac{1}{3}\leq S_{1}\leq0.34\\ -5.95+14.9S_{1}+\frac{1.48}{1-S_{1}}-\frac{11.05}{S_{1}^{2}}, & 0.34<S_{1}\leq0.64\\ -7.96+21.5\cdot S_{1}+\frac{1.48}{1-S_{1}}-13.25\cdot S_{1}^{2}, & S_{1}>0.64, \end{cases} \tag*{(26)}

while for a purely girdle fabric the κ\kappa parameter can be estimated as

κ={12⋅S3,0≤S3≤0.060.961−7.08⋅S3+0.466S3,0.06<S3≤0.323.75⋅(1−3⋅S3),0.32<S3≤13.(27)\kappa= \begin{cases} \frac{1}{2\cdot S_{3}}, & 0\leq S_{3}\leq0.06\\ 0.961-7.08\cdot S_{3}+\frac{0.466}{S_{3}}, & 0.06<S_{3}\leq0.32\\ 3.75\cdot(1-3\cdot S_{3}), & 0.32<S_{3}\leq\frac{1}{3}. \end{cases} \tag*{(27)}

For ice fabrics, the plane of girdle c axes shall intersect the poles, where the c axes of a single-pole distribution are also found. As such the directions sampled from girdle Watson distributions are rotated by 90∘90^{\circ}. Due to the underlying rotational symmetry the eigenvalues of the resulting Watson distributions follow a strict relation.

S1=S2&S3=1−2⋅S1for a girdle WatsonS2=S3&S1=1−2⋅S2for a unimodal Watson(28)\begin{aligned} S_{1}&=S_{2}\quad\&\quad S_{3}=1-2\cdot S_{1}\quad\text{for a girdle Watson}\\ S_{2}&=S_{3}\quad\&\quad S_{1}=1-2\cdot S_{2}\quad\text{for a unimodal Watson} \tag*{(28)} \end{aligned}

Obviously no single Watson distribution can describe an arbitrary set of eigenvalues with S1≠S2≠S3S_{1}\ne S_{2}\ne S_{3}. This is achieved by combining directions sampled from a girdle and a unimodal Watson distribution.

Given a sample of c axes from a girdle Watson distribution with eigenvalues SigS_{ig} and a sample of c axes from a unimodal Watson distribution with eigenvalues SiuS_{iu}, as well as a relative

fractional contribution of the girdle sample fgf_g to the total sample, the eigenvalues of the combined sample SiS_i are given by Si=fg⋅Sig+(1−fg)⋅SiuS_i=f_g\cdot S_{ig}+(1-f_g)\cdot S_{iu} The combination of fgf_g, S1gS_{1g}, and S2uS_{2u} which yields the desired eigenvalues S1S_1, S2S_2, and S3S_3 is found by solving the equation system SiS_i, which has been simplified using the relations in (28):

S1=fg⋅S1g+(1−fg)⋅(1−2⋅S2u),S2=fg⋅S1g+(1−fg)⋅S2u,S3=fg⋅(1−2⋅S1g)+(1−fg)⋅S2u.(29)\begin{aligned} S_{1}&=f_g\cdot S_{1g}+(1-f_g)\cdot(1-2\cdot S_{2u}),\\ S_{2}&=f_g\cdot S_{1g}+(1-f_g)\cdot S_{2u},\\ S_{3}&=f_g\cdot(1-2\cdot S_{1g})+(1-f_g)\cdot S_{2u}. \tag*{(29)} \end{aligned}

The third equation is not independent since S1+S2+S3=1S_1+S_2+S_3=1. Thus further information is needed to be able to constrain the variables. To fulfill the assumption that SuS_u is unimodal and single-pole and SgS_g is girdle one can further constrain 1/3<S1g<0.51/3<S_{1g}<0.5 and 0<S2u<1/30<S_{2u}<1/3. For cases in which the system is still underconstrained one can, for example, further demand that S2g=S2u=S1gS_{2g}=S_{2u}=S_{1g} (equivalent to S1u+S3g=1S_{1u}+S_{3g}=1) so that both distributions have an equal spread around the girdle plane. The solution is then given as

fg=0.5⋅(ϵ−4S1−2S2+3),S1g=2S1−ϵ+4S1+2S2+1,S2u=ϵ−2S2−12⋅(ϵ−4S1−2S2−1).(30)\begin{aligned} f_g&=0.5\cdot(\epsilon-4S_{1}-2S_{2}+3),\\ S_{1g}&=\frac{2S_{1}}{-\epsilon+4S_{1}+2S_{2}+1},\\ S_{2u}&=\frac{\epsilon-2S_{2}-1}{2\cdot(\epsilon-4S_{1}-2S_{2}-1)}. \tag*{(30)} \end{aligned}

with ϵ=16⋅S12+16S1⋅(S2−1)+(2S2+1)2\epsilon= \sqrt{16 \cdot S_1^2 + 16S_1 \cdot(S_2 - 1) + (2S_2 + 1)^2}. From these one can derive the Watson parameters κ\kappa using the approximations as given in Eqs. (A4) and (A5).

To verify and visualize the success of the presented sampling approach, cc-axis distributions according to a number of combinations of ln⁡(S1/S2)\ln(S_1/S_2) and ln⁡(S2/S3)\ln(S_2/S_3) have been generated as shown in Figure A2. The sampled cc-axis distributions yield eigenvalues which are accurate to within the approximation of the parameter estimation for the Watson distributions and sufficient for most applications.

Example c-axis distributions in Schmidt equal-area projection

Figure A2. Example c-axis distributions in Schmidt equal-area projection generated using the described method.

Note that by design the cc-axis distributions for intermediate fabric states do not contain a single elliptical distribution but a rotationally symmetric girdle and a circular single pole. This seems suitable for our application to ice fabrics. In very deep glacial ice where the fabric slowly evolves from girdle to single pole, experimental distributions such as published by [89] indeed show the described superposition and not an elliptical distribution usually sketched for these eigenvalues.

Appendix B: Sampling surface orientations from an ellipsoid

As the average grain shape deviates from a sphere, the encountered distribution of face orientations depends on the photon direction. Assuming that the face orientation of a solid, tessellated into elongated polyhedra, to be described by the surface orientation density of an ellipsoid describing the average grain shape, one can sample the distribution as follows.

The surface of an ellipsoid is defined by the equation

f(x′,y′,z′)=x′2a2+y′2b2+z′2c2=1,(31)f(x',y',z') = \frac{x'^2}{a^2} + \frac{y'^2}{b^2} + \frac{z'^2}{c^2} = 1, \tag*{(31)}

where aa, bb, and cc are the dimensions of the major and minor axes. The normal vector on any point of the surface is given by the gradient

∇f=[2⋅x′a2,2⋅y′b2,2⋅z′c2].(32)\nabla f = \left[2 \cdot\frac{x'}{a^2}, 2 \cdot\frac{y'}{b^2}, 2 \cdot\frac{z'}{c^2}\right]. \tag*{(32)}

For a given set of azimuth and zenith angles, the coordinates on a unit sphere (x,y,z)(x,y,z) and on the ellipsoid (x′,y′,z′)(x',y',z') are given as

x=sin⁡θ⋅cos⁡ϕ and x′=a⋅x,(33)x = \sin\theta\cdot\cos\phi\text{ and } x' = a \cdot x, \tag*{(33)}
y=sin⁡θ⋅sin⁡ϕ and y′=b⋅y,(34)y = \sin\theta\cdot\sin\phi\text{ and } y' = b \cdot y, \tag*{(34)}
z=cos⁡θ and z′=c⋅z.(35)z = \cos\theta\text{ and } z' = c \cdot z. \tag*{(35)}

Substituting the ellipsoid surface position into (32) the surface normal at this position is then

n=[2a⋅sin⁡θ⋅cos⁡ϕ,2b⋅sin⁡θ⋅sin⁡ϕ,2c⋅cos⁡θ].(36)\mathbf{n} = \left[\frac{2}{a} \cdot\sin\theta\cdot\cos\phi, \frac{2}{b} \cdot\sin\theta\cdot\sin\phi, \frac{2}{c} \cdot\cos\theta\right]. \tag*{(36)}

One can now sample these gradients with angles chosen to be uniform on a sphere. As the surface density per solid angle of an ellipsoid is different from a sphere, the relative surface density,

μ(x,y,z)=∥dS′∥∥dS∥=(ac⋅y)2+(ab⋅z)2+(bc⋅x)2,(37)\begin{aligned} \mu(x,y,z) = \frac{\lVert\mathrm{d}S'\rVert}{\lVert\mathrm{d}S\rVert} \\ &= \sqrt{(ac \cdot y)^2 + (ab \cdot z)^2 + (bc \cdot x)^2}, \tag*{(37)} \end{aligned}

has to be applied as a weighting factor, where the maximum weighting factor is given as

μmax⁡=max⁡(ac,ab,bc).(38)\mu_{\max} = \max(ac,ab,bc). \tag*{(38)}

Instead of weighting, one can also employ a rejection sampling with an acceptance probability of μ/μmax⁡\mu/\mu_{\max}.

In addition to the distribution of face orientations, the distribution of face orientations actually encountered by a photon can be obtained by weighting the distribution of face orientations with the scalar product of the photon’s propagation vector and each face normal vector. The probability of encountering a given plane is therefore simply the projected area relative to the incident light.

Figure B1 shows the cos⁡(θ)\cos(\theta) distribution of (encountered) face normal vectors for a spheroid with elongation 2, which has the major axis aligned with the zz axis. The distribution is compared to a crystal-like Voronoi tessellation generated with Neper and assuming the same mean elongation. Lines have been traced through the tessellation, identifying grain boundary encounters and computing their incidence angles. The distributions are found to be indistinguishable, confirming that the ensemble of polyhedra faces follows the average ellipsoid.

Ellipsoid surface sampling for an ellipsoid with unity minor axes and a major axis of two along the z axis

Figure B1. Ellipsoid surface sampling for an ellipsoid with unity minor axes and a major axes of two along the zz axis. Green: analytic cos⁡(θ)\cos(\theta) distribution of face normal vectors. Red: analytic cos⁡(θ)\cos(\theta) distribution weighted by the encounter probability, given by the scalar product with a photon propagating along zz. Blue: encounter probability as found in a Neper crystal tessellation simulation when tracing photons along vertical lines.

Code and data availability

The photon propagator software (PPC), compatible ice model configurations including the model derived in this work, and the electromagnetics code used to generate the diffusion patterns are available from https://doi.org/10.5281/zenodo.10410726 [28]. IceCube raw data, including the LED calibration data, are generally not publicly available. For specific inquiries please contact [email protected].

Author contributions

The IceCube collaboration designed, constructed, and now operates the IceCube Neutrino Observatory. Data processing and calibration, Monte Carlo simulations of the detector and of theoretical models, and data analyses were performed by a large number of collaboration members, who also discussed and approved the scientific results presented here. The IceCube collaboration acknowledges the substantial contributions to this paper from MR and DCh. The paper was reviewed by the entire collaboration before publication, and all authors approved the final version.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

The authors gratefully acknowledge support from the following agencies and institutions: USA – U.S. National Science Foundation, Office of Polar Programs; U.S. National Science Foundation, Physics Division; U.S. National Science Foundation, EPSCoR; Wisconsin Alumni Research Foundation; Center for High Throughput Computing (CHTC) at the University of Wisconsin–Madison; Open Science Grid (OSG); Extreme Science and Engineering Discovery Environment (XSEDE); Frontera computing project at the Texas Advanced Computing Center; U.S. Department of Energy, National Energy Research Scientific Computing Center; particle astrophysics research computing center at the University of Maryland; Institute for Cyber-Enabled Research at Michigan State University; and astroparticle physics computational facility at Marquette University. Belgium – Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programs, and Belgian Federal Science Policy Office (Belspo). Germany – Bundesministerium für Bildung und Forschung (BMBF), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astro-particle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and high-performance computing cluster of the RWTH Aachen. Sweden – Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation. Australia – Australian Research Council. Canada – Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Compute Canada. Denmark – Villum Fonden, Carlsberg Foundation, and European Commission. New Zealand – Marsden Fund. Japan – Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University. Korea – National Research Foundation of Korea (NRF). Switzerland – Swiss National Science Foundation (SNSF). United Kingdom – Department of Physics, University of Oxford.

Review statement

This paper was edited by Carlos Martin and reviewed by David Lilien and one anonymous referee.

References

References

  1. [1]Aab, A., Abreu, P., Aglietta, M., et al.: Features of the Energy Spectrum of Cosmic Rays above 2.5 × 10¹⁸ eV Using the Pierre Auger Observatory, Phys. Rev. Lett., 125, 121106, https://doi.org/10.1103/PhysRevLett.125.121106, 2020.
  2. [2]Aartsen, M. G., Abbasi, R., Abdou, Y., et al.: South Pole glacial climate reconstruction from multi-borehole laser particulate stratigraphy, J. Glaciol., 59, 1117–1128, https://doi.org/10.3189/2013JoG13J068, 2013a.
  3. [3]Aartsen, M. G., Abbasi, R., Abdou, Y., et al.: Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector, Science, 342, https://doi.org/10.1126/science.1242856, 2013b.
  4. [5]Aartsen, M. G., Abbasi, R., Ackermann, M., et al.: Energy Reconstruction Methods in the IceCube Neutrino Telescope, J. Instr., P03009, https://doi.org/10.1088/1748-0221/9/03/P03009, 2013d.
  5. [6]Aartsen, M. G., Ackermann, M., Adams, J., et al.: Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data, Phys. Rev. Lett., 113, 101101, https://doi.org/10.1103/PhysRevLett.113.101101, 2014.
  6. [7]Aartsen, M. G., Abraham, K., Ackermann, et al.: Search for correlations between the arrival directions of IceCube neutrino events and ultrahigh-energy cosmic rays detected by the Pierre Auger Observatory and the Telescope Array, J. Cosmol. Astroparticle Physics, https://doi.org/10.1088/1475-7516/2016/01/037, 2015.
  7. [8]Aartsen, M. G., Ackermann, M., Adams, J., et al.: The IceCube Neutrino Observatory: Instrumentation and Online Systems, J. Instr., 12, P03012, https://doi.org/10.1088/1748-0221/12/03/P03012, 2017.
  8. [9]Aartsen, M., Ackermann, M., Adams, J., et al.: Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert, Science, 361, 147–151, https://doi.org/10.1126/science.aat2890, 2018a.
  9. [10]Aartsen, M. G., Ackermann, M., Adams, J., et al.: Measurement of Atmospheric Neutrino Oscillations at 6–56 GeV with IceCube DeepCore, Phys. Rev. Lett., 120, 071801, https://doi.org/10.1103/PhysRevLett.120.071801, 2018b.
  10. [11]Abbasi, R., Ackermann, M., Adams, J., et al.: The IceCube data acquisition system: Signal capture, digitization, and timestamping, Nucl. Instrum. Meth. A, 601, 294–316, https://doi.org/10.1016/j.nima.2009.01.001, 2009.
  11. [12]Abbasi, R., Abdou, Y., Abu-Zayyad, T., et al.: Calibration and characterization of the IceCube photomultiplier tube, Nucl. Instrum. Meth. A, 618, 139–152, https://doi.org/10.1016/j.nima.2010.03.102, 2010.
  12. [14]Ackermann, M., Ahrens, J., Bai, X., et al.: Optical properties of deep glacial ice at the South Pole, J. Geophys. Res.-Atmos., 111, https://doi.org/10.1029/2005JD006687, 2006.
  13. [15]Albert, A., André, M., Anghinolfi, M., et al.: Combined search for neutrinos from dark matter self-annihilation in the Galactic Center with ANTARES and IceCube, Phys. Rev. D, 102, 082002, https://doi.org/10.1103/physrevd.102.082002, 2020.
  14. [16]Alley, R. B.: Fabrics in Polar Ice Sheets: Development and Prediction, Science, 240, 493–495, https://doi.org/10.1126/science.240.4851.493, 1988.
  15. [17]Alley, R. B., Fitzpatrick, J., Fegyveresi, J., and Voigt, D.: Physical properties of the South Pole Ice Core, SPC14, IceCube Polar Science Workshop, https://events.icecube.wisc.edu/event/128/contributions/7313/ (last access: 20 December 2023), 2021.
  16. [18]Allgaier, M., Cooper, M. G., Carlson, A. E., Cooley, S. W., Ryan, J. C., and Smith, B. J.: Direct measurement of optical properties of glacier ice using a photon-counting diffuse LiDAR, J. Glaciol., 68, 1–11, https://doi.org/10.1017/jog.2022.34, 2022.
  17. [19]Bay, R. C., Rohde, R. A., Price, P. B., and Bramall, N. E.: South Pole paleowind from automated synthesis of ice core records, J. Geophys. Res.-Atmos., 115, https://doi.org/10.1029/2009JD013741, 2010.
  18. [20]Benson, T., Cherwinka, J., Duvernois, M., Elcheikh, A., Feyzi, F., Greenler, L., Haugen, J., Karle, A., Mulligan, M., and Paulos, R.: IceCube Enhanced Hot Water Drill functional description, Ann. Glaciol., 55, 105–114, https://doi.org/10.3189/2014ao g68a032, 2014.
  19. [21]Best, D. J. and Fisher, N. I.: Goodness-of-fit and discordancy tests for samples from the Watson distribution on the sphere, Aust. J. Stat., 28, 13–31, https://doi.org/10.1111/j.1467-842X.1986.tb00580.x, 1986.
  20. [22]Casey, K. A., Fudge, T. J., Neumann, T. A., Steig, E. J., Cavitte, M. G. P., and Blankenship, D. D.: The 1500 m South Pole ice core: recovering a 40 ka environmental record, Ann. Glaciol., 55, 137–146, https://doi.org/10.3189/2014aog68a016, 2014.
  21. [23]Cherenkov, P. A.: Visible Radiation Produced by Electrons Moving in a Medium with Velocities Exceeding that of Light, Phys. Rev., 52, 378–379, https://doi.org/10.1103/physrev.52.378, 1937.
  22. [24]Chirkin, D.: Photon tracking with GPUs in IceCube, Nucl. Instrum. Meth. A, 725, 141–143, https://doi.org/10.1016/j.nima.2012.11.170, 2013a.
  23. [25]Chirkin, D.: Likelihood description for comparing data with simulation of limited statistics, arxiv, http://arxiv.org/abs/1304.0735v3 (last access: 20 December 2023), 2013b.
  24. [26]Chirkin, D.: Event reconstruction in IceCube based on direct event re-simulation, ICRC, https://ui.adsabs.harvard.edu/abs/2013ICRC...33.3402C (last access: 20 December 2023), 2013c.
  25. [27]Chirkin, D.: Evidence of optical anisotropy of the South Pole ice, in: Proceedings, 33rd International Cosmic Ray Conference (ICRC2013): Rio de Janeiro, Brazil, 2–9 July 2013, p. 0580, http://inspirehep.net/record/1412998/files/icrc2013-0580.pdf (last access: 20 December 2023), 2013d.
  26. [28]Chirkin, D. A.: icecube/ppc: Version 1.0 (V1.0), Zenodo [code], https://doi.org/10.5281/zenodo.10410726, 2023.
  27. [30]Abbasi, R., Ackermann, M., Adams, J., et al.: A calibration study of local ice and optical sensor properties in IceCube, PoS, ICRC2021, 1023, https://doi.org/10.22323/1.395.1023, 2021.
  28. [31]Cooper, M. A., Jordan, T. M., Siegert, M. J., and Bamber, J. L.: Surface Expression of Basal and Englacial Features, Properties, and Processes of the Greenland Ice Sheet, Geophys. Res. Lett., 46, 783–793, https://doi.org/10.1029/2018gl080620, 2019.
  29. [32]Cuffey, K. M.: The physics of glaciers, Butterworth-Heinemann/Elsevier, Burlington, MA, ISBN 9780123694614, 2010.
  30. [33]Dow, C. F., Karlsson, N. B., and Werder, M. A.: Limited Impact of Subglacial Supercooling Freeze-on for Greenland Ice Sheet Stratigraphy, Geophys. Res. Lett., 45, 1481–1489, https://doi.org/10.1002/2017GL076251, 2018.
  31. [34]Durand, G., Weiss, J., Lipenkov, V., Barnola, J. M., Krinner, G., Parrenin, F., Delmonte, B., Ritz, C., Duval, P., Röthlisberger, R., and Bigler, M.: Effect of impurities on grain growth in cold ice sheets, J. Geophys. Res., 111, https://doi.org/10.1029/2005jf000320, 2006.
  32. [35]Faria, S. H., Weikusat, I., and Azuma, N.: The microstructure of polar ice. Part I: Highlights from ice core research, J. Struct. Geol., 61, 2–20, https://doi.org/10.1016/j.jsg.2013.09.010, 2014a.
  33. [37]Fisher, N. I., Lewis, T., and Embleton, B. J. J.: Statistical Analysis of Spherical Data, Cambridge University Press, https://doi.org/10.1017/cbo9780511623059, 1987.
  34. [38]Fitzpatrick, J. J., Voigt, D. E., Fegyveresi, J. M., Stevens, N. T., Spencer, M. K., Cole-Dai, J., Alley, R. B., Jardine, G. E., Cravens, E. D., Wilen, L. A., Fudge, T., and McConnell, J. R.: Physical properties of the WAIS Divide ice core, J. Glaciol., 60, 1181–1198, https://doi.org/10.3189/2014jog14j100, 2014.
  35. [40]Fujita, S., Maeno, H., Uratsuka, S., Furukawa, T., Mae, S., Fujii, Y., and Watanabe, O.: Nature of radio echo layering in the Antarctic Ice Sheet detected by a two-frequency experiment, J. Geophys. Res.-Sol. Ea., 104, 13013–13024, https://doi.org/10.1029/1999JB900034, 1999.
  36. [41]Fujita, S., Maeno, H., and Matsuoka, K.: Radio-wave depolarization and scattering within ice sheets: a matrix-based model to link radar and ice-core measurements and its application, J. Glaciol., 52, 407–424, https://doi.org/10.3189/172756506781828548, 2006.
  37. [42]Gebhart, J.: Response of Single-Particle Optical Counters to Particles of irregular shape, Particle & Particle Systems Characterization, 8, 40–47, https://doi.org/10.1002/ppsc.19910080109, 1991.
  38. [43]Gödert, G. and Hutter, K.: Induced anisotropy in large ice shields: theory and its homogenization, Continuum Mechanics and Thermodynamics, 10, 293–318, https://doi.org/10.1007/s001610050095, 1998.
  39. [46]He, Y. D. and Price, P. B.: Remote sensing of dust in deep ice at the South Pole, J. Geophys. Res.-Atmos., 103, 17041–17056, https://doi.org/10.1029/98jd01643, 1998.
  40. [47]Henningsen, F., Böhmer, M., Gärtner, A., Geilen, L., Gernhäuser, R., Heggen, H., Holzapfel, K., Fruck, C., Papp, L., Rea, I. C., Resconi, E., Schmuckermaier, F., Spannfellner, C., and Traxler, M.: A self-monitoring precision calibration flight source for large-volume neutrino telescopes, J. Instr., 15, P07031, https://doi.org/10.1088/1748-0221/15/07/P07031, 2020.
  41. [48]Henyey, L. C. and Greenstein, J. L.: Diffuse radiation in the Galaxy, Astrophys. J., 93, 70, https://doi.org/10.1086/144246, 1941.
  42. [49]Hobbs, P.: Ice physics, Oxford University Press, New York, ISBN 9780199587711, 2010.
  43. [51]Ishihara, A.: The IceCube Upgrade — Design and Science Goals, PoS, ICRC2019, 1031, https://doi.org/10.22323/1.358.1031, 2021.
  44. [52]Johnson, J. A., Shturmakov, A. J., Kuhl, T. W., Mortensen, N. B., and Gibson, C. J.: Next generation of an intermediate depth drill, Ann. Glaciol., 55, 27–33, https://doi.org/10.3189/2014aog68a011, 2014.
  45. [53]Jordan, T. M., Schroeder, D. M., Castelletti, D., Li, J., and Dall, J.: A polarimetric coherence method to determine ice crystal orientation fabric from radar sounding: application to the NEEM Ice Core Region, IEEE T. Geosci. Remote, 57, 8641–8657, 2019.DOI
  46. [54]Kluskiewicz, D., Waddington, E. D., Anandakrishnan, S., Voigt, D. E., Matsuoka, K., and McCarthy, M. P.: Sonic methods for measuring crystal orientation fabric in ice, and results from the West Antarctic ice sheet (WAIS) Divide, J. Glaciol., 63, 603–617, https://doi.org/10.1017/jog.2017.20, 2017.
  47. [55]Landau, L. and Lifshitz, E.: Electrodynamics of Continuous Media, Addison-Wesley, ISBN 0080091059, 1960.
  48. [56]Langway, C. C.: Ice fabrics and the universal stage, vol. 62, Department of Defense, Department of the Army, Corps of Engineers, Snow Ice, https://hdl.handle.net/11681/6005, 1958.
  49. [60]Lilien, D. A., Fudge, T. J., Koutnik, M. R., Conway, H., Osterberg, E. C., Ferris, D. G., Waddington, E. D., and Stevens, C. M.: Holocene Ice-Flow Speedup in the Vicinity of the South Pole, Geophys. Res. Lett., 45, 6557–6565, https://doi.org/10.1029/2018GL078253, 2018.
  50. [61]Lipenkov, V., Barkov, N., Duval, P., and Pimienta, P.: Crystalline Texture of the 2083 m Ice Core at Vostok Station, Antarctica, J. Glaciol., 35, 392–398, https://doi.org/10.1017/S0022143000009321, 1989.
  51. [62]Liu, P.: A new phase function approximating to Mie scattering for radiative transport equations, Phys. Med. Biol., 39, 1025, https://doi.org/10.1088/0031-9155/39/6/008, 1994.
  52. [63]MacGregor, J. A., Fahnestock, M. A., Catania, G. A., Paden, J. D., Prasad Gogineni, S., Young, S. K., Rybarski, S. C., Mabrey, A. N., Wagman, B. M., and Morlighem, M.: Radiostratigraphy and age structure of the Greenland Ice Sheet, J. Geophys. Res.-Earth Surf., 120, 212–241, https://doi.org/10.1002/2014JF003215, 2015.
  53. [64]Margolis, S. H., Schramm, D. N., and Silberberg, R.: Ultrahigh-energy neutrino astronomy, Astrophys. J., 221, 990, https://doi.org/10.1086/156104, 1978.
  54. [65]Markov, M. A.: On high energy neutrino physics, in: Proceedings, 10th International Conference on High-Energy Physics (ICHEP 60): Rochester, NY, USA, 25 August–1 September 1960, 578–581, 1960.
  55. [66]Matsuoka, K., Furukawa, T., Fujita, S., Maeno, H., Uratsuka, S., Naruse, R., and Watanabe, O.: Crystal orientation fabrics within the Antarctic ice sheet revealed by a multipolarization plane and dual-frequency radar survey, J. Geophys. Res.-Sol. Ea., 108, https://doi.org/10.1029/2003JB002425, 2003.
  56. [67]McConnel, J. C.: On the plasticity of an ice crystal, P. Roy. Soc. Lond., 49, 323–343, https://doi.org/10.1098/rspl.1890.0099, 1891.
  57. [68]Paxman, G. J. G., Jamieson, S. S. R., Ferraccioli, F., Jordan, T. A., Bentley, M. J., Ross, N., Forsberg, R., Matsuoka, K., Steinhage, D., Eagles, G., and Casal, T. G.: Subglacial Geology and Geomorphology of the Pensacola-Pole Basin, East Antarctica, Geochem. Geophys. Geosyst., 20, 2786–2807, https://doi.org/10.1029/2018GC008126, 2019.
  58. [69]Petrenko, V. F. and Whitworth, R. W.: Physics of Ice, Oxford University Press, https://doi.org/10.1093/acprof:oso/9780198518945.001.0001, 2002.
  59. [70]Potenza, M. A. C., Albani, S., Delmonte, B., Villa, S., Sanvito, T., Paroli, B., Pullia, A., Baccolo, G., Mahowald, N., and Maggi, V.: Shape and size constraints on dust optical properties from the Dome C ice core, Antarctica, Sci. Rep., 6, 28162, https://doi.org/10.1038/srep28162, 2016.
  60. [71]Price, P. B. and Bergström, L.: Optical properties of deep ice at the South Pole: scattering, Appl. Optics, 36, 4181–4194, https://doi.org/10.1364/AO.36.004181, 1997. Price, P. B., Woschnagg, K., and Chirkin, D.: Age vs depth of glacial ice at South Pole, Geophys. Res. Lett., 27, 2129–2132, https://doi.org/10.1029/2000GL011351, 2000.
  61. [72]Quey, R., Dawson, P., and Barbe, F.: Large-scale 3D random polycrystals for the finite element method: Generation, meshing and remeshing, Comput. Method. Appl. M., 200, https://doi.org/10.1016/j.cma.2011.01.002, 2011.
  62. [73]Raman, C. V. and Viswanathan, K. S.: The theory of the propagation of light in polycrystalline media, P. Natl. Acad. Sci. India A, 41, 37–44, https://doi.org/10.1007/BF03047170, 1955.
  63. [75]Rongen, M.: Calibration of the IceCube Neutrino Observatory, PhD thesis, RWTH Aachen University, https://doi.org/10.18154/RWTH-2019-09941, 2019.
  64. [76]Rongen, M. and Chirkin, D.: Advances in IceCube ice modelling & what to expect from the Upgrade, J. Instr., 16, C09014, https://doi.org/10.1088/1748-0221/16/09/C09014, 2021.
  65. [78]Scargle, J. D.: Studies in Astronomical Time Series Analysis. V. Bayesian Blocks, a New Method to Analyze Structure in Photon Counting Data, Astrophys. J., 504, 405–418, https://doi.org/10.1086/306064, 1998.
  66. [80]Simonsen, M. F., Cremonesi, L., Baccolo, G., Bosch, S., Delmonte, B., Erhardt, T., Kjær, H. A., Potenza, M., Svensson, A., and Vallelonga, P.: Particle shape accounts for instrumental discrepancy in ice core dust size distributions, Clim. Past, 14, 601–608, https://doi.org/10.5194/cp-14-601-2018, 2018.
  67. [81]Stokstad, R.: Design and performance of the IceCube electronics, in: 11th Workshop on Electronics for LHC and Future Experiments (LECC 2005), p. 4, https://cds.cern.ch/record/920022/files/p20.pdf (last access: 20 December 2023), 2005.
  68. [82]Stoll, N., Eichler, J., Hörhold, M., Erhardt, T., Jensen, C., and Weikusat, I.: Microstructure, micro-inclusions, and mineralogy along the EGRIP ice core – Part 1: Localisation of inclusions and deformation patterns, The Cryosphere, 15, 5717–5737, https://doi.org/10.5194/tc-15-5717-2021, 2021.
  69. [84]Uchida, T., Shimada, W., Hondoh, T., Mae, S., and Barkov, N. I.: Refractive-index measurements of natural air–hydrate crystals in an Antarctic ice sheet, Appl. Optics, 34, 5746–5749, https://doi.org/10.1364/AO.34.005746, 1995.
  70. [85]Uchida, T., Miyamoto, A., Shin’yama, A., and Hondoh, T.: Crystal growth of air hydrates over 720 ka in Dome Fuji (Antarctica) ice cores: microscopic observations of morphological changes below 2000 m depth, J. Glaciol., 57, 1017–1026, https://doi.org/10.3189/002214311798843296, 2011.
  71. [86]Voigt, D. E.: c-Axis Fabric of the South Pole Ice Core, SPCI14, https://doi.org/10.15784/601057, 2017.
  72. [87]Warren, S. G. and Brandt, R. E.: Optical constants of ice from the ultraviolet to the microwave: A revised compilation, J. Geophys. Res., 113, https://doi.org/10.1029/2007jd009744, 2008.
  73. [89]Weikusat, I., Jansen, D., et al.: Physical analysis of an Antarctic ice core – towards an integration of micro- and macrodynamics of polar ice, Philos. T. Roy. Soc. A, 375, 20150347, https://doi.org/10.1098/rsta.2015.0347, 2017.
  74. [90]Westhoff, J., Stoll, N., Franke, S., Weikusat, I., Bons, P., Kerch, J., Jansen, D., Kipfstuhl, S., and Dahl-Jensen, D.: A stratigraphy-based method for reconstructing ice core orientation, Ann. Glaciol., 62, 191–202, https://doi.org/10.1017/aog.2020.76, 2021.
  75. [91]Wilen, L., Diprinzio, C., Alley, R., and Azuma, N.: Development, principles, and applications of automated ice fabric analyzers, Microscopy Research and Technique, 62, 2–18, 2003.
  76. [92]Wilks, S. S.: The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses, Ann. Math. Statist., 9, 60–62, https://doi.org/10.1214/aoms/1177732360, 1938.
  77. [93]Wilson, C. J. L., Russell-Head, D. S., and Sim, H. M.: The application of an automated fabric analyzer system to the textural evolution of folded ice layers in shear zones, Ann. Glaciol., 37, 7–17, https://doi.org/10.3189/1727564037815401, 2003.
  78. [94]Winski, D. A., Fudge, T. J., Ferris, D. G., Osterberg, E. C., Fegyveresi, J. M., Cole-Dai, J., Thundercloud, Z., Cox, T. S., Kreutz, K. J., Ortman, N., Buizert, C., Epifanio, J., Brook, E. J., Beaudette, R., Seveninghaus, J., Sowers, T., Steig, E. J., Kahle, E. C., Jones, T. R., Morris, V., Aydin, M., Nicewonger, M. R., Casey, K. A., Alley, R. B., Waddington, E. D., Iverson, N. A., Dunbar, N. W., Bay, R. C., Souney, J. M., Sigl, M., and McConnell, J. R.: The SPI9 chronology for the South Pole Ice Core – Part 1: volcanic matching and annual layer counting, Clim. Past, 15, 1793–1808, https://doi.org/10.5194/cp-15-1793-2019, 2019.
  79. [95]Woodcock, N. H.: Specification of fabric shapes using an eigenvalue method, Geol. Soc. Am. B., 88, 1231, https://doi.org/10.1130/0016-7606(1977)88<1231:sofs>2.0.co;2, 1977.
  80. [96]Young, T. J., Martín, C., Christoffersen, P., Schroeder, D. M., Tulaczyk, S. M., and Dawson, E. J.: Rapid and accurate polarimetric radar measurements of ice crystal fabric orientation at the Western Antarctic Ice Sheet (WAIS) Divide ice core site, The Cryosphere, 15, 4117–4133, https://doi.org/10.5194/tc-15-4117-2021, 2021.
  81. [97]Zhang, Z. and Caulfield, H.: Reflection and refraction by interfaces of uniaxial crystals, Opt. Laser Technol., 28, 549–553, https://doi.org/10.1016/S0030-3992(96)00022-9, 1996.

Paper details

Publication

vlnc:W_pj8e4nay
v1 · latestHistory (1)

published · Fresh — published within the last 30 days

views · 0

downloads · 0

A view counts once per browser per day, a download each time the PDF is requested.

CC BY 4.0

Original source

Taxonomy

  • Physics and Astronomy

Authors

AI contributions

Not recorded for this paper · View statement →

Looking for related papers…
Contents