Constraints on dark photons and axionlike particles from the SuperCDMS Soudan experiment

We present an analysis of electron recoils in cryogenic germanium detectors operated during the SuperCDMS Soudan experiment. The data are used to set new constraints on the axioelectric coupling of axionlike particles and the kinetic mixing parameter of dark photons, assuming the respective species constitutes all of the galactic dark matter. This study covers the mass range from 40 eV =c 2 to 500 keV =c 2 for both candidates, excluding previously untested parameter space for masses below ∼ 1 keV =c 2 . For the kinetic mixing of dark photons, values below 10 − 15 are reached for particle masses around 100 eV =c 2 ; for the axioelectric coupling of axionlike particles, values below 10 − 12 are reached for particles with masses in the range of a few-hundred eV =c 2 .


I. INTRODUCTION
Many astrophysical and cosmological observations support the existence of dark matter, which constitutes more than 80% of the matter in the universe [1,2]. While the current Standard Model of particle physics (SM) does not describe dark matter, a number of theoretical extensions predict new particles that are viable dark matter candidates. For the past three decades weakly interacting massive particles (WIMPs) have been the primary candidate of interest [3]. Recently however, low-mass candidates (OðkeV=c 2 Þ or below) such as axions or axionlike particles (ALPs) and dark photons have gained traction [3]; these particles occur naturally in many proposed models, and have been invoked to explain various experimental and observational anomalies [4][5][6]. Both axions/ALPs and dark photons may be produced in the early universe and thus may constitute a significant fraction (if not all) of the dark matter in the universe [7][8][9][10]. In theories that include ALPs or dark photons, a coupling between the new particle and a SM particle may lead to a process by which the new particle is absorbed by an atom and an electron is ejected, carrying away the excess energy [11]. This process is analogous to the photoelectric effect and is henceforth referred to as dark absorption. As galactic dark matter is nonrelativistic, the observable signature from relic ALPs or dark photons would be a peak in the recoil spectrum at the rest mass energy of the particle.
The Super Cryogenic Dark Matter Search (SuperCDMS) experiment [12], located in the Soudan Underground Laboratory in Northern Minnesota, used cryogenic germanium detectors to search for signals produced by dark matter. Particle interactions within the detector produce phonon and ionization signals, which are measured using superconducting transition edge sensors (TES) and charge electrodes, respectively. The ratio of the two signals is different for interactions with nuclei or electrons, providing an efficient discrimination tool [12]. The interleaved layout of the charge and phonon sensors allows for the further discrimination of events near the electrodes, where they can suffer from reduced charge collection [13], giving rise to the name interleaved Z-sensitive Ionization and Phonon (iZIP) detectors. In the CDMS Low Ionization Threshold Experiment (CDMSlite), a much higher bias voltage is applied across the detector to make use of the Neganov-Trofimov-Luke (NTL) effect [14,15], in which additional phonons are created in proportion to the number of drifting charges and the magnitude of the bias voltage. This effect leads to a sensitivity to considerably lower-energy interactions and thus lighter dark matter particles. However, the discrimination between electron and nuclear recoils is lost. Additionally, the increased amplification causes saturation effects to occur at lower energies than in the regular iZIP operating mode, lowering the upper end of the usable energy range.
Traditionally, SuperCDMS dark matter analyses have focused on searches for WIMPs [16] scattering off detector nuclei. In this paper, electron recoil data from SuperCDMS Soudan are analyzed to constrain the parameters that describe the absorption of ALPs and dark photons. The expected signal from this process is a peak in the energy spectrum at the energy corresponding to the mass of the dark matter particle, with a width that is given by the resolution of the detector at that energy. For this analysis, we do not model or subtract the background; thus, only upper limits on the rates of dark absorption can be set. Conservative limits are placed on the coupling of ALPs to electrons and the kinetic mixing between the dark and SM photons for particle masses between 40 eV=c 2 and 500 keV=c 2 , where CDMSlite data cover the lower and iZIP data the higher masses.
An overview of the two dark matter candidate models, the expected signals in the detector, and the assumptions used to determine limits on the couplings are given in Sec. II. In Sec. III the experimental setup of SuperCDMS Soudan is described. Section IV describes the key steps in the data analysis and the method used to derive limits on the coupling parameters, with the results presented in Sec. V and discussed in the concluding Sec. VI.

A. Axions and ALPs
The axion was originally proposed to account for the apparent fine-tuning associated with the lack of chargeparity (CP) violation in the strong interaction [4]. This phenomenon could be naturally explained with a new, spontaneously broken, approximate global Uð1Þ symmetry; the axion would be the pseudo-Goldstone boson associated with this spontaneous symmetry breaking [4]. Other spontaneously broken global symmetries appear in many extensions to the SM (such as string theories [5,17]), giving rise to axion-like particles (ALPs). Axions and ALPs would both feature an effective coupling to the SM photon [17,18] and to electrons [19]. However, the mass of the canonical axion (associated with the strong CP problem) is limited to be less than ∼10 −2 eV=c 2 by the duration of the neutrino signal from SN 1987A [20] and the cooling of neutron stars [21], meaning it is too light to excite an electron in a SuperCDMS detector through the absorption process. These constraints are dependent on the required coupling of axions to nucleons; since ALPs do not require such a coupling, the constraints do not apply. Therefore, only ALPs (and not axions) will henceforth be discussed.

B. Dark photons
Dark photons are hypothetical vector bosons that would mediate a new dark force in models with a new local Uð1Þ symmetry. Kinetic mixing with the SM photon [22] enables the dark photon to interact with electrically charged particles. Although the dark photon does not necessarily need to be accompanied by other particles outside of the SM, it is usually invoked as part of a dark sector, where it serves as a mediator between the dark sector and the SM. The kinetic mixing is then generated by loop-level interactions of much heavier dark particles, charged under both SM electromagnetism and the new Uð1Þ charge. There is no general constraint on the mass of dark photons in these models. However, in order to be a viable dark matter candidate its mass must be less than twice the electron mass; otherwise, its relic abundance would be depleted by decays into electron-positron pairs [11].

C. Dark absorption
The absorption cross sections for ALPs and dark photons relate the observed rate to physical quantities relevant to the particular candidate: namely the axioelectric coupling in the case of ALPs, and the kinetic mixing between the dark and SM photons. These cross sections depend on the photoelectric cross section of the target atom (in this case Ge).
We assume the dark matter candidate is nonrelativistic and constitutes all of the galactic dark matter. The galactic dark matter flux Φ ¼ ρ χ v χ =m χ depends on the local dark matter density ρ χ , mass m χ , and velocity v χ . Throughout this work a local dark matter density of ρ χ ¼ 0.3 GeV=ðc 2 cm 3 Þ [23] is assumed. As we will see, the cross sections in question are inversely proportional to the dark matter particle's velocity and therefore cancel the velocity dependence in the flux, making this search velocity independent.
It should be noted that the predicted dark matter signal rates are an approximation under the assumption that the dark matter Compton wavelength is much larger than the size of the atom [24]. This is a common assumption in present literature for the calculation of sensitivity limits [24][25][26][27][28] and limit projections [3,11]. This assumption requires corrections at dark matter masses above ∼10-100 keV=c 2 [24,29]. More accurate rate predictions taking multi-body and relativistic effects into account require dedicated calculations of dark matter absorption by the atom [30]. Such calculations do not yet exist for dark photons. For ALP absorption in Ge they exist only for masses up to 100 keV=c 2 [19], depending on the target material, which does not cover the full energy range for our search. For this reason, and for consistency with existing publications, the presented analysis applies this imperfect assumption.

Axioelectric effect
The effective ALP-electron interaction is quantified by the axioelectric coupling g ae . Dark absorption of an ALP via the axioelectric effect would eject a bound electron from an atom. The expected cross section σ a is proportional to the photoelectric cross section σ pe [31,32] of the target where E a is the ALP's total energy, β a ¼ v a =c is its relativistic beta factor with velocity v a and speed of light c, α is the fine structure constant, and m e is the mass of the electron. The axioelectric coupling g ae is the parameter on which we set a limit. Under the assumption of nonrelativistic dark matter ALPs (E a ¼ m a c 2 and β a ≪ 1) constituting all galactic dark matter, the event rate for the axioelectric effect [23] can be expressed as a function of the axioelectric coupling g ae :

Dark photon absorption
The kinetic mixing of dark photons to SM photons enables an effective coupling to electrons, and with it the absorption of dark photons by atoms. The cross section for this process [11] is given by where the index V denotes the dark photon, ϵ is the kinetic mixing parameter and β V ¼ v V =c is the dark photon's relativistic beta factor. Under the same assumptions as above (nonrelativistic dark photons which constitute all of the galactic dark matter), the event rate as a function of the kinetic mixing parameter ϵ [11] is given by: While in-medium effects can alter the effective kinetic mixing parameter that is probed through the dark absorption channel, this correction is only necessary for dark photon masses ≲20 eV=c 2 [3], which is below the mass range considered in this analysis.

Photoelectric cross section
The cross sections of interest depend on the photoelectric cross section σ pe of the target, Ge. The σ pe data in the analysis range of 40 eV to 500 keVare expected to be approximately independent of temperature [3]. Discrepancies were found in the literature [33][34][35][36][37][38][39][40][41] across the analysis range, mostly for energies below 1 keV. For this reason we use both a nominal and a conservative photoelectric cross section curve; the nominal curve is used to calculate our main results while we use the conservative one in the estimate of the associated uncertainty. The construction of the nominal σ pe curve follows the approach taken in Ref. [3], with data from Refs. [33,34,35] for photon energies below 1 keV, from (1-20) keV, and above 20 keV respectively. The conservative σ pe curve was determined by using the smallest σ pe values found in the literature search, which results in the largest (most conservative) implied values of ϵ or g ae for a given measured rate. As much of the data found in the literature search were presented in plots, the data were extracted using a digitization web tool [42]. Both the nominal and the conservative photoelectric cross sections are shown in Fig. 1.

III. SUPERCDMS SOUDAN SETUP
SuperCDMS searches for dark matter interactions in cryogenic semiconductor detectors. The SuperCDMS Soudan detectors consisted of cylindrical Ge crystals, 25 mm thick with a diameter of 76 mm and a mass of ∼0.6 kg. The high-purity crystals were instrumented on their top and bottom faces with hundred of tungsten transition edge sensors (TES). The TES arrays were interleaved with biased electrodes, used to collect charge carriers liberated by particle interactions in the detector substrates. The detectors were grouped in 5 stacks (towers, T1 to T5) with 3 detectors each (labeled Z1 through Z3).
The detector towers were located in the innermost of a set of nested copper cans for thermal shielding, surrounded by layers of polyethylene and lead shielding against environmental radioactivity, and a layer of scintillator panels to identify and discard interactions caused by residual cosmogenic radiation; see Ref. [43] for details.
The total measured phonon energy from both primary recoil phonons and NTL phonons (the additional phonons produced by charges moving through the detector in the presence of an applied electric field) is given by where E r is the primary recoil energy, e is magnitude of the electron charge, V bias is the applied bias voltage and ε eh is the average energy required to produce an electron-hole pair (ε eh ¼ 3.0 eV for electron-recoil interactions in Ge [44,45]). In turn, the recoil energy can be expressed in terms of the measured phonon energy E t and the measured charge signal E q as For iZIP detectors a 4 V bias voltage across the detector was used. The sensor layout of the iZIP detectors [13] made it possible to identify and discard events near the top and bottom surfaces, where signals could suffer from reduced charge collection efficiency.
The data used for the iZIP analysis were taken between October 2012 and July 2013, from four of the original seven detectors included in the low-mass WIMP search, described in Ref. [46]. The reasons for this selection are described in Sec. IV B.
For CDMSlite [47], the charge and phonon sensors on one side of a detector were set to a voltage bias of ∼70 V while the sensors on the opposite side were grounded; only the phonon sensors on the grounded side were read out. For the chosen bias voltage, an amplification factor of more than 20 is achieved for the phonon energy from electronrecoil events. The large intrinsic amplification causes the measured signal to saturate for energies exceeding approximately 25-30 keV.
Beginning in 2012, SuperCDMS Soudan operated individual detectors in CDMSlite mode. In total, three temporally separated datasets, referred to as runs, were acquired: Run 1 was a proof of principle with a single detector [47], Run 2 used the same detector for an extended period of time to yield an improved dark matter search result [48], and Run 3-with slightly less exposure than Run 2-was performed with a different detector [49]. During Run 3, a change in operating conditions caused the phonon noise performance to worsen, motivating the decision to separate the analysis of Run 3 into two parts, referred to as Period 1 and Period 2.
Data from Run 2 [48] and Run 3 [49] were used for the CDMSlite part of the analysis discussed here. CDMSlite Run 2 has a lower threshold, a larger exposure and a moderately lower background than Run 3 (see Ref. [49] for a discussion of the difference in backgrounds between the two runs). The main limitation for this analysis is the background, which leads to a similar sensitivity for both runs.
Interspersed throughout the dark matter search, calibration data were taken using 133 Ba and 252 Cf sources. Neutrons from the 252 Cf source activated the detectors, producing 71 Ge that decays via electron capture with a half-life of 11.43 days [50]. The resulting K-, L-, and M-capture lines at 10.37 keV, 1.30 kev and 160 eV are used for energy calibrations, as are the gamma absorption lines from the 133 Ba source.
More details of the experimental setup, the different operating modes, and past analyses for SuperCDMS Soudan can be found in Refs. [12,43,46,47,49,51].

IV. ANALYSIS
To relate the observed event rate in the region of interest to an interaction rate, the detection and event selection efficiencies must be determined. In addition, for dark absorption of nonrelativistic particles, the primary signal is a fixed energy deposition; the expected signature in the measured energy spectrum is a Gaussian peak at the energy corresponding to the candidate particle's mass, with a width given by the resolution of the detector. Therefore, it is also necessary to characterize each detector's energy resolution.
The dark matter mass range under consideration for the CDMSlite data is 40 eV=c 2 to 25 keV=c 2 , where the lower limit is motivated by the limit setting method (see Sec. IV C). In the iZIP analysis the considered masses range from 3 to 500 keV=c 2 , where the lower limit is chosen to avoid the rapidly dropping efficiency at lower energies. For CDMSlite we use the same selection criteria used for the Run 2 and Run 3 WIMP searches, for which the resolution model and detection efficiencies are already published [48,49,52]. Section IVA summarizes these results. Electron recoils in iZIP detectors have not previously been the focus of a dark matter analysis. As such, a reanalysis of iZIP electron recoils was necessary. Section IV B describes the details of the iZIP analysis, including event selection criteria and their efficiencies, the energy scale calibration, and the resolution model.
In Sec. IV C we motivate the definition of the analysis window size, the selection method for detectors to be included in the analysis (for a particular dark matter mass), and the technique used to combine data from different detectors to produce an upper limit on the rate.

A. CDMSlite
CDMSlite Run 2 data were acquired in 2014 between February and November. CDMSlite Run 3 Period 1 and Period 2 data were acquired in 2015 from February to the end of March and from April to May, respectively. The total detector exposures-defined as the product of detector live time and mass-are shown in Table I for each run.

Event selection and signal efficiency
The event selection criteria, and therefore the signal efficiency and associated uncertainties, are the same as those used for the Run 2 and Run 3 WIMP searches [48,49]. The signal efficiency shown in Fig. 2 describes the fraction of all recorded detector events that met all of the data selection criteria at a particular measured energy.
At low energy the efficiency is mostly determined by the trigger [53]. For Run 2, a trigger threshold as low as 56 eV was reached. In Run 3, the trigger rate at low energies was dominated by noise-induced events. These events were removed in the analysis based on their pulse shape, which lowered the efficiency and raised the effective threshold to 70 eV. For events above approximately 100 eV in Run 2 and 200 eV in Run 3, the reduction in efficiency has little energy dependence and largely results from the radial fiducialization which removes events near the edge of the detector where inhomogeneities in the electric field lead to reduced NTL amplification. For a detailed discussion of the systematic and statistical uncertainties on the efficiency, see Refs. [49,51]. The original WIMP search analyses of CDMSlite Run 2 and 3 extended up to 2 keV and the efficiency curves that were only derived up to this energy do not cover the full range needed for this analysis. For CDMSlite Run 2, the efficiency curve was extended up to 30 keV for a background study [52]. This was accomplished by linearly interpolating the efficiency between the values measured at 2 keV and the 71 Ge peak at 10.37 keV; above this energy, electron recoils from 133 Ba calibration data were compared to Monte Carlo simulations, showing a drop in efficiency starting at ∼20 keV, which is attributed to saturation in one of the phonon sensors (see Ref. [52] for details). For this analysis, the efficiency for Run 3 was extended in the same manner. However, the saturation effect observed in Run 2 occurs at a higher energy in the Run 3 detector, leading to a constant efficiency above the Ge K-shell line and below the upper analysis threshold. For both runs, the uncertainty on the efficiency is extended beyond 2 keV in the same manner as the efficiency itself.

Resolution model
We use the resolution model developed in the original CDMSlite analyses [49,51], given by where E is the measured energy, σ E is the baseline noise resolution, σ F describes the impact of the electron-hole pair statistics (accounting for the Fano factor [54]), and σ PD contributes a term to the resolution that is linear in energy and accounts for factors such as position dependence in the detector. The latter two quantities are energy dependent and are parametrized as σ F ¼ ffiffiffiffiffiffi ffi BE p and σ PD ¼ AE, where A and B are constants. The resulting three free parameters of the resolution model, A, B and σ E , are constrained by the observed resolution of the easily identifiable K-, L-, and M-shell electron capture peaks of Ge, as well as the baseline width that results from electronic noise. The resolution of each of these peaks is extracted from Gaussian fits; Eq. (7) is then fit to the widths of these Gaussians (weighted by their uncertainties) at their respective energies, to extract the model parameters and their uncertainties from the fit. These are listed in Table II.
The resulting resolution model is plotted in Fig. 3, separately for CDMSlite Run 2, Run 3 Period 1, and Run 3 Period 2. While the highest energy point used in the fits is the Ge K-shell capture line at 10.37 keV, we assume that the fitted model is accurate for the entire analysis range. More details on the resolution model and a discussion of the uncertainties on the fit parameters can be found in Refs. [49,51]. The upper (lower) uncertainty band is formed by evaluating the resolution model with the bestfit parameters plus (minus) their uncertainties. The resulting 1σ uncertainty band is more conservative than the published bands in Refs. [49,51].

B. iZIP
The iZIP portion of this analysis draws data from the SuperCDMS low-mass WIMP search in Ref. [46], acquired with seven detectors between Oct 2012 and July 2013 [46]. The same live-time selection is used here, which excludes periods directly following the 252 Cf calibrations. Here we exclude data from two of the seven detectors due to shorts on one or more readout channels (phonon and charge), and we exclude data from a third detector that shows evidence  of incomplete charge collection on one side. The four detectors selected for this analysis are listed in Table III along with their exposures.

Event selection and signal efficiency
To select events, we require that they fulfill a set of data quality criteria (data quality cut), that the detector in question issued a trigger (trigger requirement) while no other detector did (single-scatter cut), that the event was not coincident with a signal in the muon veto detector (muon veto cut), and that they pass a fiducial volume cut. The first four criteria are identical to those described in Ref. [46] and have a combined efficiency of around 95% for each detector, with slight energy dependence near the analysis threshold.
The role of the fiducial volume cut is to reject the 210 Pb surface-event background and to remove events near the edge and close to the flat surfaces of the detector where the distorted electric field may lead to a reduced NTL gain. The fiducial volume cut developed for nuclear recoils, and its efficiency, was determined using 252 Cf calibration data in the low-mass WIMP search of Ref. [46]. This cut definition uses the charge signals, which differ significantly between nuclear and electron recoils, and its efficiency is energy dependent. As the signal of interest here consists of single-scatter (and thus essentially pointlike) electronrecoil events, the signal efficiency should be reevaluated using this type of events. Most electron recoils originate from gamma interactions which often scatter multiple times within the detector and thus cannot be used as proxy for the signal events. However, there is one identifiable sample of single scatter events in our dataset: the Ge K-shell capture events at E K ¼ 10.37 keV. Since these events only appear at a fixed energy, they cannot be used to measure the energy dependence of the efficiency of this cut. Therefore, a new fiducial volume cut was developed (based exclusively on the charge signal distribution and therefore referred to as the "charge fiducial volume" or QFV cut) with the goal to make its efficiency largely energy independent so that it can be measured using the Ge K-shell events. This cut consists of two components, a radial charge cut to remove events near the cylindrical surface and a charge symmetry cut to remove events near the flat surfaces. Below we describe each of these cuts and how we assess the energy dependence of their efficiency, before we discuss how the overall efficiency of the combined cut is determined.
Radial charge cut.-The radial charge cut removes events where the charge collected in the outer electrodes exceeds a certain fraction of the total charge collected. 1 This fraction is determined at the Ge K-shell peak in the total charge spectrum (inner plus outer electrode). The signal in the outer electrode is required to be less than three times the baseline resolution of that electrode, as measured using the events that appear at E K in the inner electrode. This new definition of the radial cut fulfills the requirement of a constant efficiency above the K-shell energy, since the energy distribution between inner and outer signal is energy-independent for pointlike events occurring at a given position in the detector. Due to the noise in the measurement, events at lower energy with very little or no charge collected on the outer electrode still have a significant chance to have a reconstructed amplitude above the cut limit and thus fail the cut (see Fig. 4). This effect is quantified by modeling the signal distribution between the inner and outer charge electrodes using the Ge K-shell events. This model is then scaled to lower energies and convolved with the applicable charge resolution before using it to determine the radial charge cut efficiency.  [51] and Run 3 [49]. The best-fit curves (solid lines) and 1σ uncertainty bands (shaded regions) are shown with the measured widths (points, with 1σ error bars too small to see on this scale) of the three Ge electron capture peaks and the baseline noise resolution. The top, middle, and bottom panels show the models for Run 2, Run 3 Period 1, and Run 3 Period 2, respectively. Note that this definition differs from that of the radial charge cut developed in Ref. [46], where events with any discernible signal in the outer electrodes were removed. can thus be removed by requiring that the signal amplitude is similar between the two sides ("charge symmetry"). In Ref. [46] the charge symmetry cut was defined using only the inner charge signals on both sides of the detector. With our new definition of the radial charge cut, the total charge signal on each side of the detector may include a contribution from the outer electrode. Therefore, the charge symmetry cut is redefined for this analysis based on the total charge signal on each side. A cut parameter is defined as the ratio of the difference between the charge collected on each side to the sum of the total charge collected on both sides. For detector-bulk events at a given energy, the distribution of this parameter is roughly Gaussian, centered near zero. The distribution is widest at low energy, before narrowing to an approximately constant width at energies above ∼10 keV. From 1 keV to the K-shell electron capture peak energy, the mean μ and width σ of the distribution are measured in 1 keV energy bins, and events outside of μ AE 3σ are removed. In the range above E K the cut is set to stay constant at the value determined at E K (see Fig. 5).
Below E K , the total number of events just outside the cut boundary is very small, so the final spectrum does not change significantly even if the cut position is loosened far beyond possible uncertainties. This guarantees that the assumption of a constant efficiency below E K is conservative. At higher energies, the event density near the bulk distribution is higher overall so that the exact choice of the cut position has a greater influence on the efficiency, and small non-Gaussianities of the distribution become more relevant. The constant cut value, together with the still slightly narrowing distribution suggests a moderately increasing efficiency. However, in absence of a method to measure the precise efficiency value, we assume a constant efficiency which leads to a conservative limit in the final analysis.
QFV cut efficiency at E K .-The final step is to determine the efficiency of the new QFV cut for the (single-scatter) Ge K-shell events by exploiting the short 11.43 day half-life of 71 Ge [50]. For this we use the data with a particularly high 71 Ge decay rate acquired over 10 or 20 day periods directly after the 252 Cf neutron calibrations that were excluded from the dark matter search. The livetime, data quality, and veto criteria are identical to those used in the dark matter analysis. The approach used is similar to that used in [48]. Events with recoil energy within a AE4σ window of the Ge K-shell electron capture peak are selected, where σ is the energy resolution at E K . The recoil energy is determined using Eq. (6) to ensure that we select all true K-shell events, including those with reduced NTL gain that are removed by the QFV cuts. The selected events are divided into two categories: (1) those passing the QFV cut and with total phonon energy within the signal region (E K AE 1σ, see Sec. IV C below for a discussion on the choice of the signal region), and (2) those failing the former criteria. FIG. 4. Signal amplitude in the outer vs the inner charge channel for side 2 of detector T2Z2 (Q 2;outer vs Q 2;inner ). The cut line that removes events with an outer signal above a given fraction of the total signal is shown. At energies below E K (10.37 keV) the line cuts into the noise distribution, removing some of the events with a real outer fraction below the threshold. This effect is modeled and accounted for in the final efficiency calculation.
FIG. 5. Charge symmetry cut parameter vs total charge signal for detector T2Z2, where Q 1 is the total charge collected on side 1, and Q 2 is the total charge collected on side 2. Bulk events are localized in a narrow band around zero that widens toward low energy. A Gaussian distribution is fit to the band in 1 keV bins at low energy and the position of three times the standard deviation is indicated by black crosses in the figure. The cut line is an exponential function fit to these points below E K with a smooth transition to a horizontal line above E K . The distribution continues to narrow above E K , so a constant cut value together with the assumption of a constant efficiency (as function of energy) will lead to a conservative upper limit on the extracted rate, compensating for any uncertainty that is introduced by the higher density of events near the cut line at high energies.
The number of K-shell events in each category is extracted via a likelihood fit to the time distribution of events, where time is measured from the most recent 252 Cf calibration period. The model probability distribution function Pðt; rÞ is the sum of an exponential (decaying with the 71 Ge lifetime τ) and a constant background component: where t is the time since the last 252 Cf calibration period and r is the ratio between the number of 71 Ge events that are represented by the exponential component and all other events that are represented by the flat component. The efficiency is then the number of 71 Ge events in the signal region divided by the sum of the 71 Ge events in both categories. The uncertainty on the fit result is determined by calculating the likelihood as a function of r and extracting the 1σ confidence interval from the resulting distribution. Note that this efficiency is the combined efficiency of the QFV cut and the signal window selection. This ranges from 30%-36%, depending on the detector.
Combined efficiency.-The combined efficiency of all selection criteria is constant above ∼5 keV for all selected detectors; the uncertainty on this (approximately 3% for all detectors) is dominated by the statistical uncertainty on the QFV efficiency from the likelihood fit. As an example, the total combined efficiency for detector T2Z2 is shown in Fig. 6; all detectors exhibit similar behaviour.

Energy calibration and resolution model
With the selection criteria applied, the events with reduced NTL amplification are removed, and we can calibrate the energy scale for all events using a simple scaling function for the conversion between the measured total phonon energy E t and the recoil energy E r instead of applying an event-by-event correction using Eq. (5). A quadratic function is fit to the ratio of measured phonon signal to true energy for peaks that are used for the iZIP calibration: the Ge K-shell peak at 10.37 keV, the 66.7 keV peak that appears in 252 Cf calibration data as consequence of inelastic neutron scattering, the 356 keV peak from 133 Ba calibration data, and the 511 keV peak from positron annihilation that is observed in dark matter search data. The L-and M-shell peaks used for fitting the resolution model in the CDMSlite analysis are below the analysis threshold of 3 keV imposed for the iZIPs in this analysis (see Sec. IV C below). Also, the baseline noise peak is not used, as at the true energy of 0 keV the ratio of measured to true energy is undefined. Defining the energy scale in this manner accounts for detector saturation at higher energies.
The functional form of the resolution model used for the iZIPs is the same as the one used for CDMSlite [see Eq. (7)]. The model is fit to the resolution of five peaks weighted by their uncertainties: the four peaks used for calibrating the energy scale and the baseline noise peak. The model parameters and the 1σ uncertainty band on the resolution are determined for each detector individually. The upper (lower) edge of the band is formed by taking the upper (lower) value of the 1σ confidence interval for each parameter determined by the fit. The resulting model parameters and their statistical uncertainties from the fit (systematic uncertainties are comparatively negligible) are shown in Table IV. The model for T2Z2, including the measured peak energies and widths to which the model is fit, is shown as an example in Fig. 7.

C. Limit setting
An upper limit on the rate was calculated for the three CDMSlite and four iZIP datasets. This section describes the calculation of the limit on each dataset, and the method by which the limits are combined. No background modeling is performed in this analysis, therefore, only an upper limit on the signal strength can be extracted. The signal model, a Gaussian centered at the dark matter mass with the width determined by the detector's resolution at that energy, is the same for the interactions of ALPs and dark photons. To set a limit on the axioelectric coupling and the dark photon kinetic mixing parameter, we first set a limit on the observed rate of signal events. Limits on the physical quantities of interest are then extracted using Eqs. (2) and (4), for g ae and ϵ respectively. To determine a limit on the interaction rate for a particular dark matter mass, the events in a window around the mass equivalent energy in the spectrum are counted and a Poisson upper limit on that number is calculated, which is then divided by the efficiency-corrected exposure to convert it to a limit on the rate. The process is repeated for each dataset with its  corresponding event spectrum, resolution model, and efficiency, before the different datasets are combined. The choice of window size is a compromise between maximizing the signal efficiency (large window) and minimizing the included background (small window). The optimal window size depends on the observed event numbers. A window size of AE1σ was selected which is close to the optimum for numbers on the order of one to a few tens of events per σ (where σ is the resolution at the respective mass-equivalent energy), which are the typical values found in our datasets.
However, the event selection efficiency drops quickly near the trigger threshold and the expected signal shape (the Gaussian times the efficiency function) is no longer centered about the mass-equivalent energy, but skewed toward higher energies. In this case, setting the cut at þ1σ may remove the dominant fraction of the recorded signal. Therefore, for energies close to the trigger threshold, the upper limit of the window is chosen based on the expected signal shape rather than the primary Gaussian shape. The cut value is the þ1σ-equivalent, cutting the same 15.9% that þ1σ would cut from a Gaussian. The lower edge of the window is kept at −1σ. This change in window choice is implemented for energy depositions in the CDMSlite detector below 100 eV for Run 2 and below 200 eV for Run 3.
Dark matter masses below the trigger threshold can be studied due to the positive tail of the expected event distribution. Since the efficiency estimate for masses far below the trigger threshold is strongly impacted by potential non-Gaussian tails and the uncertainty on the resolution, the lower bound of the CDMSlite analysis range is chosen to be roughly 2σ below the trigger threshold.
For the iZIP detectors we keep the AE1σ window throughout, but we limit the dark matter analysis to an energy range above 3 keV, which is well above the trigger threshold (around 1 keV) so that the discussed effect due to rapidly decreasing signal efficiency is small. Since the best sensitivity at low masses is expected to be derived from CDMSlite data anyway (due to the superior energy resolution), this imposed threshold is not expected to limit the sensitivity of this analysis.
The upper edge of the dark matter analysis range of 25 keV for the CDMSlite datasets is set to avoid saturation effects caused by the large intrinsic signal amplification. For the iZIP detectors, the upper limit of 500 keV is just below the energy of the highest available peak (511 keV) used for calibration and to measure the resolution.
Limits on the event rate are calculated for a set of narrowly spaced discrete dark matter masses. The limits are calculated separately for each dataset which includes the energy range corresponding to each given dark matter mass.
Differences in background rates between detectors call for a method that enables the calculation of a combined limit that is not unnecessarily weakened by datasets with background rates that are truly higher, while still minimizing selection bias from statistical fluctuations when only including low-rate datasets.
For a given dark matter mass the following procedure is implemented: (1) Select the dataset with the lowest observed rate.
(2) Discard any dataset where the difference between its observed rate and the lowest rate exceeds three times the uncertainty on this difference (3σ). (3) Include any remaining datasets in the combined analysis. Given the small number of datasets, a 3σ deviation in the rates will rarely occur if all detectors truly measure the same rate and differences are only due to statistical fluctuations. However, discarding datasets still may introduce a small bias, so a Monte Carlo simulation was performed to quantify this bias. For the case that all detectors indeed measure the same rate, the bias was found to be less than one percent, given the typical observed rates and the number of detectors involved. Since the actual rate is not the same in all detectors, the true bias is lower, and can thus be considered negligible.
The measured event numbers of all included datasets (within the energy window for a particular dark matter mass) are summed and the statistical 90% upper limit on that number is calculated. This upper limit is then divided by the combined efficiency-weighted exposure to determine the 90% upper limit on the measured event rate. The actual number of datasets that are included for a particular dark matter mass varies as a function of energy due to variation of background contributions in each detector; in FIG. 7. The fitted resolution model used in the iZIP analysis, for the example of detector T2Z2. The measured peak widths (points) are used to determine the best-fit curve (solid line) and 1σ uncertainty band (shaded region). The inset plot shows the same data, zoomed in to the region between 0 and 70 keV. some energy ranges only one detector is kept, whereas for other ranges all detectors are included.
With the bias of the selection of datasets to be included in the analysis shown to be negligible, the main remaining systematic uncertainty in the data analysis as described here is the uncertainty on the efficiency of the charge symmetry cut for the iZIP detectors. As discussed above, the uncertainty of the charge symmetry cut efficiency is negligible at low energy compared to the uncertainties already explicitly included in the uncertainty band on the efficiency curve (particularly the statistical uncertainty from the method of determining the efficiency at the Ge K-shell peak, see Sec. IV B 1); at high energy, the choice made in the assumption of the efficiency accounts for possible systematics and leads to a conservative limit. In the CDMSlite analysis the dominant uncertainty is also the uncertainty on the efficiency, where the main contributing factor is again the fiducial volume cut. The uncertainty on the energy resolution is subdominant. Further details of the uncertainties in the CDMSlite analysis can be found in [48,49]. We also assessed the impact of the assumption that single-scatter events are pointlike. Using data from the NIST online database ESTAR [55], we estimated the effect of a finite extension of the highest energy events considered in the analysis to be less than 2% and thus negligible compared to other uncertainties. So, for the final uncertainty on the rate limit, we only include the uncertainty bands on the resolution and efficiency curves as shown in the respective figures above.
For a given mass, the upper (lower) bound for the rate limit is generally determined using the upper (lower) limit on the resolution and the lower (upper) limit on the efficiency. However, the alternative method of determining the analysis window for low masses in CDMSlite can lead to a situation where a wider resolution leads to a better sensitivity. For these masses the upper and lower uncertainty bands are determined by calculating the limits on the rate for all combinations of the upper and lower limits of the resolution and efficiency, and picking the lowest and highest rate limits, respectively.
Three combined rate limits were produced: (1) a combined iZIP limit from the four iZIP detectors, (2) a combined CDMSlite limit, and, in the energy range from 3-25 keV=c 2 , where both CDMSlite and iZIP detectors are used, (3) an overall combined limit using all seven datasets (CDMSlite Run 2, Run 3 Periods 1 and 2, and four iZIP detectors). The calculated rate limits and their uncertainties are shown in Fig. 8.
In the overlap between the CDMSlite and iZIP analysis range (3-25 keV=c 2 ), the CDMSlite rate limit is generally stronger than the iZIP limit. This is expected if background rates are comparable, since in CDMSlite mode the detectors have a better energy resolution (about a factor of 2 at 10 keV, compare Figs. 3 and 7). Additionally, the iZIP data analyzed here were acquired during the first measurement period (late 2012 to mid 2013), while the CDMSlite data were acquired toward the end of SuperCDMS Soudan (early 2014-mid 2015); therefore, some of the cosmogenic activity (e.g., 65 Zn at 9 keV or 55 Fe around 6 keV) has decayed to a significantly lower rate. However, the difference in rate is small enough that the method of selecting detectors for the combined limit still includes iZIP detectors, so the combined limit is typically between that from CDMSlite and iZIPs. There is a data point near 7 keV=c 2 where the combined limit is stronger than the iZIP or CDMSlite limits separately. Here the rate of all detectors included in the combination is similar and the combined limit benefits from the improved statistics.

V. RESULTS
The limits on the rate are converted to limits on the relevant physical quantities, g ae and ϵ, using Eqs. (2) and (4), respectively. The upper (lower) bounds of the uncertainty bands are calculated by combining the upper (lower) bound on the limit on the rate with the conservative (nominal) photoelectric cross section.
An additional systematic uncertainty stems from the measurement of ρ χ , which enters linearly in the calculation of the final limits. However, the results in literature to which we compare our limits (see below) use the same assumption of ρ χ ¼ 0.3 GeV=ðc 2 cm 3 Þ; thus we chose not to reflect this uncertainty in our final results.
The results on the search for ALPs and dark photons are shown in Figs. 9 and 10 in the form of exclusion limits on FIG. 8. Calculated rates from combining data from the CDMSlite runs (light red dotted line), the iZIP detectors (dark red dashed line), and all datasets (black solid line). The shaded region around each line corresponds to the uncertainty on the rate limits. The elevation in the limits near 10 keV, 1 keV, 160 eV and in iZIP data around 9 keV results from the higher measured rates due to the Ge K-, L-and M-captures and the K-capture of 65 Zn respectively. the axioelectric coupling g ae and the dark photon kinetic mixing parameter ϵ, respectively. Limits from other directdetection experiments and astrophysical constraints from models of stellar cooling are also shown for comparison.

VI. CONCLUSION
This analysis sets the strongest laboratory constraint on galactic dark matter ALPs in the mass range ð40-186Þ eV=c 2 . The absolute values of the laboratory limits depend on the assumption that the respective species constitutes all of the galactic dark matter with a local density of 0.3 GeV=ðc 2 cm 3 Þ. Astrophysical constraints on the observed cooling of white dwarf [11,62] and red giant [61,62] stars set the strongest exclusion limits below 1 keV. It should be noted though that a different set of assumptions are used to produce the astrophysical constraints. For example, the production rates of ALPs in a stellar environment requires a precise understanding of the energy levels and occupation levels for each nucleus in the star. In practice this is done with state of the art models of the radiative opacities in the Sun [66].
World leading or competitive limits are also set on the kinetic mixing of dark photons in the mass range of ð40-186Þ eV=c 2 . Astrophysical limits [24] from horizontal branch stars, red giants, and the Sun are weaker than those from direct detection experiments below 1 keV=c 2 .