In-situ calibration of the single-photoelectron charge response of the IceCube photomultiplier tubes

We describe an improved in-situ calibration of the single-photoelectron charge distributions for each of the in-ice Hamamatsu Photonics R7081-02[MOD] photomultiplier tubes in the IceCube Neutrino Observatory. The characterization of the individual PMT charge distributions is important for PMT calibration, data and Monte Carlo simulation agreement, and understanding the effect of hardware differences within the detector. We discuss the single photoelectron identification procedure and how we extract the single-photoelectron charge distribution using a deconvolution of the multiple-photoelectron charge distribution.

: We describe an improved in-situ calibration of the single-photoelectron charge distributions for each of the in-ice Hamamatsu Photonics R7081-02 [MOD] photomultiplier tubes in the IceCube Neutrino Observatory. The characterization of the individual PMT charge distributions is important for PMT calibration, data and Monte Carlo simulation agreement, and understanding the effect of hardware differences within the detector. We discuss the single photoelectron identification procedure and how we extract the single-photoelectron charge distribution using a deconvolution of the multiple-photoelectron charge distribution.

Introduction
The IceCube Neutrino Observatory [1,2] is a cubic-kilometer-sized array of 5,160 photomultiplier tubes (PMTs) buried in the Antarctic ice sheet, designed to observe high-energy neutrinos interacting with the ice [3]. In 2011, the IceCube Collaboration completed the installation of 86 vertical strings of PMT modules, eight of which were arranged in a denser configuration known as the DeepCore sub-array [4]. Each string in IceCube contains 60 digital optical modules (DOMs), which contain a single PMT each, as well as all required electronics [5]. The primary 78 strings (excluding DeepCore) are spaced 125 m apart in a hexagonal grid, with the DOMs extending from 1450 m to 2450 m below the surface of the ice sheet. The additional DeepCore strings (79-86) are positioned between the centermost strings in the detector, reducing the horizontal DOM-to-DOM distance in this region to between 42 m and 72 m. The lower 50 DOMs on these strings are located in the deepest 350 m of the detector surrounded by the cleanest ice [6], while the upper ten provide a cosmic ray veto extending down from 1900 m to 2000 m below the surface. Above the in-ice detectors, there exists a surface array, IceTop [7], consisting of 81 stations located just above the in-ice IceCube strings. The PMTs located in IceTop DOMs operate at a lower gain and the data from these PMTs was not included in the current analysis; however, the IceTop PMTs are calibrated The version of AC coupling, old toroids (dark green) and new toroids (light green). DOMs that have been removed from service (OTS) are shown in white. These are primarily DOMs that resulted in failures during deployment and freeze-in, and failures during subsequent operation, as described in Ref. [9].
to single photoelectron charge distribution in a similar way as the in-ice PMTs (see Sec. 5.1 in Ref. [7]).
Each DOM consists of a 0.5"-thick spherical glass pressure vessel that houses a single downfacing 10" PMT from Hamamatsu Photonics. The PMT is coupled to the glass housing with optical gel and is surrounded by a wire mesh to reduce the effect of the Earth's ambient magnetic field. The glass housing is transparent to wavelengths of 350 nm and above [8].
Of the 5,160 DOMs, 4,762 house a R7081-02 Hamamatsu Photonics PMT, sensitive to wavelengths ranging from 300 nm to 650 nm, with a peak quantum efficiency of 25% near 390 nm. These are classified as Standard Quantum Efficiency (Standard QE) DOMs. The remaining 398 DOMs are equipped with the Hamamatsu R7081-02MOD PMTs, which, having a peak quantum efficiency of 34% near 390 nm (36% higher efficiency than the Standard QE DOMs), are classified as High Quantum Efficiency (HQE) DOMs [4]. These DOMs are primarily located in DeepCore and on strings 36 and 43, as shown in the left side of Fig. 1.
The R7081-02 and R7081-02MOD PMTs have 10 dynode stages and are operated with a nominal gain of 10 7 and achieved with high voltages ranging from approximately 1215 ± 83 V and 1309 ± 72 V, respectively. A typical amplified single photoelectron generates a 5.2 ± 0.3 mV peak voltage after digitization with a full width at half maximum of 13 ± 1 ns. The PMTs are operated with the anodes at high voltage, so the signal is AC coupled to the amplifiers (front-end amplifiers). There are two versions of AC coupling in the detectors, referred to as the new and old toroids, both of which use custom-designed wideband bifilar wound 1:1 toroidal transformers1. The locations of DOMs with the different versions of AC-coupling are shown on the right side of Fig. 1. The DOMs with the old toroids were designed with an impedance of 43 Ω, while the new toroids are 50 Ω [9]. All HQE DOMs are instrumented with the new toroids.
IceCube relies on two observables per DOM to reconstruct events: the total number of detected 1The toroidal transformer effectively acts as a high-pass filter with good signal fidelity at high frequencies and offers a higher level of reliability than capacitive coupling. Conventional AC-coupling high-voltage ceramic capacitors can also produce undesirable noise from leakage currents and are impractical given the signal droop and undershoot requirements [8].
photons and their timing distribution. Both the timing and the number of photons are extracted from the digitized waveforms. This is accomplished by deconvolving the digitized waveforms [10] into a series of scaled single photoelectron pulses (so-called pulse series), and the integral of the individual pulses divided by the load resistance defines the observed charge. It will often be expressed in units of PE, or photoelectrons, which further divides the measured charge by the charge of a single electron times the nominal gain.
When one or more photoelectrons produce a voltage at the anode sufficient to trigger the onboard discriminator, the signal acquisition process is triggered. The discriminator threshold is set to approximately 1.2 mV, or equivalently to ∼0.23 PE, via a digital-to-analog converter (DAC). The signal is presented to four parallel channels for digitization. Three channels pass through a 75 ns delay loop in order to capture the waveform leading up to the rising edge of the triggering pulse, and are then subject to different levels of amplification prior to being digitized at 300 million samples per second (MSPS) for 128 samples using a 10-bit Analog Transient Waveform Digitizer (ATWD). The high-gain channel has a nominal amplification of 16 and is most suitable for single photon detection. Two ATWD chips are present on the DOM Mainboard (MB) and alternate digitization between waveforms to remove dead time associated with the readout. The signal to the fourth parallel channel is first shaped and amplified, then fed into a 10-bit fast analog-to-digital converter (fADC) operating at a sampling rate of 40 MSPS. Further detail regarding the description of the DOM electronics can be found in Refs. [5,11].
This article discusses a method for determining the in-situ single-photoelectron charge distributions of individual PMTs, which can be used to improve calibration and the overall detector description in Monte Carlo (MC) simulation. The SPE charge distribution refers to the charge probability density function of an individual PMT generated by the amplification of a pure sample of single photoelectrons. The measured shape of the SPE charge distributions is shown to be useful for examining hardware differences and assessing long term stability of the detector. This was made possible with the development of two pieces of software: 1. A specially-designed unbiased pulse selection developed to reduce the multiple photoelectron (MPE) contamination while accounting for other physical phenomena (e.g. late pulses, afterpulses, pre-pulses, and baseline shifts) and software-related effects (e.g. pulse splitting). This is further described in Sec. 2.1.

2.
A fitting procedure developed to separate the remaining MPE contamination from the SPE charge distribution by deconvolving the measured charged distribution. This is further described in Sec. 2.3.
By using in-situ data to determine the SPE charge distributions, we accurately represent the individual PMT response as a function of time, environmental conditions, software version and hardware differences, and realistic photocathode illumination conditions. This is beneficial since it also allows us to inspect the stability and long-term behavior of the individual DOMs, verify previous calibration, and correlate features with specific DOM hardware.

Single-photoelectron charge distributions
Ideally, a single photon produces a single photoelectron, which is then amplified by a known amount, and the measured charge corresponds to 1 PE. However, there are many physical processes that create structure in the measured charge distributions. For example: • Statistical fluctuation due to cascade multiplication [12]. At every stage of dynode amplification, the number of emitted electrons that make it to the next dynode is randomly distributed. This in turn causes a smearing in the measured charge after the gain stage of the PMT.
• Photoelectron trajectory. Some electrons may deviate from the favorable trajectory, reducing the number of secondaries produced at a dynode or the efficiency to collect them on the following dynode. This can occur at any stage, but it has the largest effect on the multiplication at the first dynode [13]. The trajectory of a photoelectron striking the first dynode will depend on many things, including where on the photocathode it was emitted, the uniformity of the electric field, the size and shape of the dynodes [12], and the ambient magnetic field [14,15].
• Late or delayed pulses. A photoelectron can elastically or inelastically backscatter off the first dynode. The scattered electron can then be re-accelerated to the dynode, creating a second pulse. The difference in time between the initial pulse and the re-accelerated pulse in the R7081-02 PMT was previously measured to be up to 70 ns [8,16]. Elastically backscattered photoelectrons will carry the full energy and are thus expected to produce similar charge to a non-backscattered photoelectron, albeit with a time offset. The mean measured charge of an ineleastic backscattered photoelectron, by contrast, is expected to be smaller than a nominal photoelectron [17].
• Afterpulses. When photoelectrons or the secondary electrons produced during the electron cascade gain sufficient energy to ionize residual gas in the PMT, the resulting positively charged ionized gas will be accelerated in the electric field towards the photocathode. Upon impact with the photocathode, multiple electrons can be released from the photocathode, creating what is called an afterpulse. For the R7081-02 PMTs used in IceCube, the timescale for afterpulses was measured to occur from 0.3 to 11 µs after the initial pulse, with the first prominent afterpulse peak occurring at approximately 600 ns [8]. The spread in the afterpulse time depends on the position of photocathode, the charge-to-mass ratio of the ion produced, and the electric potential distribution [18], whereas the size of the afterpulse is related to the momentum and species of the ionized gas and composition of the photocathode [19].
• Pre-pulses. If an incident photon passes through the photocathode without interaction and strikes the first dynode, it can eject an electron that is only amplified by the subsequent stages, resulting in a lower measured charge (lower by a factor of approximately 20). For the IceCube PMTs, the prepulses have been found to arrive approximately 30 ns before the signal from photoelectrons from the photocathode [8].
• MPE contamination. When multiple photoelectrons arrive at the first dynodes within few nanoseconds of each other, they can be reconstructed by the software as a single MPE pulse.
• Dark noise. Photoelectron emission, not initiated from an external event, can be attributed to thermionic emission from the low work function photocathode and the dynodes, Cherenkov radiations initiated from radioactive decay within the DOM, and field emission from the electrodes. Dark noise originating from thermionic emission from the dynodes is shown in Ref. [20] to populate the low-charge region.
• Electronic noise. This refers to the combined fluctuations caused by noise generated from the analog-frontend and the analog-to-digital converters (ATWDs and fADC). When integrated over a time window the resulting charge is generally small and centered around zero, thus only leading to a small broadening in the low charge region. The standard deviation of the electronic noise was found to be approximately ±0.11 mV.
Beyond the physical phenomena above that modify the measured charge distribution, there is also a lower limit on the smallest charge that can be extracted. For IceCube, the discriminator only triggers for peak voltages above the threshold and subsequent pulses in the readout window are subject to a threshold defined in the software. This software threshold was set conservatively to avoid extracting pulses that originated from electronic noise. It can be modified to gain access to lower charge pulses and will be discussed in Sec. 2.2.
The standard SPE charge distribution used for all DOMs in IceCube, known as the TA0003 distribution [8], models the above effects as the sum of an exponential plus a Gaussian. The TA0003 distribution is the average SPE charge distribution extracted from a lab measurement of 118 Hamamatsu R7081-02 PMTs. The measurement was performed in a -32 • C freezer using a pulsed UV LED centered along the axis of the PMT, directly in front of the photocathode.
In 2013, IceCube has made several lab measurements of the SPE charge distribution of R7081-02 PMTs using single photons generated from synchronized short duration laser pulses. The coincident charge distribution generated by the laser pulses was found to include a steeply falling low-charge component in the region below the discriminator threshold. To account for this, a new functional form including a second exponential was introduced. This form of the normalized charge probability distribution f (q) SPE = Exp 1 + Exp 2 + Gaussian, is referred to as the SPE charge template in this article. Explicitly, it is: where q represents the measured charge; w 1 and w 2 are the exponential decay widths; and µ, σ are the Gaussian mean and width, respectively. The coefficients P e1 , P e2 , and 1 − P e1 − P e2 correspond to the probability of a photoelectron contributing to each component of the SPE template. The Erfc function used to normalize the Gaussian represents the complementary error function. Eq. 1.1 is the assumed functional shape of the SPE charge distributions, and the components of Eq. 1.1 are determined in this article for all in-ice DOMs. IceCube has chosen to define 1 PE as the location of the Gaussian mean (µ) and calibrates the gain of the individual PMTs prior to the start of each season to meet this definition. Any overall bias in the total observed charge can be absorbed into an efficiency term, such as the quantum efficiency. This is valid since the linearity between the instantaneous total charge collected and the number of incident photons is satisfied up to ∼2 V [9], or approximately 560 PE. That is, the average charge collected from N photons is N times the average charge of the SPE charge distribution, and the average charge of the SPE charge distribution is always a set fraction of the Gaussian mean.

IceCube datasets and software definitions
The amount of observed light depends on the local properties of the ice [6]. Short term climate variations from volcanoes and longer-term variations from atmospheric dust affect the optical properties of the ice, producing nearly horizontal layers. This layered structure affects how much light the DOMs observe, and, with it, the trigger rate. The largest contribution to the IceCube trigger rate comes from downward-going muons produced in cosmic ray-induced showers [21]. Cosmic ray muons stopping in the detector cause the individual trigger rates to decrease at lower depths. If a DOM and its nearest or next-to-nearest neighbor observe a discriminator threshold crossing within a set time window, a Hard Local Coincidence (HLC) is initiated, and the corresponding waveforms are sampled and read out on the three ATWD channels. Thermionic emission induced dark noise can be present in the readout, however it is suppressed at lower temperatures and is unlikely to trigger an HLC event.
After waveform digitization, there is a correction applied to remove measured baseline offsets. Distortions to the waveform, such as from droop and undershoot [8] introduced by the toroidal transformer AC coupling are compensated for in software during waveform calibration by adding the expected reaction voltage of the distortion to the calibrated waveform. If the undershoot voltage drops below 0 ADC counts, the ADC values are zeroed and then compensated for once the waveform is above the minimum ADC input. For each version of the AC coupling, scaled single photoelectron pulse shapes are then fit to the digitized waveforms using software referred to as "WaveDeform" (waveform unfolding process), which determines the individual pulse times and charges and populates a pulse series.
The pulse series used in this analysis come from two datasets: 1. The MinBias dataset. This dataset preserves the full waveform readout of randomly-selected HLC events, collecting on average 1:1000 events. The largest contribution to this dataset comes from downward-going muons produced in cosmic-ray-induced showers. The average event for this sample is approximately 26 PE bright and distributed over an average of 16 triggered DOMs. The full waveform of these events allows us to extract the raw information about the individual pulses. This dataset will be used to measure the individual PMT charge distributions.
2. The BeaconLaunch dataset. This dataset is populated with digitized waveforms that are initiated by the electronics (forced-triggered) of a channel that has not gone above the threshold. The forced triggered waveforms are typically used to monitor the individual DOM baselines and thus includes the full ATWD waveform readout. Since this dataset is forced-triggered, the majority of these waveforms represent electronic noise with minimal contamination from random accidental coincidence SPEs. This dataset will be used to examine the noise contribution to the charge distributions.
When using this dataset, the weight of every pulse is multiplied by a factor of 28.4 to account for the livetime difference between the MinBias dataset and the BeaconLaunch dataset. Weight, in this context, refers to the number of photons in the MinBias dataset proportional to one statistical photon in the BeaconLaunch dataset for which both datasets have the same equivalent livetime.
This analysis uses the full MinBias and BeaconLaunch datasets from IceCube seasons 2011 to 2016 [22], subsequently referred to as IC86.2011 to IC86.2016. Seasons in IceCube typically start in May of the labeled year and end approximately one year later. Calibration is performed before the start of each season.

Extracting the SPE charge templates 2.1 Single photoelectron pulse selection
The pulse selection is the method used to extract candidate, unbiased, single photoelectron pulses from high-gain ATWD channel while minimizing the MPE contamination. The pulse selection was designed such that it avoids collecting afterpulses, does not include late pulses from the trigger, accounts for the discriminator threshold, reduces the effect of signal droop and undershoot, and gives sufficient statistics to perform a season-to-season measurement. An illustrative diagram of the pulse selection is shown in the left side of Fig. 2, while a description of the procedure is detailed below.
We restrict the pulse selection to only extract information from waveforms in which the trigger pulse does not exceed 10 mV (∼2 PE) and no subsequent part of the waveform exceeds 20 mV (∼4 PE). This reduces the effect of the baseline undershoot due to the AC coupling or other artifacts from large pulses.
In order to trigger a DOM, the input to the front-end amplifiers must exceed the discriminator threshold. To avoid the selection bias of the discriminator trigger (i.e. only selecting pulses greater than the discriminator threshold), we ignore the trigger pulse as well as the entire first 100 ns of the time window. Ignoring the first 100 ns removes late pulses that could be attributed to the triggering pulse, which occurs approximately 4% of the time [8]. To ensure we are not accepting afterpulses into the selection, we also enforce the constraint that the pulse of interest (POI) is within the first 375 ns of the ATWD time window. This also allows us to examine the waveform up to 50 ns after the POI. In the vicinity of the POI, we ensure that WaveDeform did not reconstruct any pulses up to 50 ns prior to the POI, or 100 to 150 ns after the POI (the light gray region of Fig. 2 (left)). This latter constraint is to reduce the probability of accidentally splitting a late pulse in the summation window.
If a pulse is reconstructed between 100 and 375 ns after the start of the waveform and the voltage criteria are met, it is accepted as a candidate photoelectron and several checks are performed on the waveform prior to and after the pulse. The first check is to ensure that the waveform is near the baseline just before the rising edge of the POI. This is accomplished by ensuring that the waveform does not exceed 1 mV, 50 to 20 ns prior to the POI, and eliminates cases where the POI is a late pulse. We also ensure the waveform returns to the baseline by checking that no ADC measurement exceeds 1 mV, 100 to 150 ns after the POI. These constraints are illustrated as the horizontal red dotted lines and black arrows in the left side of Fig. 2.
If all the above criteria are met, we sum the reconstructed charges from the POI time, given by WaveDeform, to +100 ns (the dark gray area in Fig. 2 (left)). This ensures that any nearby pulses are either fully separated or fully added. This is important since WaveDeform may occasionally split an SPE pulse into multiple smaller pulses, therefore it is always critical to perform a summation of the charge within a time window. The 100 ns summation also implies that the pulse selection will occasionally accept MPE events. We chose 100 ns window for the summation to ensure that we collect the charge of the late pulse (recall that late pulses were measured up to 70 ns after the main pulse), should it be there, while minimizing the MPE contamination. We estimate that there is on average a 6.5% probability of the summation time window includes two or more photons.

Characterizing the low-charge region
This analysis aims to describe the full SPE charge distribution for each DOM. This is required by the IceCube simulation. However, we cannot extract charge to arbitrary low PE before electronic noise starts dominating. The aim of this section is to describe how we extract information in the low-charge region (below 0.25 PE) to guide the full fit. Fig. 2 (right) shows the charge distributions of the selected pulses that pass the single photoelectron pulse selection for string 1, optical module 1, DOM(1,1). In the low-charge region, we see a second threshold at approximately 0.13 PE, i.e. the charge distribution terminates. This threshold arises from a termination condition in WaveDeform, in which the pulses that are smaller than predefined criteria are rejected. The threshold was set to avoid electronic noise being interpreted as PMT pulses and contaminating the low-charge region.
The steeply falling component of the region from 0.13 PE to 0.25 PE is in agreement with the laser measurements mentioned in Sec. 1.1 and emphasizes the importance of collecting data below the discriminator threshold. This section will assess the noise contribution to this region and examine the effect on the charge distribution and noise contribution by lowering the WaveDeform threshold.   Fig. 3 (left) shows the charge distributions for the MinBias (black) and the BeaconLaunch (red) datasets using the default settings of WaveDeform (standard WaveDeform). As mentioned in Sec. 1.2, occasionally a photoelectron will be coincident with the forced BeaconLaunch time window. These charges populate a SPE charge distribution. Subtracting the shape of the MinBias charge distribution from the BeaconLaunch dataset yields an estimate of the amount of electronic noise contamination (blue). The bin in the MinBias data with the lowest signal-to-noise ratio (SNR) above 0.1 PE was found to have a SNR of 744.7. The SNR for the full distribution was found to be 1.98×10 5 . Fig. 3 (right) shows the same data after lowering the WaveDeform threshold (modified WaveDeform), and is found to have SNR of 57.9 in the bin with the largest contamination and the total SNR was found to be 0.69×10 5 .
The modified WaveDeform datasets show a minimal increase in the contribution of noise to the low-charge region. From this, however, we are able to extract charge information down to approximately 0.10 PE and improve the overall description of the charge distribution below the discriminator. This will help constrain the values of the steeply falling exponential, defined with Exp 1 .

Fitting procedure
We would now like to fit the charge distribution to extract the SPE charge templates (the components of Eq. 1.1) for all DOMs.
Contamination from two-photon events is suppressed by the pulse selection, but can not be entirely avoided. To minimize potential biases by the charge entries resulting from two photons, the one and two photon contribution to the charge distributions is fitted at the same time, using something we call a convolutional fitter. It assumes that the charge distribution resulting from two photons is the SPE charge distribution convolved with itself [23]. In each step of the minimizer the convolution is updated given the current set of SPE parameters to be evaluated and the relative one and two photon contributions is determined.
We do not account for the three-photon contribution, which is justified by the lack of statistics in the 3 PE region as well as the significant rate difference between the 1 PE and 2 PE region, as shown in Fig. 2 (right).
Pulses that fall below the WaveDeform threshold and are not reconstructed contribute to an inefficiency in the individual DOMs. That is, the shape below the WaveDeform software threshold does not have a significant impact, but the relative area of the SPE charge template below compared to above this threshold changes the efficiency of the DOM. This analysis assumes the same shape of the steeply falling exponential component (Exp 1 ) for all DOMs in the detector to avoid large fluctuations in the DOM-to-DOM efficiencies. The modified WaveDeform data will strictly be used to determine the Exp 1 component. Specifically, using the aggregate of the entire ensemble of DOMs with the modified WaveDeform dataset, we background-subtract the BeaconLaunch distribution from the MinBias data, fit the resulting distribution to determine the components of Eq. 1.1, and use only the measured shape and normalization of Exp 1 in all subsequent standard WaveDeform fits.
As described in Sec. 1.1, the Gaussian mean (µ) is used to determine the gain setting for each PMT. Therefore, it is particularly important that the fit quality in the peak region accurately describes the data. While fitting to the full charge distribution improves the overall fit agreement, the mismatch between the chosen functional form (Eq. 1.1) and a true SPE charge distribution can cause the Gaussian component to pull away from its ideal location. To compensate for this, the fitting algorithm prioritizes fitting to the data around the Gaussian mean. This is accomplished by first fitting to the full distribution to get an estimate of the Gaussian mean location. Then, the data in the region ±0.15 PE around the original estimated Gaussian mean is weighted to have a higher impact on the fit, and the distribution is re-fitted.
Upon fitting the MinBias dataset with the predetermined values for Exp 1 , the residual of each fit is calculated by measuring the percentage difference between the fit and the data. The average residual is then used as a global scaling factor for all SPE charge templates to account for the difference between the chosen model (Eq. 1.1) and the actual data.

SPE charge template fit results
We now present the results of the fits and describe the correlations of the fit parameters with hardware differences, and time variations in the next section. Using the background-subtracted modified WaveDeform dataset, the Exp 1 component was determined by fitting the aggregate distribution from 0.1 PE to 3.5 PE. The result of the fit yielded P e1 = 0.186 ± 0.041 and w 1 = 0.027 ± 0.002 PE. This shape of Exp 1 is now used to describe the low-PE charge region for all subsequent standard WaveDeform fits.  Using the MinBias dataset with the measured values of Exp 1 , the SPE charge templates are extracted for every DOM, separately for each IceCube season from IC86.2011 to IC86.2016. The fit range for Exp 2 and the Gaussian components is selected to be between 0.15 PE and 3.5 PE. An average fit was also performed on the cumulative charge distribution, in which all the data for a given DOM was summed (labeled as "AVG").
All the DOMs with "failed fits" are not included in this analysis. A DOM is classified as having a failed fit if it does not pass one of the validity checks on the data requirements (e.g. the number of valid pulses) or goodness of fit. Over the seasons considered, between 107 and 111 DOMs have been excluded from service (the white elements shown in Fig. 1) and represent the majority of the failed fits. The remaining 6 DOMs that failed the AVG fits are known to have various issues. In the IceCube MC simulation chain, these DOMs are assigned the average SPE charge template.
We can divide the DOMs into subsets with the following hardware differences: the HQE DOMs with the new toroids, the Standard QE DOMs with the new toroids, and the Standard QE DOMs with the old toroids. The mean value and standard error of the IC86.AVG fit parameters, excluding Exp 1 , for the subset of hardware differences are listed in Table 1. The residual, averaged over all DOMs, from 0 to 1 PE is shown in Fig. 4 An example fit is shown in Fig. 5 for the cumulative MinBias charge distribution for DOM (1,1). The collected charge distribution is shown in the black histogram, while the fit to the data is shown as the red line. The extracted SPE charge template from the fit is shown in blue. Both the fit and extracted SPE charge template have been scaled by the average residual shown in Fig. 4.

Correlations between fit parameters and DOM hardware differences
It is evident from the data in Table 1 that the average shape of the SPE charge templates depends on the DOM hardware. These differences can also be seen in the measured peak-to-valley ratios and average charge of the SPE charge template (see Fig. 6). When we examine the subset of DOMs instrumented with the new toroids, the average HQE DOM were found to have a 20.2 ± 0.6% larger P e2 component and 4.7 ± 0.4% smaller Gaussian width. Consequently, the average HQE peak-tovalley ratio is measured to be 2.255 ± 0.012, or 13.67 ± 0.01% lower than the average for Standard QE DOMs. Also, the average charge of the average HQE DOM was found to be 3.17 ± 0.01% lower than that of the Standard QE DOMs. The average charge is calculated by integrating over the full SPE charge template including the residual correction. The values shown in Fig. 6 (right) are found to be below 1 PE due to the low-PE contribution from Exp 1 and Exp 2 , whose physical description can be found in Sec. 1.1.
IceCube compensates for the change in the mean measured charge in simulation, by increasing the HQE DOM efficiency by the equivalent amount. This ensures that the total amount of charge collected by the HQE DOMs remains the same prior to, and after, inserting the SPE charge templates into simulation.
Similarly, using only the subset of Standard QE DOMs, the SPE charge templates were found to have measurably different shapes when comparing the different methods of AC coupling. The average Gaussian width for the DOMs instrumented with the old toroids were found to be 8.1 ± 0.1% narrower, albeit with a very similar probability of the photoelectron populating the Gaussian component. With these differences, we find a peak-to-valley ratio of 2.612 ± 0.007 for the new toroid DOMs and 2.987 ± 0.012 for the old toroid DOMs. The average Gaussian mean of the fit for the DOMs with the old toroids was also found to be 1.6 ± 0.1% lower than those with the new toroids. This corresponds proportionally to a change in the expected gain. The average charge, however, between these two hardware configurations is very similar, 0.012 ± 0.001%.
Although the DOMs instrumented with the old toroids were deployed into the ice earlier than those with the new toroids, the differences above remain when examining individual deployment years; therefore, the shape differences are not attributed to the change in the DOM behavior over time. However, the DOMs with the old toroids were the first PMTs to be manufactured by Hamamatsu. A gradual change of the fit parameters was observed when ordering the PMTs according to their PMT serial number (i.e. their manufacturing order). Fig. 7 shows the change in the measured peak-to-valley ratio as a function of PMT serial number for the standard QE DOMs (blue) and HQE PMTs (red). Here, each data point represents a single PMT and the blue (red) data points indicate a PMT instrumented with the new (old) toroid. This is compelling evidence that the observed differences between the new and old toroids is due to a change in the PMT production procedure rather than version of AC coupling. Fig. 8 illustrates the average shape differences in the extracted SPE charge templates between the HQE DOM with the new toroids (solid white line), Standard QE with the new toroids (dotted white line), Standard QE with the old toroids (dashed white line), compared to the spread in the measured SPE charge templates for all DOMs in the detector (dark blue contours). The figure also shows how the previous default SPE charge distribution, the TA0003 distribution, compares to this recent measurement. All curves in this figure have been normalized such that the area above 0.25 PE is the same. The observable shape differences from the TA0003 are attributed to a better understanding of the low-charge region, the difference in functional form (described in Section 1.1), and the fact that the SPE charge templates were generated using a realistic photocathode illumination.

Fitting parameters variation over time
The SPE charge templates were extracted for each IceCube season independently so to investigate the time dependence of the fit parameters. For every DOM in the detector, the change over time of each fit parameter (excluding Exp 1 ) was calculated. Fig. 9 shows the change in a given fit parameter, relative to the mean value, per year. The measured distributions were found to be consistent with randomly scrambling the order of each seasonal measurement. The largest deviation was found to be average of each fit parameters are found to deviate less than 0.1%, which is in agreement with the stability checks performed in Ref. [9]. This observation holds for the individual subsets of DOMs with different hardware configurations as well.

Quantifying observable changes when modifying the PMT charge distributions
Changing the assumed gain response in simulation has different implications depending on the typical illumination level present in different analyses. These differences are outlined in the following discussion.
The PMT response is described by a combination of a "bare" efficiency, η 0 , and a normalized charge response function, f (q). The bare efficiency represents the fraction of arriving photons that result in any nonzero charge response, including those below the discriminator threshold. The normalization condition is: Generally, f (q) and η 0 have to be adjusted together to maintain agreement with a quantity known from lab or in-ice measurements, such as the predicted number of pulses above threshold for a dim source.
Dim source measurements Where light levels are low enough, the low occupancy ensures that sub-discriminator pulses do not contribute to any observed charge as they do not satisfy the trigger threshold. Given some independent way of knowing the number of arriving photons, a lab or in-ice measurement determines the trigger fraction above threshold η 0.25 and/or the average charge over threshold Q 0.25 , either of which can be used to constrain the model as follows: Here, the discriminator threshold is assumed to be 0.25 times the peak position q pk . It is also useful to multiply observed charges by q pk , since we set each PMT gain by such a reference, and then a measurement constraint would be stated in terms of Q 0.25 /q pk .

Semi-bright source measurements
For semi-bright sources, pulses that arrive after the readout time window is opened are not subject to the discriminator threshold. WaveDeform introduces a software termination condition at ∼0.13 PE (described at the end of Section 2.1). The average charge of an individual pulse that arrives within the time window is: Bright source measurements For light levels that are large, the trigger is satisfied regardless of the response to individual photons, and the total charge per arriving photon therefore includes contributions below both the discriminator and the WaveDeform thresholds: As such, the total charge is directly proportional to the average charge of the SPE charge template.

Model comparison
A natural question to ask is whether or not a change in f (q) would cause observable changes in the bright-to-dim ratios. That is, when we change the SPE charge distribution in simulation, should we expect the charge collected by bright events compared to dim events to change? When the charge distribution model is changed in a way that preserves agreement with the measured η 0.25 or Q 0.25 /q pk , i.e. η 0 is adjusted properly for changes in f (q), the physical effect can be summarized by the change in the bright-to-dim ratios Q 0 /Q 0.25 , and Q 0 /Q 0.10 . Conveniently, these ratios depend only on the shape of f (q). Table 2 compares these ratios in terms of the TA0003 charge distribution and the SPE charge templates described here. It is shown that there are sub-percent level differences in the physically-observable bright-to-dim ratios. The largest difference in the shape between the SPE charge templates and the TA0003 distribution is in the low-charge region, particularly below ∼0.10 PE. Charge from this region can only inflate bright events. That is, these pulses are small to trigger the discriminator or be reconstructed by WaveDeform, however they can reside on top of other pulses, inflating them. Since these pulses by definition contain little charge, they do not tend to inflate the measured charge by a noticeable amount, as shown by the Q 0 /Q 0.25 measurements in Table 2.

Model
Detector

SPE charge templates for calibration
The gain setting on each PMT is calibrated prior to the beginning of each season such that the Gaussian mean of the charge distribution corresponds to a gain of 10 7 , or equivalently 1 PE. This gain calibration method, run directly on the DOMs, uses waveform integration for charge determination instead of WaveDeform unfolding, resulting in a small systematic shift in gain. This systematic shift was determined for every PMT. The mean shift obtained over all DOM was found to be 1.47 ± 0.04% with a standard deviation of 2.62%, corresponding to an overestimation of the measured charge in the detector. The correction to the systematic shift in the measured charge can be implemented retroactively by dividing the reported charge from WaveDeform by the corresponding offset for a given DOM. Alternatively, we can account for this by simply inserting SPE charge templates, measured in this analysis, into simulation such that the corresponding systematic shift is also modeled in simulation. This will be performed in the following subsection.

SPE charge templates in simulation
To model the IceCube instrument, we must implement the PMT response in simulation. The IceCube MC simulation chain assigns a charge to every photoelectron generated at the surface of the photocathode. The charge is determined by sampling from a normalized charge distribution probability density function (PDF). A comparison to data between describing the charge distribution PDF using the SPE charge templates and the TA0003 distribution follows.
Two simulation sets consisting of the same events were processed through the IceCube Monte Carlo simulation chain to the final analysis level of an update to the IC86.2011 sterile neutrino analysis [24]. Here, the events that pass the cuts are >99.9% upward-going (a trajectory oriented upwards relative to the horizon) secondary muons produced by charged current muon neutrino/antineutrino interactions. The muon reconstructed energy range of this event selection is between approximately 500 GeV and 10 TeV. Fig. 10 (left) shows the distribution of the total measured charge during each event per DOM (data points). The simulation set using the TA0003 charge distribution is shown in orange, and that using the SPE charge templates is shown in blue. The data is shown for the full IC86.2012 season but is statistically equivalent to any of the other seasons. Fig. 10 (right) shows the distribution of the total measured charge of an event divided by the number of channels (NChan), or DOMs, that participated in the event. Both plots in Fig. 10 have been normalized such that the area under the histograms is the same.
The SPE charge templates clearly improve the overall MC description of these two variables. This update may be useful for analyses that rely on low-occupancy events (low-energy or dim events) in which average charge per channels is below 1.5 PE, and will be investigated further within IceCube.

Conclusion
This article outlines the procedure used to extract the SPE charge templates for all in-ice DOMs in the IceCube detector using in-situ data from IC86.2011 to IC86.2016. The SPE charge templates represent an update to the modeled single photoelectron charge distribution previous used by IceCube, the TA0003 distribution. The result of this measurement was shown to be useful for improving the overall data/MC agreement as well as calibration of the individual PMTs. It also prompted a comparison between the shape of the SPE charge templates for a variety of hardware configurations and time dependent correlations.
The subset of HQE DOMs were found to have a smaller peak-to-valley ratio relative to the Standard QE DOMs, as well as an overall 3.17 ± 0.01% lower average charge. It was also found that the DOMs instrumented with the old toroids used for AC coupling (the first PMTs to be manufactured by Hamamatsu) had narrower Gaussian component corresponding in an average increased peak-tovalley ratio of 14.36 ± 0.01%. By assuming a relationship between the production time and serial number, this difference was shown to be likely due to a change in the manufacturing over time rather than the actual AC coupling method. No significant time dependence in any of the fitted parameters associated with the SPE charge templates over the investigated seasons was observed. A reassessment of the PMT gain settings found a systematic bias of 1.47 ± 0.04% with a standard deviation of 2.62%.
The SPE charge templates were inserted into the MC simulation and the results were compared to the default TA0003 distribution. A significant improvement in the description of the variables total charge per DOM and total charge over the number of channels was shown. Analyses which rely on low-light occupancy measurements, will benefit from this update. As shown in the bright-to-dim ratios, the average charge for various light levels will not be affected by this update.