Improved modeling of in-ice particle showers for IceCube event reconstruction

: The IceCube Neutrino Observatory relies on an array of photomultiplier tubes to detect Cherenkov light produced by charged particles in the South Pole ice. IceCube data analyses depend on an in-depth characterization of the glacial ice, and on novel approaches in event reconstruction that utilize fast approximations of photoelectron yields. Here, a more accurate model is derived for event reconstruction that better captures our current knowledge of ice optical properties. When evaluated on a Monte Carlo simulation set, the median angular resolution for in-ice particle showers improves by over a factor of three compared to a reconstruction based on a simplified model of the ice. The most substantial improvement is obtained when including effects of birefringence due to the polycrystalline structure of the ice. When evaluated on data classified as particle showers in the high-energy starting events sample, a significantly improved description of the events is observed.


Introduction
The IceCube Neutrino Observatory detects neutrinos interacting with nucleons and electrons in the South Pole ice via Cherenkov radiation produced by charged secondaries.The ice is instrumented with 5,160 digital optical modules (DOMs), each with a single downward-facing photomultiplier tube (PMT), arrayed across a cubic-kilometer volume below the surface.The DOMs are attached to 86 strings-cables installed in the ice that provide mechanical and electrical support.The DOMs are vertically spaced 17 m apart on standard IceCube strings and 7 m apart on DeepCore strings, a denser infill region of the detector.Standard IceCube strings are spaced approximately 125 m apart [1].
IceCube covers a rich and diverse physics program at a wide energy range, and can detect neutrinos with energies spanning from about 5 GeV to above 10 PeV, as well as bursts of MeV neutrinos from sufficiently nearby sources.Highlights include the discovery and subsequent confirmation of a diffuse flux of astrophysical neutrinos [2,3], the first identification of an astrophysical neutrino source, TXS 0506+056 [4], arising from IceCube's real-time program [5], and the detection of the first astrophysical neutrino interaction at the Glashow resonance [6].More recently, the nearby Seyfert galaxy NGC 1068 [7] and the Milky Way itself [8] have been identified as steady sources of astrophysical neutrinos.These and future results rely on, and will continue to benefit from, refined calibration of the detector and improved event reconstruction.
In IceCube there are, broadly speaking, three general event categories: tracks, cascades, and double cascades.In this paper, these terms refer to the hypothesized model underlying event reconstructions, which is related to but conceptually distinct from the actual physical process of particle showers or muon propagation through matter [9].High-energy muons produced in muon neutrino charged-current (CC) interactions can travel large distances through the detector and appear track-like.Electromagnetic (EM) particle showers induced by electron neutrino CC interactions or hadronic showers from neutral-current (NC) interactions of all flavors appear as cascades, a roughly spherical deposition of light in the detector.Tau neutrino CC interactions can appear as a double cascade for certain tau lepton decay channels and energies.For the purposes of event reconstruction, high-energy muons can be approximated as a series of cascades due to muon stochastic losses [10], and other events can be approximated using a single-or two-cascade model [11].
The glacial ice serves as a natural Cherenkov medium for detection of the byproducts of neutrino interactions.From the surface to the bedrock, ice isochrons have formed over geological time scales, each with different scattering and absorption lengths that affect the propagation of Cherenkov photons.Calibration LEDs on board each DOM have been used to provide a detailed measurement of these optical properties as a function of depth [12].The propagation length for a typical Cherenkov photon at 400 nm is shown in red in Fig. 1, and illustrates the depth-dependent optical properties of the natural ice sheet.The propagation length is defined in Ref. [13], and is proportional to the geometric mean of the absorption length and effective scattering length.For this illustration, perfectly flat ice isochrons are assumed and no additional ice anisotropies are present besides the depth-dependent scattering and absorption.In reality, we now know that the ice is more complex and exhibits additional anisotropies.These include an axial dependence now attributed to birefrigence of the polycrystals in the ice, which leads to a higher photon yield along the axis of glacial ice flow [14].Elevation changes of ice isochrons, or ice layer undulations, have also been mapped out in detail throughout the detector volume [15].Inclusion of these effects into a model usable for event reconstruction is our primary focus.
Accurate event reconstruction requires accurate modeling of photoelectron arrival time distributions for a given physics hypothesis.In this paper we focus on energies most relevant for high-energy astrophysical neutrinos, above roughly 10 TeV.As an example, Fig. 1 provides a visualization of the expected Cherenkov photoelectron yield arising from an EM shower.The shower is placed at various depths in the detector,  s , with  s = 0 corresponding to a depth of 1,948 m below the surface, arriving from a direction Θ s = 90 • relative to vertical.Photoelectron yields are shown for receivers, which are idealized representations of IceCube's downward facing PMTs, placed at distances 100 m in front of (black, solid) and behind (black, dashed) the shower.Note the impact of the propagation length (red) on the photoelecton yields.Differences between yields at the two receivers illustrate how shower directionality can be reconstructed, while the strong correlation to photon propagation length highlights the importance of accurate ice modeling.
The rest of this paper is organized as follows.In Section 2, we provide an overview of some of the challenging aspects of shower reconstruction in IceCube.Section 3 describes improvements when fitting tabulated photoelectron distributions from Monte Carlo (MC) simulations with tensor-product B-splines, which reduce observed artifacts in the reconstructed zenith distribution while improving The solid (dashed) black line indicates the number of expected photoelectrons for receivers placed 100 m in front of (behind) the shower.The red line shows the depth dependence of a 400 nm photon's propagation length (right axis), which is proportional to the geometric mean of the absorption length and effective scattering length as defined in Ref. [13].Since the ice is a natural medium formed over epochs, large variations in optical properties that appear over large commensurate time periods are observed [16].A region of heightened dust concentration is highlighted in gray.The instrumented depths that IceCube spans are highlighted in blue, and include the dust layer.Clear differences for the two receivers illustrate how shower directionality can be reconstructed, while the strong correlation between yield and propagation length highlights the importance of accurate ice modeling.
the median angular resolution for a benchmark simulation set by about 1 • at 1 PeV from an original resolution of 12 • .Section 4 details an effective treatment to include ice birefringence without increasing the dimensionality of the model used in reconstruction, further lowering the median angular resolution by 5.3 • at 1 PeV.Section 5 describes a correction to account for ice layer undulations that yields another 1.7 • improvement at 1 PeV.Section 6 touches on an approximation of shower longitudinal extensions with a two-cascade model, which yields an additional 0.5 • improvement at 1 PeV, resulting in a final median angular resolution of 3.5 • , and includes a brief discussion of systematic uncertainties.In Section 7, we provide a summary of the results and conclude the discussion.

Challenges in shower reconstruction
Figure 2 shows an event view of a high-energy shower in IceCube.For events contained within a fiducial volume such as this, the energy reconstruction is relatively well constrained by calorimetry [10].Further, at the energies relevant in IceCube, particles produced by neutrino interactions in the ice are, to good approximation, colinear with the neutrino direction, and thus their reconstruction points to the neutrino arrival direction.However, due to the large length scales between IceCube sensors, the directional reconstruction of showers can be a challenge and depends critically on the ice model.Anisotropies in the ice can induce differences in detected Cherenkov photoelectron yields that mimic those caused by a shift in direction.
The most realistic model of photon propagation in the ice requires a full MC simulation.The parallel nature of photon propagation makes such a task well suited for graphics processing units (GPUs) [17,18].Such an approach has allowed for a refined description of the optical properties of the ice by simulating LED calibration devices and fitting to the observed in situ calibration data [12,14].Physics simulations of particle showers also rely on MC-based photon propagation.
Processing speeds on GPUs may be sufficient for simulation, but event reconstruction often requires testing orders of magnitude more points in the physics parameter space for each event.As of this writing, utilizing GPUs for event reconstructions where each Cherenkov photon is fully resimulated for every tested hypothesis is only feasible for  (100) events [19,20].At the expense of some accuracy, fast approximations of Cherenkov photon yields have been developed to address this.Initial approximations relied on look-up tables [21], which were improved to allow for gradient-based minimization using a tensor product of B-splines [10,22].More recently, neural network (NN) [23][24][25] models have been employed as well.All such models are approximations, being derived from full MC simulations of particle showers or minimum ionizing muons.
There are then two main challenges that arise in event reconstruction.The first is to accurately model the physics and instrumentation to the best extent possible within the full MC simulation chain.To a large extent this includes the modeling of the optical properties of the ice, where significant progress has been made over the past few years that led to an improved description of the observed ice anisotropy based on birefrigence [14] and a more realistic mapping of layer undulations [15].Including these effects has resulted in much better agreement with in situ calibration data.
The second challenge is to quickly and accurately model the photoelectron yields derived from full MC simulation.Anisotropies in the ice due to birefringence and layer undulations bring further complexity by breaking symmetry and thus introducing additional dimensionality.A full parametrization of the photoelectron yield for a cascade would span nine dimensions: ( s ,  s ,  s ) to define its position, (Θ s , Φ s ) to define its arrival direction, (, , ) to define the DOM position relative to the cascade, and , the photon arrival time at the DOM.This increased complexity makes a tensor product of B-splines intractable to evaluate over all relevant regions of the parameter space; memory requirements over the full nine dimensions would exceed 500 TB.In addition, the simulation needed to construct approximators becomes prohibitively expensive to produce even on the Open Science Grid [26], likely exceeding 10 10 CPU hours.While NNs are not bound by these dimensionality constraints, they can be subject to unknown inaccuracies in interpolation and extrapolation while also requiring a large training dataset.Depending on the approach, for example if a fixed set of geometry, ice properties and detector settings are assumed in the simulation, the model might need to be retrained for each successive update [23][24][25].Furthermore, training an NN may provide less insight into where limitations in simulation and reconstruction lie.If MC simulation is taken as an absolute ground truth, then we are subject to any existing mismodeling there.A more iterative approach, with cross checking between simulation and expected results from event reconstruction, can help resolve issues that exist either in the MC or in the construction of any approximate model.

Impact of ice anisotropies on shower reconstruction
To confront the asymmetries introduced by a refined understanding of the ice, a set of corrections to approximate birefringence and layer undulations can be employed.As will be shown, these corrections do not increase the dimensionality of the model; six dimensions is enough, with (Θ s ,  s ) to fix the zenith angle and depth of the shower and (, , , ) to define the arrival position and time at a receiving DOM [21].The impact of these corrections can be evaluated on a benchmark, high-energy starting events (HESE) MC set, consisting of simulated showers from electron-neutrino and electron-antineutrino interactions contained within a fiducial region of the detector [3].Electron (anti)neutrinos are sampled from an  −1.5 spectrum at the Earth's surface, propagated through the Earth to reach IceCube, whereupon they may interact and produce signatures that pass the HESE selection criteria.After cuts, the events range in energies from 10 TeV to 10 PeV, a regime where the astrophysical neutrino flux is expected to transition above the atmospheric neutrino background.The simulation relies on GPUs for photon propagation [18] and incorporates our most up-to-date knowledge of the ice, namely birefringence [14] and recent updates to ice layer undulations [15].The left panel of Fig. 3 shows the distribution of the angle between the true direction and reconstructed direction, Ψ.Including the birefringence correction (green, dashed) improves upon the event reconstruction that employs no corrections (red, dotted).Including both birefrigence and layer undulation corrections (orange) exhibits even better performance.The implementation of these corrections is the subject of Section 4 and Section 5. ι (1) − ι (2)   Figure 3.These distributions highlight the improvements in Ψ and  when corrections that approximate the discussed ice anisotropies are included in reconstruction.The left panel compares Ψ distributions with no corrections (red, dotted), with birefringence (green, dashed) and with both layer undulation and birefrigence (orange, solid).Each correction successively reduces the space angle between the true and reconstructed direction.The right panel shows two distributions across Δ relative to a model that includes all corrections.Lower  values correspond to improved agreement with simulated events in the benchmark MC.The reduced negative log-likelihood for a model that includes corrections due to both layer undulations and birefringence is given by  (2) , including only the birefringence correction  (1) , and including no corrections  (0) .The improved description with  (2) is evidenced by the skew towards positive values, which is more apparent for  (0) −  (2)  (red, solid) but still visible for  (1) −  (2) (green, dashed).The vertical dashed line at Δ = 0 serves as a visual guide to highlight the positive skew.
Another metric for comparison is the negatived log-likelihood per degree of freedom (reduced negative log-likelihood), , obtained by comparing expected photoelectron yields and time profiles from the best-fit cascade to the simulated data.Here we use an effective Poisson-based likelihood [27], modified to use a fixed relative uncertainty across all bins to capture residual error in the model.The right panel of Fig. 3 illustrates the  improvement when all corrections are applied.Shown are two Δ distributions,  (1) −  (2) and  (0) −  (2) , where  (0) ,  (1) and  (2) refer to reduced negative log-likelihoods for event reconstructions that assume no corrections, birefrigence correction only, and all ice-associated corrections, respectively.As lower values correspond to an improved event description compared to the simulated data, the positive skew observed for both distributions indicates the improvement of  (2) over  (1) and  (0) .

Improvements in model construction
As mentioned in Section 2, the photoelectron-yield model that is the focus of this paper derives from GPU-based MC photon propagation [18].The convention used to describe source and receiver coordinates, introduced in Ref. [21], is shown in Fig. 4. For directional sources of Cherenkov photons, a simplified MC where birefringence (see Section 4) and ice layer undulations (see Section 5) are not included is used to tabulate yields, so that azimuthal symmetry is preserved about the ẑ axis and translational symmetry is preserved over ( s ,  s ).Further, in ice with depth-dependent optical properties (see Fig. 1) the axial symmetry of Cherenkov photon emission about the shower axis, p, which holds to good approximation for the colinear particles produced in high-energy particle showers, is reduced to a bilateral symmetry where  is degenerate with −.Thus, the source can be completely described by (Θ s ,  s ), and the receiver  ranges from 0 The convention used to describe source and receiver coordinates [21].Source and receiver positions are indicated as filled and empty circles, respectively.The unit vector p denotes the shower momentum direction, while (Θ s , Φ s ) describes the arrival direction, − p, in the standard spherical coordinate system with polar axis aligned with ẑ.The receiver coordinates, (, , ), are given in a spherical coordinate system centered on the source position and defined such that the polar axis aligns with p and  = 0 corresponds to where the projection of r onto the plane perpendicular to p (shaded) is maximally vertical.In a media with symmetry about the ẑ axis, Φ s is degenerate.If, additionally, photon emission is axially symmetric about p,  is degenerate with −.
A 1 GeV EM shower, the photon source, is repeatedly simulated along a grid of points, (Θ s, ,  s,  ), indexed by  and .For each source, the photoelectron yield for receiver  at a location r  = (  ,   ,   ) relative to the source is tabulated as a function of where  is the group index of refraction in ice and  is the absolute photon arrival time [21].The receiver position is defined using a spherical coordinate system aligned with the shower principal axis, as shown in Fig. 4. Once tabulated, the total number of photoelectrons, or amplitude, as well as the cumulative density function (CDF) in  res is fitted across (r  , [ res ]) for each simulated cascade in (Θ s, ,  s,  ) using a tensor product of B-splines.Then each fitted surface over the receiver coordinates is stacked across (Θ s, ,  s,  ) [22], resulting in a function of the form where  is the amplitude and  is the CDF.As the number of Cherenkov photons scales linearly with shower energy [10], Eq. (3.2) gives, within a constant factor, the number of photoelectrons detected by a DOM for an EM shower.

Interpolation
The first step towards improved modeling for shower reconstruction required resolving a long-standing issue that was observed in the reconstructed Θs distribution.Figure 5 illustrates the problem for our benchmark MC.The right panel shows MC statistics after applying HESE cuts [3] as a function of the true EM-equivalent deposited energy,  dep .Sometimes referred to as visible energy [10],  dep is calculated for each simulated event by summing over the energies of all EM and hadronic shower components, with a correction that scales down the true energy of hadronic showers to match the energy of an EM shower that, on average, produces an equivalent number of Cherenkov photons.Gray dashed and dotted lines show the breakdown for upgoing and downgoing events, respectively.The left panel compares the true cos Θ s distribution (black) against the reconstructed distribution using a now-outdated model (gray) [10], which was derived from the best knowledge of the ice in 2013 [12].Following the terminology in Ref. [12], we refer to this as the "Mie model".At the time, ice layers in the Mie model were assumed to be flat for purposes of event reconstruction, and no anisotropies beyond the depth dependence shown in Fig. 1 were included.In the left panel, dashed vertical lines indicate B-spline knot positions in cos Θ s used in the Mie model, and a clear ringing pattern arising from an artificial pull away from knot locations emerges.The solution came in the form of improved interpolation with additional support points, the topic of this section and illustrated by the cos Θs distribution in red.Multidimensional B-spline surfaces constructed with a tensor product can have interpolation artifacts along diagonals across coordinate axes.This is due to defining basis functions on Cartesian grids [28].In the absence of dense support points, smooth features along diagonals are not preserved and interpolation exhibits an oscillatory behavior.In particular, this effect led to the ringing effect shown by the gray distribution in the left panel of Fig. 5. Including the updates described in this section results in the reconstructed cos Θs distribution shown in red, without ringing and in better agreement with the true cos Θ s distribution (black).The birefrigence and layer undulation corrections, which will be discussed in Section 4 and Section 5 respectively, were not included in the updated model (red) shown here.Including these effects would further improve agreement to the true distribution in the most upgoing region (cos Θ s ≈ −1).
Figure 6, left panel, illustrates the interpolation artifacts observed in the amplitude, (r, Θ s ,  s ), for the Mie model at the coordinates given in the title.In this case,  s = 30 m corresponds to a source placed in a region of clear ice (see Fig. 1).As shown in Fig. 4, the receiver coordinate system is spherical and aligned with respect to Θ s .For example, at Θ s = 90 • a receiver at  = 0 • would be placed directly in front of the shower and hence in the same ice layer.As Θ s sweeps through its parameter space, a compensation in  of equal magnitude ensures that the receiver is placed at the same depth as the source.The result is that for clear ice, the amplitude will be greater along that diagonal than away from it, where a photon would have traversed through regions of higher absorption and scattering.However, due to the low number of support points in Θ s (every 10 • ) over which the model is stacked, interpolation artifacts appear.

N
All events cos Θ s < 0 cos Θ s ≥ 0 Figure 5.The left panel shows the number of events per bin, , distributed over cos Θ s for the benchmark MC sample consisting of electron neutrino and antineutrino interactions contained within a fiducial volume [3].The effects of birefrigence [14] and ice layer undulations [15] are included in the simulation, but are not included in either of the event reconstructions shown.A ringing effect is visible when the Mie model is used (gray), which disappears when the updates discussed in Section 3 are included (red).Vertical dashed lines indicate the placement of knot positions for the Mie B-splines along the cos Θ s dimension.The MC truth (black) is shown for comparison.To give a sense of the statistics and energies of this MC set, the right panel shows the overall  dep distribution in black.Gray dashed and dotted lines show the breakdown for upgoing and downgoing events, respectively.See text for more details.
A simple, but effective, solution was to double the number of support points in Θ s , tabulating at 5 • intervals instead of 10 • .The number of knot locations in Θ s roughly doubles as well, the cause of which is described in Ref. [22].With this change, the amplitude model yields smoother behavior along diagonals, as shown in the right panel of Fig. 6.Updated values of the depth-dependent ice optical properties are included in the simulation used to fit the model [15].However, similar to the construction of the Mie model, the effects of birefringence and ice layer undulations are turned off for reasons of simplification discussed in Section 2 and at the beginning of this section.In both panels, the dashed white lines bound regions in the phase space where receiver locations fall within ±10 m of  s in depth.Note that while the visualization here is for the amplitude model, the denser support points are applied to the CDF model in  res as well.
Figure 6 only presents a single two-dimensional slice; other parameters are fixed at the coordinates indicated in the title.The full effect across all dimensions is more complex, and an example of what happens when  = 90 • is shown in Fig. 7.This also explains why, although a simple coordinate transformation could resolve the issues in the slice shown in Fig. 6, a generalized coordinate transformation over the additional  dimension does not exist in simple form and hence we resorted to increasing the number of support points in Θ s .

Extrapolation
Extrapolation failures can also occur when a parameter's domain extends beyond its last support point.These failures led to sharp discontinuities in the photoelectron yields across boundaries, and they occurred in the spherical coordinates system defined by  and , near the poles where cos  = ±1 and near the boundary of  ∈ {0 • , 180 • }.The former required the construction of additional data points at cos  = ±1 based on linear interpolation, and the latter was resolved by reflection across the boundary in .By including these additional support points, extrapolation artifacts were largely reduced.
The discontinuity at cos  = ±1 is again related to the tensor product construction of the model, which does not map naturally onto a spherical coordinate system.Simulation data is linearly binned in cos , with the abscissa of each data point taken at the bin centers.Thus, at exactly cos  = ±1 no data point can be tabulated.To reach cos  = ±1, the Mie model extrapolated beyond the closest bin center.Furthermore, at cos  = ±1 the model would extrapolate to different values as a function of , leading to the discontinuity visible in the gray lines in Fig. 8, which shows the time integrated amplitude as a function of Θ s for two receivers placed at Θ r = 60 • and 80 • .Here, we have chosen r to lie in the same plane as p and ẑ for simplification, but these examples are representative of behavior across other dimensional slices.As shown in Fig. 4, when all three vectors lie in the same In a spherical coordinate system, the azimuthal angle  becomes degenerate at the poles, and there the model prediction should not change as a function of .Unfortunately, this cannot be guaranteed in the tensor product B-spline construction.The extrapolation in the Mie model, combined with the flip in  when  = 0 • , cause the discontinuities in the gray lines in Fig. 8.The discontinuities can be largely reduced by the addition of a set of support points at cos  = ±1 such that the to-be-fitted data is the same across all  at the poles.These data values were computed by first performing a linear interpolation of the tabulated data to  = 90 • , then a linear extrapolation to cos  = ±1 along the  = 90 • curve.By including these values in the B-spline fit, the polar discontinuity was reduced, as shown in the red lines in Fig. 8.
A final fix was applied in the  dimension when extrapolating towards its bounds.Equally spaced bins ranging from 0 • to 180 • are constructed along , and receiver coordinates   lie at the bin centers.The Mie model extrapolated beyond the first and last bin center to reach  = 0 • and 180 • .A simple solution to improve modeling near these boundaries was to reflect the nearest four data points across the boundary, thereby extending the support points in  to beyond its domain, and perform the B-spline fit with those additional data points.This meant that the previous extrapolation to the bounds was replaced by interpolation between well-defined, tabulated values.
Figure 9 illustrates the impact on the angular resolution going from the Mie model that is subject to interpolation and extrapolation artifacts (gray) to a newer ice model that includes the fixes described in this section.The left panel shows quartiles (25-50-75 percentile levels) of Ψ distributions for the benchmark MC as a function of  dep , obtained using the binning shown in the right panel of Fig. 5.The right panel shows distributions of Ψ in two different energy slices as described in the legend.We see that, in addition to the much improved zenith distributions shown in Fig. 5, the angular resolution is consistently worse with the Mie model, which exhibits additional degradation going towards higher energies.

Approximation of ice anisotropies due to birefringence
Accurate characterization of physics quantities across a sparse array of PMTs requires accurate calibration of ice and instrument.This has been accomplished with calibration LEDs attached to a dedicated flasher board on each IceCube DOM [12] and by using large samples of downgoing minimum ionizing muons to set the global energy scale [10].The propagation of light in the glacial ice sheet is complex at Cherenkov wavelengths, with dependencies on depth that reflect Earth's climate across geological time scales as shown in Fig. 1.Absorption and scattering coefficients that describe the mean free path of a photon are fitted using calibration LED data over a range -12 - that encompasses all instrumented depths [12,14].In general, the Cherenkov photon absorption (scattering) length in clear regions of the South Pole ice is longer (shorter) than in natural bodies of liquid water where neutrino detectors have been constructed or proposed [29].
IceCube also discovered a directional dependence in the light propagation, or ice anisotropy, with maximal effect along the ice flow axis [30].The flow axis is the direction that the bulk ice moves, which is at a rate of about 10 m/yr [31].A microscopic explanation of the ice anisotropy due to birefringence of polycrystals was given in Ref. [14,32] and leads to improved agreement with calibration data over the previous phenomenological model [30].However, as mentioned in Section 2, this breaks the azimuthal symmetry and naively would require another dimension, Φ s , to be included in the model.
An alternative approach is to separate ice anisotropies from shower physics and take advantage of the symmetries inherent in both.The photoelectron yields from a particle shower can be evaluated in a bulk ice with no birefringence and flat, horizontal ice layers thus preserving azimuthal symmetry [21].The impact of ice birefringence can be estimated by simulating an isotropic, point-like light source and tabulating and fitting the corresponding photoelectron yields using two different ice models, one with birefringence and one without, including the improvements described in Section 3.2.An isotropic source simplifies the problem by removing the dependence on Θ s .Comparisons of the two in amplitude, as done in Refs.[33,34], and in shape, as discussed in Section 4.2, allow for simple coordinate transformations that encode the impact of birefringence on Cherenkov photons.
We emphasize that the key idea is the switch from a directional source to an isotropic source and that the birefringence corrections are computed using an isotropic source, which are later applied to directional sources.This is an approximation, but one which substantially improves the model and the angular resolution of showers, as will be shown.

Amplitude
The time-integrated photoelectron yield, or amplitude, at a receiver DOM differs depending on whether the ice is modeled with or without birefringence.Since the amplitude decreases monotonically with , a translation to a larger (smaller) effective distance,   , can correct for too high (low) amplitude observed in the simplified ice model.Birefringence applies globally across the detector, independent of ( s ,  s ) [14].Due to layer-by-layer differences in the ice, the  s -dependence is kept.The goal, then, is to construct a function   (r,  s ) that can be applied to correct amplitudes in the simplified model.Two simulation sets are produced.Each consists of identical, isotropic, point-like light sources placed at a series of depths,  * s,  , where  is an index and the asterisk distinguishes this as referring to an isotropic source.One set assumes the simplified ice model, which does not include birefringence or layer undulations.The other includes birefringence but no layer undulations (see Section 5).The time-integrated amplitudes are tabulated at receivers placed in spherical coordinates r  = (  ,   ,   ) from the isotropic source, with the polar angle  now defined relative to ẑ.A tensor-product B-spline is then fit to the tabulated data in order to construct a smooth approximator for the isotropic source Once evaluated across the grid of receiver positions, the discrete set of   (r  ,  * s,  ) are fitted again with tensor-product B-splines for each  * s,  and then stacked across  * s,  to give   (r,  * s ).The obtained   is presented for a slice across cos  = 0 and  * s = 0 m (center of IceCube) in the left panel of Fig. 10.Two circles are drawn as solid lines at  = 60 m (blue) and 120 m (orange).The corresponding   at those distances are shown as dashed lines in blue and orange, respectively.The direction of ice flow is indicated by the arrow.The right panel includes a similar visualization, only now in two planes that lie perpendicular to the horizontal plane shown in the left panel.The dashed (dotted) line shows   in the vertical plane that lies along (perpendicular to) the flow axis, at [14].Thus, the right panel displays two -planes, where  is the distance along the direction given by û = cos  x + sin  ŷ.The effect of birefringence is to shrink the distance along the ice flow axis, leading to a larger amplitude, while slightly increasing the distance perpendicular to the flow axis, leading to a smaller amplitude.These effects are consistent with observations from calibration data [14].

Time probability density function
The   discussed in Section 4.1 corrects for amplitude differences due to birefringence but does not account for differences in the time profile.These photon arrival time distributions, i.e. the probability density function (PDF) in  res , add timing information to the model (see Eq. (3.2)).
The impact of birefringence on the PDF is shown in Fig. 11 for two azimuthal directions as given in the titles.Their relative positions to the source are indicated by the star and dot in the left panel of Fig. 10.The impact of birefringence on the time PDF is illustrated going from the red, dashed to the black, solid lines, with a larger difference visible in the left panel, when the source-receiver direction lies along the flow axis.Notably, along the ice flow axis, the time PDF is squeezed narrower while preserving the mode of the distribution.It is not possible to capture such an effect with a sole translation in .In order to obtain a narrower distribution, a smaller source-receiver distance is needed, which in turn causes a shift in the mode.
A natural generalization to model shape differences is to introduce a translation in  in addition to .That is, instead of a single-parameter correction   (r,  * s ) we seek a two-parameter correction   (r,  * s ) and   (r,  * s ), which can be applied simultaneously to better approximate the  res PDF.Using the same two simulation sets described in Section 2.1, and extending the construction to include timing information, we generate CDFs  bfr (r,  res ,  * s,  ) and  simple (r,  res ,  * s,  ).The corrections are obtained by Nelder-Mead [35]   fit to optimally describe the PDF over the region where most of the density is concentrated.The result is illustrated by the green lines in Fig. 11, which show  simple (  ,  res −   ,  * s,  = 0) at the best-fit values for   and   obtained by minimizing Eq. (4.2), for the receiver coordinates indicated in the panel titles.While slight differences still exist compared to the  res PDFs with birefringence (black lines), the agreement is significantly improved over the raw  res PDF without birefringence (red, dashed lines).For each  * s,  , the obtained discrete sets of   (r  ,  * s,  ) and   (r  ,  * s,  ) are separately fitted across receivers  with tensor-product B-splines.Finally, these surfaces are stacked across  * s,  to give   (r,  * s ) and   (r,  * s ).The functions   (r,  * s ),   (r,  * s ) and   (r,  * s ) can be jointly used to approximate the effects of ice birefringence on shower photoelectron yields and time profiles.In order to substitute  * s →  s , we need to assume that the corrections obtained with an isotropic source generalize to a directional source.To justify this assumption, note that these transformations are translations, rather than rotations.A modification of the radial and time coordinates passed to the amplitude and CDF given in Eq. (3.2) does not destroy information about the directionality of Cherenkov photon emission from particle showers.Thus, Eq. (3.2) becomes Here we have used the notation r  to refer to the source-receiver vector in the standard spherical coordinate system with polar axis ẑ, in contrast to (, ), which are defined in the spherical coordinate system with polar axis p as shown in Fig. 4. Finally, since   ,   and   are differentiable functions, the Jacobian can be computed for gradient-based minimization.
-16 -  9 but now comparing the angular resolution using Eq.(4.3) (green, solid) to Eq. (3.2) (red, dashed), where both include the fixes described in Section 3. In the left panel, a substantial improvement is observed when birefringence corrections are included.Note the red intervals are identical to those in Fig. 9.The right panels show distributions of Ψ in two different energy slices, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.
To validate that the corrections described in this section lead to an improved description of particle showers in IceCube, we compared the directional reconstruction performance of the benchmark MC set introduced in Section 2.1.Figure 12 shows the angular resolution, Ψ, with (green, solid) and without (red, dashed) birefringence correction.The left panel shows the median angular resolution (solid and dashed lines) as a function of  dep .The bands correspond to the 25 % and 75 % intervals in the Ψ distribution.The right panels show the full Ψ distribution in two energy regimes, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.From these comparisons, it is evident that a substantial improvement in the directional reconstruction is achieved when using Eq.(4.3) over Eq. (3.2), one which far exceeds the improvement from additional photon statistics going from lower to higher energies.

Approximation of ice layer undulations
Polar ice stratigraphy exhibits a strong dependence on depth, broadly reflecting changes in Earth's climate over time.Layers of relatively homogeneous ice isochrons-ice layers that formed at around the same time period-are each modeled with unique scattering and absorption coefficients.The elevation change for each layer was initially calibrated using data from a laser dust logger deployed in seven IceCube boreholes [16].More recently, a reevaluation of ice optical properties with in situ calibration data revealed that ice layers exhibit undulations, believed to be formed by geological structures in the bedrock, across the detector array rather than solely along a single axis [15].As layer-by-layer ice property differences can result in significant differences in the detected photoelectron yields (see Fig. 1), an accurate description is important for simulation and reconstruction.

Parameterization of the ice depth dependence
Since ice isochrons are not perfectly flat, optical properties of the ice are defined as a function of depth (or  in detector coordinates) at a fixed reference point, P, near the origin in the -plane and shown as the red dot in the left panel of Fig. 13 [12,15].In other words, the ice model is parameterized as a function of what we will refer to as the P-depth.A mapping of ice layer undulations is needed to convert a physical position, (, , ), back to the appropriate P-depth,  P (, , ), to retrieve the corresponding ice optical properties at that position.Such a mapping can then be used to apply the correct scattering and absorption lengths for any position throughout the detector.The left panel of Fig. 13 shows the ice layer elevation change across the -plane for the isochron at a P-depth of 2,248 m ( P = −300 m in detector coordinates).Specifically, the color map shows instead to obtain  P ( s ,  s ,  s ) =  P (X s ).Substituting into Eq.( 4.3) we obtain    r  ,  P (X s ) , , , Θ s ,  P (X s )    r  ,  P (X s ) , , ,  res −   r  ,  P (X s ) , Θ s ,  P (X s ) , (5.1) which gives the functional form of the photoelectron yield correcting for both birefringence and ice layer undulations.As  P is a function constructed by linear interpolation between fixed points in position space, it is not smooth everywhere [15].In practice, this does not pose an issue for computing gradients, since discontinuities in its derivatives are small and only rarely occur.Thus, we can compute the Jacobian terms  P /   , where   ∈ { s ,  s ,  s }, for use in gradient-based minimizers.
Factorizing out ice layer undulations from the B-spline model itself makes it extremely simple to switch to newer  P (, , ) models as well.If, instead, one were to generate photon-yield expectations without this factorization, then for each update to  P (, , ) a massive resimulation campaign would be needed to generate the raw data to construct an updated model.The dimensionality would also need to be increased to include ( s ,  s ).The cost of the factorization is that using only the source position of the cascade is an approximation, as the path between source and receiver DOM can traverse several layers of ice with varying amounts of undulations.This can be seen in the right panel of Fig. 13.Fortunately, the ice layer undulations are gradual and this approximation improves for shorter source-receiver distances, which is also where most of the statistics used in event reconstruction are expected.To validate Eq. (5.1), we again use the same benchmark MC as described in Section 2.1.The impact of ice layers is easily seen in reconstructed zenith, Θs , as a function of Ẑs .Figure 14  Using the benchmark MC, the left panel shows the distribution of the reconstructed distance between two cascades for three different energy regimes.There is a trend towards slightly longer distances as the true EM-equivalent deposited energy,  dep , increases, which follows expectations from shower physics [36].The right panel shows distributions of the pulls on reconstructed energy, denoted as Êdep , for the same energy slices as in the left panel, illustrating that its resolution is on the order of a few percent.two-cascade routine.As a final step, the likelihood of both single-and two-cascade reconstructions are compared and best fit is chosen.In the vast majority of cases, the two-cascade routine is preferred.The benchmark MC introduced in Section 2.1, which, like all standard IceCube MC, includes a simulation of the average longitudinal profile of particle showers [36] but not their shower-to-shower fluctuations aside from those due to the Landau-Pomeranchuk-Migdal effect [9,37] for EM showers with energies above 1 PeV, is used for validation.Figure 16 shows cross checks of the reconstructed separation distance distributions (left panel) and pulls on  dep , the EM-equivalent deposited energy (right panel).In both panels, three different  dep slices are shown.The reconstructed distance distribution is pulled to larger values with increasing  dep , though all peak below 10 m, consistent with expectations from   simulations.
Figure 17 shows the improvement in the angular resolution of the benchmark MC.Distributions of Ψ with the two-cascade, single-cascade and original Mie model are shown in solid blue, dashed orange and dotted gray, respectively.The left panel shows the median angular resolution (solid, dashed and dotted lines) as a function of  dep .Again, bands correspond to the 25 % and 75 % intervals in the Ψ distribution.The right panels show the full Ψ distribution in two energy regimes, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.We see that including a two-cascade reconstruction leads to some further improvement in the directional reconstruction, since it better describes the longitudinal extension of particle showers.It is worth emphasizing that the corrections described in Sections 3 to 5 are prerequisite; a two-cascade reconstruction exhibits improvements only when an accurate single cascade model is used.

Ice systematic uncertainties
One ice systematic that can affect shower directional reconstruction is the bubble column, or hole ice, that formed as part of the drill-hole refreezing process [38].Based on camera footage and in situ calibration data, the hole ice is known as a centrally located region of heightened scattering and absorption.As its optical properties are less understood than the bulk ice, it has traditionally been modeled as a global modification in the DOM acceptance as a function of the incident photon direction [38].When the forward scattering region is strongly modified, a degradation of the angular resolution is observed on the order of 0.5 • to 1 • .As the angular sensitivity curves modify photon acceptance along the (downward-facing) PMT axis, mismodeling of the hole ice can also pull the reconstructed Θs by up to 2 • to 3 • for events arriving horizontally, with smaller pulls for non-horizontal events.
The optimal results obtained in this work, shown in Fig. 17, rely on B-spline surfaces fitted to an ice model similar to the one used in the benchmark simulation; the hole ice model is identical, and only minor differences-in bulk ice optical properties and layer undulations-exist between the model used to construct Eq. (5.1) and that used in the MC.Exclusive of hole ice, a variance of the other ice optical properties at the percent level, which is on the order of the current uncertainty envelope, was found to have a negligible impact on angular resolution.
Further, there are effects for which robust quantification of uncertainty currently does not exist, such as birefringence and the extrapolated layer undulation.Given the accuracy with which calibration data is now described [14], any residual systematic mismodeling of the anisotropy along the ice flow axis is likely to be much smaller than what is shown in Fig. 12. Mismodeling of ice layer undulations outside the instrumented footprint could also have an impact on the reconstruction of events that occurred outside the detector.These events tend to be more difficult to reconstruct, with worsening resolution as the interaction vertex moves further away from the instrumented region due to increasingly limited arrival photon statistics, and inline with what might be expected of the accuracy of the extrapolation itself.Thus, the extrapolation systematic should be a subdominant effect.Finally, there are potentially undiscovered systematics, the impact of which cannot be evaluated based on simulation alone and may require data-driven approaches to quantify.Such effects lie beyond the scope of this work, though the general techniques described here should be applicable as future improvements are discovered.

Summary
Much progress has been made over the past few years towards an improved understanding of the South Pole ice, including a microscopic description of the observed anisotropy along the ice flow axis attributed to polycrystalline birefrigence [14], and a detailed mapping of ice isochron undulations across the detector [15].These refinements to the ice model can be incorporated into the reconstruction of in-ice particle showers, either with neural networks [23][24][25] or via a series of physically motivated corrections, as discussed in Section 4 and Section 5. Using Eq. (5.1), the shower longitudinal extension can be approximated with a two-cascade model, as highlighted in Section 6.1.Table 1 gives an overview of the median angular resolutions obtained for our benchmark MC, which is simulated with recent updates to ice modeling that include birefringence [14] and ice layer undulations [15], at four different energies as cumulative model improvements are included in reconstruction.At energies above 1 PeV, a 3.5 • median angular resolution is achieved when all corrections are included, which is over a factor of three improvement compared to using only Eq. (3.2).
Table 1.The median angular resolution at the energies listed for the given models.A detailed description of improvements in model construction can be found in Section 3, the birefringence correction in Section 4, the ice layer undulations in Section 5, and the approximation of shower extension in Section 6.1.Compared to models that do not include any corrections, the angular resolution is improved by more than a factor of three at 1 PeV.Ice model systematics, in particular bubble columns that formed as part of the drill-hole refreezing process [38], can lead to some degradation in resolution as discussed in Section 6.2.For the variations tested, a worsening in the median angular resolution of up to 1 • is seen when differences in the forward-scattering region are large.A bias in the reconstructed zenith distribution is also observed.Note, however, that with smoother B-spline surfaces, the more artificial zenith biases have been mitigated using the methods discussed in Section 3. Additionally, the IceCube Upgrade, which will be a dense infill extension at the center of IceCube, is currently planned for installation in the next couple years [39].In addition to new optical modules, next-generation calibration devices will be deployed that should provide additional information about the optical properties of the bubble columns [40].ι (1) − ι (3)   Figure 18.Distributions of Δ for a sample of HESE data events that were reconstructed as showers in Ref. [20].The left panel shows a comparison similar to the one in Fig. 3, where the reduced negative log-likelihood for a model that includes corrections due to both layer undulations and birefringence is given by  (2) .Turning off the layer undulation correction yields  (1) , and turning off all corrections  (0) .Since smaller  values corresponds to improved agreement against data, the improved description with  (2) is evidenced by the skew towards positive values.The right panel shows a similar comparison but changing the reference to  (3) , which is obtained with the additional two-cascade approximation as described in Section 6.1.In both panels, the vertical dashed line at Δ = 0 is a visual guide to highlight the positive skew.
As a final comparison, a check on data was performed using 109 HESE events that were reconstructed as showers in Ref. [20]. Figure 18 shows improvements in the reduced negative log-likelihood, , similar to that shown for the benchmark MC in Fig. 3.The left panel shows Δ distributions relative to  (2) , where  (2) is obtained using the best-fit from Eq. (5.1),  (1) using Eq.(4.3), and  (0) using Eq.(3.2), all with the fixes described in Section 3. As smaller  values correspond to a better description of data, the positive skews in the Δ distributions illustrate the improvement in describing data when using Eq.(5.1).The skew is larger for  (0) −  (2) than for  (1) −  (2) because Eq. (4.3) includes the birefringence correction.The right panel shows similar comparisons except now relative to  (3) , which is the best-fit obtained using Eq.(5.1) and the two-cascade approximation.We see that Δ shifts to more positive values, indicating that the two-cascade model yields a further improvement in the description of data.
While it may seem that increases in ice model complexity pose a challenge for the accurate reconstruction of particle showers in IceCube, in reality the problem can be naturally broken down into a few pieces, thus restoring some symmetry and reducing dimensionality.Birefringence is a global property that can be approximated with a few coordinate transformations, while ice layer undulations can be approximated by shifting the shower's physical depth to its corresponding P-depth.The shower profile itself can then be approximated using a two-cascade model.Implementing these features leads to improvement in the angular resolution by more than a factor of three.Hopefully, this work will aid future analyses and inform technical progress towards reaching the statistical limit for particle shower reconstruction in the South Pole ice.

Figure 2 .
Figure 2. Event view of an approximately 2 PeV particle shower detected by IceCube on December 4, 2012.Each string is shown as a thin line extending from top to bottom; small dots correspond to unhit DOMs.Modules that detected a photoelectron are indicated by the colored spheres.The size of each sphere corresponds to the total charge detected, and its color indicates the timing of hits, with earlier to later going from red to green.The IceCube-centered coordinate axes are indicated by the three colored arrows, with the blue arrow indicating the positive  direction.The length of each arrow is set to 100 m and gives a sense of scale.Note that noise cleaning has been applied for visualization purposes [1].

Figure 6 .
Figure 6.Visualization of the time-integrated photoelectron (P.E.) yield as a function of  and Θ s , sliced across the other dimensions at the values indicated in the title.The left (right) panel shows the Mie (updated) model expectations.Dashed white lines bound regions in the phase space where receiver locations fall within ±10 m of  s in depth, corresponding to a region of clear ice where the source is placed.The ringing in the left panel is due to the tensor product B-spline construction, which has difficulties interpolating along diagonals when support points are sparse.This is mitigated by increasing the number of support points in Θ s , the result of which is shown in the right panel.

Figure 7 .
Figure 7. Same as Fig. 6 except at a different slice,  = 90 • .The dashed white lines bound regions in the phase space where receiver locations fall within ±10 m of  s in depth, corresponding to a layer of clear ice where the source is placed.Here, large amplitude regions no longer lie along diagonals in (, Θ s ) space but instead align with the coordinate axes.As a result, interpolation artifacts do not appear in either the Mie model (left panel) or the updated model with denser support points (right panel).

Figure 8 .
Figure8.The time-integrated photoelectron (P.E.) yield as a function of Θ s for two receivers placed at a fixed distance  = 100 m from the source.The receivers are placed at two different directions, Θ r , of Θ r = 60 • (solid) and Θ r = 80 • (dashed), where Θ r is the angle between source-receiver vector and ẑ as shown in Fig.4.A discontinuity is visible in the Mie model (gray), which is smoothed out with the updates discussed in this work (red).The interpolation artifacts are also visible as small fluctuations in the dashed lines, and are reduced by the updates discussed in Section 3.1.

1Figure 9 .
Figure 9.The left panel shows quartiles of the distribution between reconstructed and true directions, Ψ, as a function of deposited EM-equivalent energy for the benchmark MC sample described in the text.Results obtained using the Mie model are shown in gray.Results obtained with a newer ice model that includes the fixes described in Section 3 are shown in red.Neither include corrections to approximate the additional ice effects to be discussed in Section 4 and Section 5.The solid red and dashed gray lines indicate their respective median angular resolutions.Lower (upper) cut offs for each color band show the respective 25 (75) percentile levels of the Ψ distribution.The right panels show distributions of Ψ in two different energy slices, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.

Figure 10 .
Figure 10.A visualization of the effect of birefringence on   for an isotropic light source centered at the origin of the detector.Solid lines are circles with radii of 60 m (blue) and 120 m (orange), centered at the origin.In the left panel, the dashed lines show the corresponding   across the -plane.Notice that   shrinks closer to the origin along the ice flow axis (black arrow), and shifts away from the origin perpendicular to the flow axis.The star (dot) indicates the location of the receiver in the left (right) panel of Fig. 11.In the right panel, a similar effect is observed along two planes perpendicular to the -plane, intersecting at the origin.The dashed (dotted) lines show   in the -plane where û = cos  x + sin  ŷ for  = 130 • (40 • ) denotes the unit vector parallel (perpendicular) to the ice flow axis shown in the left panel.

Figure 11 .
Figure11.Impact of birefringence on the arrival time PDF from an isotropic point-like light source, at a source depth of  s = 0 and the receiver coordinates indicated in the titles.For reference, the location of the receiver in the left (right) panel is indicated by the star (dot) in the left panel of Fig.10.The black, solid (red, dashed) line shows the time PDF with (without) birefringence simulated.The green, solid line shows the result of a two-parameter translation in  and , which when applied to the time PDF without birefringence gives the best agreement with the black line.

1Figure 12 .
Figure12.Same as Fig.9but now comparing the angular resolution using Eq.(4.3) (green, solid) to Eq. (3.2) (red, dashed), where both include the fixes described in Section 3. In the left panel, a substantial improvement is observed when birefringence corrections are included.Note the red intervals are identical to those in Fig.9.The right panels show distributions of Ψ in two different energy slices, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.

Figure 13 .
Figure 13.The left panel shows the elevation change for an ice isochron at a P-depth of 2,248 m, where P is the reference point indicated by the red dot.Black squares corresponds to the (, ) positions of IceCube strings, and the color map shows the elevation change relative to  P .The region outside the instrumented footprint is extrapolated [15].The dashed line indicates the -axis shown in the right panel, which lies perpendicular to the ice flow direction.The right panel shows depths of ice isochrons along û = cos  x + sin  ŷ, where  = 40 • , and the dashed line corresponds to the layer shown in the left panel.The visualized plane is perpendicular to the ice flow axis and intersects P at  = 0 m.Note that the relative change in each layer is stronger at the deeper regions of the instrumented volume.
I −  P (, ,  I ), where  I is defined such that  P (, ,  I ) = −300 m.Black squares indicate the location of the 86 strings and the color map shows the ice layer elevation change at different (, ) locations.The black arrow indicates the direction of ice flow, and the dashed line indicates the -axis shown in the right panel.The right panel of Fig. 13 shows the depths of ice isochrons in the plane normal to the flow direction and intersecting P at  = 0 m, where û = cos  x + sin  ŷ and  = 40 • .Note that the relative change in each layer increases towards at the deeper regions of the detector.

1Figure 15 .
Figure15.Same as Fig.12but now comparing the angular resolution using Eq.(5.1) (orange, solid) to Eq. (4.3) (green, dashed).In the left panel, an improvement is observed in the Ψ distribution quartiles when ice layer undulations are additionally included.Note the green intervals are identical to those in Fig.12.The right panels show distributions of Ψ in two different energy slices, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.

1
Figure16.Using the benchmark MC, the left panel shows the distribution of the reconstructed distance between two cascades for three different energy regimes.There is a trend towards slightly longer distances as the true EM-equivalent deposited energy,  dep , increases, which follows expectations from shower physics[36].The right panel shows distributions of the pulls on reconstructed energy, denoted as Êdep , for the same energy slices as in the left panel, illustrating that its resolution is on the order of a few percent.

1Figure 17 .
Figure 17.Same as Fig. 15 but now comparing the angular resolution using Eq.(5.1) with a two-cascade model (blue, solid) to a single cascade model (orange, dashed).For reference, results obtained with the Mie model without corrections are shown in dotted gray.The left panel shows that the most accurate Ψ quartiles are obtained with a two-cascade model.Note the orange intervals are identical to those in Fig. 15, and the gray to those in Fig. 9.The right panel shows distributions of Ψ in two different energy slices, between 10 TeV to 100 TeV (top) and 1 PeV to 10 PeV (bottom), with line colors and styles matching those of the left panel.