Improved WIMP-search reach of the CDMS II germanium data
R. Agnese, A.J. Anderson, M. Asai, D. Balakishiyeva, D. Barker, R. Basu Thakur, D.A. Bauer, J. Billard, A. Borgland, M.A. Bowles, D. Brandt, P.L. Brink, R. Bunker, B. Cabrera, D.O. Caldwell, R. Calkins, D.G. Cerdeño, H. Chagani, Y. Chen, J. Cooley, B. Cornell, C.H. Crewdson, P. Cushman, M. Daal, P.C.F. Di Stefano, T. Doughty, L. Esteban, S. Fallows, E. Figueroa-Feliciano, G.L. Godfrey, S.R. Golwala, J. Hall, H.R. Harris, S.A. Hertel, T. Hofer, D. Holmgren, L. Hsu, M.E. Huber, D. Jardin, A. Jastram, O. Kamaev, B. Kara, M.H. Kelsey, A. Kennedy, M. Kiveni, K. Koch, A. Leder, B. Loer, E. Lopez Asamar, P. Lukens, R. Mahapatra, V. Mandic, K.A. McCarthy, N. Mirabolfathi, R.A. Moffatt, S.M. Oser, K. Page, W.A. Page, R. Partridge, M. Pepin, A. Phipps, K. Prasad, M. Pyle, H. Qiu, W. Rau, P. Redl, A. Reisetter, Y. Ricci, H.E. Rogers, T. Saab, B. Sadoulet, J. Sander, K. Schneck, R.W. Schnee, S. Scorza, B. Serfass, B. Shank, D. Speller, D. Toback, S. Upadhyayula, A.N. Villano, B. Welliver, J.S. Wilson, D.H. Wright, X. Yang, S. Yellin, J.J. Yen, B.A. Young, J. Zhang
CCDMS October 15, 2015
Improved WIMP-search reach of the CDMS II germanium data
R. Agnese, A.J. Anderson, M. Asai, D. Balakishiyeva, D. Barker, R. Basu Thakur,
29, 30
D.A. Bauer, J. Billard, A. Borgland, M.A. Bowles, D. Brandt, P.L. Brink, R. Bunker, B. Cabrera, D.O. Caldwell, R. Calkins, D.G. Cerde˜no, H. Chagani, Y. Chen, J. Cooley, B. Cornell, C.H. Crewdson, P. Cushman, M. Daal, P.C.F. Di Stefano, T. Doughty, L. Esteban, S. Fallows, E. Figueroa-Feliciano, G.L. Godfrey, S.R. Golwala, J. Hall, H.R. Harris, S.A. Hertel, T. Hofer, D. Holmgren, L. Hsu, M.E. Huber, D. Jardin, A. Jastram, O. Kamaev, B. Kara, M.H. Kelsey, A. Kennedy, M. Kiveni, K. Koch, A. Leder, B. Loer, E. Lopez Asamar, P. Lukens, R. Mahapatra, V. Mandic, K.A. McCarthy, N. Mirabolfathi, R.A. Moffatt, S.M. Oser, K. Page, W.A. Page, R. Partridge, M. Pepin, A. Phipps, K. Prasad, M. Pyle, H. Qiu, W. Rau, P. Redl, A. Reisetter, Y. Ricci, H.E. Rogers, T. Saab, B. Sadoulet,
39, 45
J. Sander, K. Schneck, R.W. Schnee, S. Scorza, B. Serfass, B. Shank, D. Speller, D. Toback, S. Upadhyayula, A.N. Villano, ∗ B. Welliver, J.S. Wilson, D.H. Wright, X. Yang, S. Yellin, J.J. Yen, B.A. Young, and J. Zhang (SuperCDMS Collaboration) Division of Physics, Mathematics, & Astronomy,California Institute of Technology, Pasadena, California 91125, USA Institute for Particle Physics Phenomenology, Department of Physics, Durham University, Durham, United Kingdom Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA Pacific Northwest National Laboratory, Richland, Washington 99352, USA Department of Physics, Queen’s University, Kingston Ontario K7L 3N6, Canada Department of Physics, Santa Clara University, Santa Clara, California 95053, USA SLAC National Accelerator Laboratory/Kavli Institute for ParticleAstrophysics and Cosmology, Menlo Park, California 94025, USA Department of Physics, South Dakota School of Mines and Technology, Rapid City, South Dakota 57701, USA Department of Physics, Southern Methodist University, Dallas, Texas 75275, USA Department of Physics, Stanford University, Stanford, California 94305, USA Department of Physics, Syracuse University, Syracuse, New York 13244, USA Department of Physics and Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,Texas A&M University, College Station, Texas 77843, USA Departamento de F´ısica Te´orica and Instituto de F´ısica Te´orica UAM/CSIC,Universidad Aut´onoma de Madrid, 28049 Madrid, Spain Department of Physics & Astronomy, University of British Columbia, Vancouver, British Columbia V6T 1Z1, Canada Department of Physics, University of California, Berkeley, California 94720, USA Department of Physics, University of California, Santa Barbara, California 93106, USA Department of Physics, University of Colorado Denver, Denver, Colorado 80217, USA Department of Physics, University of Evansville, Evansville, Indiana 47722, USA Department of Physics, University of Florida, Gainesville, Florida 32611, USA Department of Physics, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801, USA School of Physics & Astronomy, University of Minnesota, Minneapolis, Minnesota 55455, USA Department of Physics, University of South Dakota, Vermillion, South Dakota 57069, USA Department of Physics, University of Florida, Gainesville, FL 32611, USA Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA SLAC National Accelerator Laboratory/Kavli Institute for Particle Astrophysics and Cosmology, Menlo Park 94025, CA School of Physics & Astronomy, University of Minnesota, Minneapolis, MN 55455, USA Fermi National Accelerator Laboratory, Batavia, IL 60510, USA Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Department of Physics, Syracuse University, Syracuse, NY 13244, USA Department of Physics, South Dakota School of Mines and Technology, Rapid City, SD 57701, USA Department of Physics, Stanford University, Stanford, CA 94305, USA Department of Physics, University of California, Santa Barbara, CA 93106, USA Department of Physics, Southern Methodist University, Dallas, TX 75275, USA Institute for Particle Physics Phenomenology, Department of Physics, Durham University, Durham, UK Division of Physics, Mathematics, & Astronomy,California Institute of Technology, Pasadena, CA 91125, USA Department of Physics, Queen’s University, Kingston ON, Canada K7L 3N6 Department of Physics, University of California, Berkeley, CA 94720, USA Pacific Northwest National Laboratory, Richland, WA 99352, USA a r X i v : . [ h e p - e x ] O c t Department of Physics and Astronomy, and the Mitchell Institute for Fundamental Physics and Astronomy,Texas A&M University, College Station, TX 77843, USA Department of Physics, University of Colorado Denver, Denver, CO 80217, USA Department of Physics & Astronomy, University of British Columbia, Vancouver, BC V6T 1Z1, Canada Department of Physics, University of Evansville, Evansville, IN 47722, USA Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA Department of Physics, University of South Dakota, Vermillion, SD 57069, USA Department of Physics, Santa Clara University, Santa Clara, CA 95053, USA (Dated: October 15, 2015)CDMS II data from the five-tower runs at the Soudan Underground Laboratory were reprocessedwith an improved charge-pulse fitting algorithm. Two new analysis techniques to reject surface-event backgrounds were applied to the 612 kg days germanium-detector weakly interacting massiveparticle (WIMP)-search exposure. An extended analysis was also completed by decreasing the10 keV analysis threshold to ∼ c . Afterunblinding, there were zero candidate events above a deposited energy of 10 keV and six events in thelower-threshold analysis. This yielded minimum WIMP-nucleon spin-independent scattering cross-section limits of 1.8 × − and 1.18 × − cm at 90% confidence for 60 and 8.6 GeV/ c WIMPs,respectively. This improves the previous CDMS II result by a factor of 2.4 (2.7) for 60 (8.6) GeV/ c WIMPs.
PACS numbers: 95.35.+d, 95.30.Cq, 85.25.Oj, 29.40.Wk
I. INTRODUCTION
The mass balance of the Universe is the subject of in-tense research and debate. Discrepancies between grav-itationally determined galaxy-cluster masses and theirobserved luminosities provided the earliest motivationfor dark matter [1, 2]. Modern galactic rotation curvessharpen the argument [3–5], as do more recent studies ofgalaxy cluster dynamics [6, 7]. Likewise, spectroscopy ofintergalactic x-ray-emitting gas [6, 7] and gravitationallensing [8–12] elevate these presumed mass discrepanciesto the level of a crisis.To bring the data sets into consistency requires ei-ther modifications of gravity [13, 14] and/or large quan-tities of nonluminous matter. Observations of collid-ing clusters [15, 16] provide evidence for excess nonlu-minous matter, though alternate models of gravity andsome nonluminous matter can apparently also reproducesuch results [17, 18]. Nonrelativistic (cold) relic particledark matter alone could resolve these discrepancies andis also considered an essential ingredient in gravitationalsimulations of the large-scale structure of the Universe.For example, the Via Lactea and Millennium simulationsshow excellent agreement with the observed large-scalestructure of our Universe when cold dark matter is in-cluded [19–21].While the observed galactic dynamics and large-scalestructure naturally lead us to consider cold particle darkmatter, cosmological measurements have been importantin constructing a consistent model for the evolution ofthe Universe using such cold dark matter (CDM). Theaccelerating expansion of the Universe [22–24], big bangnucleosynthesis [25], baryon acoustic oscillations [26], ∗ Corresponding author: [email protected] and the cosmic microwave background [27, 28] supporta cosmology whose dominant components are dark en-ergy (which could correspond to a cosmological constantΛ) and nonbaryonic CDM. When interpreted within theframework of the ΛCDM model, these cosmological mea-surements enable precise determination of the CDM anddark-energy content of the Universe [29].Particle physics provides clues as to the possible iden-tity of nonbaryonic CDM. It was realized early on thatweakly interacting massive particles (WIMPs) with GeV-to TeV-scale masses could thermally freeze out in theearly Universe to give the correct relic density [30]. Su-persymmetry provides a WIMP candidate in the light-est supersymmetric particle and has many other benefits,like solving the hierarchy problem [31]. Still other parti-cle physics models considered more recently [32] providemotivation for light WIMPs and linkages to the mat-ter/antimatter asymmetry of the Universe.WIMPs can be searched for directly through their scat-tering off nuclei in a terrestrial detector. Since these in-teractions are expected to be rare it is important for adark-matter detector to have a low threshold (10 keV orbelow) and excellent background rejection capabilities.The Cryogenic Dark Matter Search (CDMS) collabora-tion has developed cryogenic semiconductor detectors fo-cusing on those properties for the purpose of measuringWIMP-scattering events.This paper explores new analysis techniques using datataken by the CDMS experiment for the direct detectionof WIMPs during the CDMS II running period [33–35].There are four data periods associated with the data setused in this work, varying from 1 to 6 months in du-ration. The experiment used Ge and Si detectors andwas located in the Soudan Underground Laboratory at ashielding depth of 2090 meters water equivalent (m.w.e.).The largest payload consisted of a total of 19 Ge ( ∼
240 geach) and 11 Si ( ∼
110 g each) detectors. All detectorsare cylindrical, 7.6 cm in diameter, ∼ x Z y where x (1–5) is the tower number and y (1–6) indicates the position within the stack (from topto bottom). FIG. 1. (Color online) Isometric representation of theCDMS II detector arrangement and tower occupation, withdirection of north indicated. The ZIP detectors in each towerare numbered 1–6 and are either silicon (yellow) or germa-nium (brown). The identification numbers of each detectorwere given to track the raw material the detectors were pro-duced from and the processing to which they were subjected.For example the detector “G9” is a germanium detector withidentification number 9. Green stars indicate detectors thatwere used in this data analysis (see Sec. II).
In each detector, particle interactions produceelectron-hole pairs (ionization) together with phonons.The charge carriers are drifted by a small electric field(3 V/cm for Ge and 4 V/cm for Si) and collected on con-centric aluminum electrodes deposited on one of the flatfaces of the crystal. On the opposite face, supercon-ducting transition-edge sensors (TESs) arranged in quad-rants collect phonons before thermalization. Charge andphonons are measured independently for the purpose ofevent-by-event discrimination: WIMP signal events willproduce a nuclear recoil (NR) in the detector, whereasmost background processes produce an electron recoil(ER). The ratio of ionization to phonon signal amplitudesallows discrimination of NRs from the far greater numberof ERs with a rejection factor > [33] in the 10–100 keVregion. Lead shielding further reduced gamma-inducedbackgrounds, and neutron-induced NR backgrounds were reduced with polyethylene shielding. The shielding issurrounded by an active scintillator veto to tag eventsinduced by cosmic-ray muon showers.Recoils within ∼ µ m of the detector surface can havereduced ionization signals because of poor charge col-lection, which can be sufficient to misclassify a surfaceER as an NR. These reduced-charge ERs originate fromthree major sources: electrons produced by β emitterson or near the detector surfaces, electrons ejected fromphotons scattering in nearby material, and photons thatinteract near the detector surfaces. All of these mecha-nisms can lead to events in which an energy depositionis observed in only a single detector, producing “surfaceevents” that can be mistaken for single-interaction NRs.The long-lived Rn daughter
Pb, implanted in thedetector surfaces and their copper housings, is the pri-mary source of such surface events. This class of eventsconstitutes the dominant background for the CDMS IIWIMP search and presents a considerable challenge [37].Fortunately, recoils near detector surfaces are character-ized by prompt phonon absorption in the TESs. Con-sequently, the phonon signals for surface ERs are (onaverage) faster than for NRs in the detector’s bulk, en-abling surface-event background discrimination based onphonon-pulse timing. The event selection criteria derivedfrom phonon timing (called “timing cuts”) presented hereare essential to obtaining optimal WIMP sensitivity forthe CDMS II data set because they mitigate the surface-event background, improving the overall ER rejection to > .Nonsignal NR events have two known sources: Ge re-coils induced by scattering neutrons and Pb nucleifrom the decay of
Po near the detector surfaces. Theexpected neutron background for CDMS II is roughlyequal parts cosmic-ray muon-induced neutrons and radio-genic neutrons from trace contaminants in the shieldingand detector. Radiogenic neutrons produced outside theshielding have too little energy to cause a detectable NRafter penetrating the polyethylene. By combining simu-lations with in situ data, we estimate the neutron back-ground to be subdominant compared to the surface-eventbackground (see Secs. VII B and VII C). The
Pb back-ground contribution is also subdominant in this analysis.The data set used in this work comes from the Ge de-tectors in the final CDMS II five-tower exposure and wasacquired between July 2007 and September 2008. Thetotal raw exposure for this running period was approxi-mately 612 kg days. The original analysis of this data setprovided world-leading sensitivity to spin-independentelastic WIMP-nucleon scattering in 2010 when it was firstpublished [33]. Here, we reevaluate this data set using im-proved data reduction algorithms and surface-event rejec-tion methods that reduce the expected surface-event con-tamination (“leakage”) in the WIMP signal region com-pared with the original analysis. In the 2010 analysis a10 keV [38] analysis threshold was used to limit the ex-pected background to less than one event for the entireexposure. Intriguing results at low WIMP mass [39–43]motivated us additionally to examine the data with re-duced thresholds.
II. CDMS II DETECTOR PROPERTIES
The CDMS II Z-sensitive ionization and phonon de-tectors (ZIPs) are operated at a temperature of ∼
50 mKand feature six readout channels: two charge electrodeson one side and four phonon sensors on the oppositeside [44]. The analysis undertaken here includes only Gedetectors, for which the ionization channels are biasedwith a 3 V potential across the crystal. Furthermore, fiveof the 19 Ge detectors were omitted from the analysisbecause of readout channel failures across all data sets.This gives a direct analogy to the original analysis of thisdata set, which used the same detector subset. Ioniza-tion channels are read out by a low-noise junction field-effect transistor circuit with an operating temperature of ∼
150 K [45]. Phonons are detected by quasiparticle-trap-assisted electrothermal-feedback transition-edge sensors(QETs). The QET signal is amplified by supercon-ducting quantum-interference devices (SQUIDs) thatare thermally coupled to the cryostat’s 600 mK coldstage [45]. Figure 2 shows a schematic of the ZIP chan-nel layout. The outer charge electrode acts as a vetoagainst events that deposit energy near the sidewalls,where higher background is expected and charge collec-tion is more likely to be incomplete. !"
FIG. 2. (Color online) a) Inner and outer electrodes for mea-suring the ionization signal. The inner electrode extends fromthe center to a radius of 34.5 mm and the outer electrode ex-tends from a radius of 35.5 mm to just before the edge (0.5–2 mm; the diameter of the detectors is 76 mm). b) Phononchannels (QETs) A–D arranged in quadrants on the oppositeside of the detector. Each phonon channel is composed ofapproximately 1000 TESs formed from thin films depositedonto the crystal surface and then photolithographically struc-tured. The TESs are connected in parallel to form the QETsand each QET is read out by its own circuit containing aSQUID array.
Data taken during the CDMS II experiment are com-posed of calibration (using a
Ba gamma source, or a
Cf neutron source) and WIMP search. The
Ba cali-bration data were used to determine the energy scale andto study various systematic effects. The
Cf calibrationdata were used to study detector response to NRs and actas a proxy for determining WIMP acceptance.The CDMS II front-end electronics issued experimentaltriggers in response to activity in either the surroundingscintillator veto or the ZIPs. For veto-triggered events,the front-end electronics required two or more veto pan-els to have coincident signals in excess of their hardwarethresholds. ZIP-triggered events occurred when a com-posite phonon signal for any ZIP exceeded its preset dis-criminator threshold. The composite phonon signal wasthe analog sum of the four phonon signals from the ZIPwith a 900–18000 Hz band-pass filter applied to reducenoise-like fluctuations. During
Ba calibration runs,selective readout was employed because of the high rate;i.e. only detectors in the tower in which a ZIP triggeroccurred were read out. During
Cf calibration andWIMP-search runs, all detectors were read out in re-sponse to either veto or ZIP triggers. Veto-triggeredWIMP-search events were recorded for use in studies ofcosmogenically induced neutron backgrounds. Each trig-gered event includes two charge and four phonon traces(per ZIP readout). Each trace consists of 2048 digitizedamplitudes acquired at a rate of 1.25 MHz, with 512 sam-ples prior to the trigger time.Charge carriers moving in the electric field of the de-tector generate phonons with a total energy proportionalto the potential difference traversed (Neganov-Luke ef-fect [46, 47]). These “Luke phonons” are typically bal-listic and add to the phonon signal generated by the pri-mary recoil. The recoil energy E r was constructed fromthe total phonon energy E p by subtracting a term pro-portional to the ionization-derived recoil energy E q (or“ionization energy”) to account for the Neganov-Luke ef-fect: E r = E p − eV(cid:15) E q , (1)where e is the elementary charge, V is the absolute valueof the operating potential of the detector, and (cid:15) is the av-erage electron-hole pair creation energy. In Ge the valueof (cid:15) is approximately 3.0 eV/pair for ERs [48]. The “ion-ization yield” (or yield) is defined as the ratio of ion-ization energy to recoil energy ( E q /E r ) and we requirethis quantity to be unity for ERs (see below). An NR willproduce less ionization than an ER of equal recoil energy.This effect is well known [49] and provides the basis forthe ZIP detector’s primary method of ER/NR discrim-ination. For recoil energies between 10 and 100 keV inGe, the NR ionization yield varies between 0.2 and 0.3.Using the Ba calibration data, the charge andphonon amplitudes were calibrated so that the energiesof known gamma lines are reproduced and the yield isunity for ERs. Specifically, the charge channels were cal-ibrated using the 356 keV line after a ∼
10% correction tothe ionization amplitudes to account for small systematicvariations with interaction location (within the detector).We do not fully understand the origins of the systematicvariation with interaction location, though it is empir-ically robust. One possibility is that the detectors arenot neutralized uniformly by the infrared (940 nm) LEDthat is activated between data-taking periods to removetrapped charge. The relative phonon-channel calibrationwas performed by scaling the phonon-channel amplitudefractions–the amplitudes of individual channels dividedby the sum of the amplitudes–to have equivalent distri-butions. Finally, the summed phonon energy was cali-brated by requiring the ionization yield to be unity onaverage for ER events in the 65–100 keV region.
III. RAW DATA REDUCTION AND THEREPROCESSING
Our analysis parameters are calculated from digitizedcharge and phonon pulses for each ZIP detector usinga set of pulse-processing algorithms. These algorithmsdistill timing and amplitude information from the fourphonon and two charge traces (per ZIP). Surface eventsare rejected using timing quantities derived from theraw data, such as the rise time τ and the delay t del ofthe largest of the four phonon pulses with respect tothe faster ionization signal. The amplitude of the innercharge-channel pulse gives a measurement of E q and theamplitude of the analog sum of the four phonon pulsesgives a measurement of E p after calibration. These en-ergy variables are used to construct the ionization yield( y ) and recoil energy ( E r ) described in Eq. (1).The relative charge-channel amplitudes (charge “par-tition”), relative phonon-channel amplitudes (phonon“partition”), and relative phonon pulse timing provideinformation about an event’s position within the detec-tor. These additional parameters are also calculated us-ing the output of the pulse-processing algorithms and areused for the phonon event-position-based correction andcharge-derived fiducial-volume restrictions.This section gives an overview of the most importantalgorithms for constructing these analysis parameters. Italso explains the upgrade to the charge-pulse processingthat motivated the reprocessing of this data set. A. Parameter extraction
Timing parameters . Timing estimators are derivedusing an algorithm that steps along a low-pass-filteredtrace to identify the pulse rise time (RT) and fall time(FT). We call this the RT-FT-walk algorithm. The tracesare first filtered with a low-pass Butterworth filter [50].Two sets of parameters are produced; one where thecutoff frequency is 50 kHz and one where it depends onthe signal-to-noise of the trace. The filter removes high-frequency noise and effectively smooths the pulse for animproved determination of the RT and FT at various percentages of the pulse maximum. RT and FT informa-tion is determined by “walking” along the filtered tracestarting at the maximum and identifying the times atwhich the respective threshold levels are reached. Forexample, it infers the time at which the rising pulse edgereaches 20% of the pulse’s maximum amplitude. The RT-FT-walk algorithm was applied to all charge and phonontraces. The rise time τ is computed as the 10–40% timespan along the rising edge of the largest of a detector’sfour phonon pulses using the RT-FT-walk algorithm. Optimal filtering . A pulse-template optimal filter(OF) [51] is used to produce the best resolved energyquantities. The OF has superior energy resolution com-pared to pulse-integral quantities and is well suited to theanalysis of small recoil energies [52]. A template for theexpected pulse shape is fit to the pulse in Fourier space,deweighting the frequency bins with high noise. Thefrequency deweighting is done for each individual data“series” – a data-taking block normally lasting between10 and 12 h–using noise power spectral densities (PSDs).The PSD is constructed from randomly triggered tracestaken before the series to sample the noise environ-ment. Phonon pulse templates are two-exponential func-tional forms, with rise- and fall-time parameters tuned tomatch individual detectors by fitting to an average pulse.Charge pulse templates were produced empirically by av-eraging normalized data traces. This fitting procedure isdone for all possible pulse template delays and the best-fit delay is chosen.A single-pulse-template OF is performed on eachphonon pulse separately and on the sum of a detector’sfour phonon traces. For charge pulses, to account forcrosstalk between the inner and outer channels, the twopulses are fit simultaneously with a crosstalk-correctingOF (OFX) described in more detail below. During theOF fit an array of best-fit amplitudes is produced thatcorresponds to each possible delay of the template withrespect to the experimental pulse. The maximum suchamplitude in a preselected delay window is chosen andits time defines the delay quantity. For an OF usinga single pulse template this maximum amplitude choicealso has the lowest χ in the preselected window. TheOF amplitude for the summed phonon trace gives E p af-ter calibration and the inner charge amplitude from theOFX procedure gives E q after calibration. The charge-to-phonon delay t del is determined from the differencebetween the OFX charge delay and the 20% crossing ofthe largest phonon signal as determined by the RT-FT-walk algorithm. The delay is computed as follows: t del = t p − t OFX , (2)with t OFX being the OFX ionization delay relative to theglobal trigger time and t p being the 20% crossing timeof the largest phonon signal relative to the global trigger. B. Parameter corrections
Phonon position correction . The phonon sensor re-sponse is highly position dependent [53]. In order to opti-mize the derived phonon variables, a position correctionwas used to modify energy, timing, and yield quantitiesby normalizing to the mean of nearby events. We firstdefined averaging neighborhoods using a five-dimensionalmetric that included total phonon energy, two relativephonon delay parameters, and two phonon energy par-tition parameters [54]. For each parameter to be cor-rected, a lookup table containing the neighborhood aver-ages for each five-dimensional bin was then constructedusing bulk ER events confirmed to be of good qualityfrom
Ba calibration. To apply the correction we cre-ated a corrected version ( v c ) of a given parameter ( v ) asin Eq. (3): v c = v (cid:104) v (cid:105) bin (cid:104) v (cid:105) global , (3)where the subscript “global” refers to averaging over thewhole calibration data set, and “bin” refers to averag-ing over one neighborhood. After the lookup table wasapplied to the Ba calibration data, the selection ofbulk ERs was refined to create the final correction ta-ble, which was then applied to the WIMP-search data.The final energy, timing and yield quantities were mademore uniform across the detector for bulk events by us-ing this procedure, thus improving the rejection of surfaceevents.
Charge crosstalk . Charge signal quantities use theOFX algorithm to account for crosstalk. Using separatepulse-template OFs for the inner and outer electrodes,a single delay was chosen. In the original analysis ofour data set [33], this delay corresponded to the delaythat made the sum of the two OF amplitudes maximal.While this is likely to be very close to the delay thatminimizes the χ of the fit, it is not guaranteed becauseof the combined fitting of the inner and outer electrodes. C. Reanalysis motivation
After finalizing the first analysis of this data set in2009 [33] we noticed that one of the WIMP candidateevents had an OFX delay value that did not correspondto a global minimum in the χ of the fit (see Fig. 3). Thisissue was one of the main reasons for proceeding with areanalysis.The reason for implementing the maximum-amplitudemethod instead of a full χ minimization in the origi-nal analysis was to save computing time; a global triggerissued by one detector forces the readout of all detec-tors, giving a substantial number of traces with near-zeroamplitudes that need not be processed with the time-consuming χ minimization. The result of using thismaximum-amplitude method is to smear the ionization-to-phonon delay distributions slightly, which can cause s] µ delay [ − − − do f / n χ f i t FIG. 3. (Color online) The OFX χ as a function of delayfor the T3Z4 WIMP candidate from the original analysis [33].The delay chosen by the original ionization fit to this event(red dashed line) is not a global minimum in χ , rather it isdisplaced ∼ µs from the global minimum (green dot-dashedline). While this small error in the inferred delay for this andother events did not affect the measured energy significantly(a ∼
2% decrease), it caused an error in the surface-event-rejection timing parameter for low-energy events. additional background events to leak past the timingcuts. This delay smearing affects lower-ionization eventsmore than higher-ionization events. In particular, NRsthat deposit energy in a single detector (“single scat-ters”) are more susceptible than the (higher-ionization)low-yield ERs used to quantify the expected leakage ofsurface events into the signal region.The OFX inaccuracy made only a small contribution tothe previously published limit, and background estimateswere adjusted upward in the original result to accountfor the resulting larger uncertainties at low energies [33].Nevertheless, the number and character of the candidateevents provides useful information on the possibility of asignal. It is very likely that the candidate event featuredin Fig. 3 would have been removed by the timing cut inthe original work had the ionization delay at the globalminimum of the χ been chosen, but using the improvedOFX procedure also has the potential to cause previouslyexcluded events to become WIMP candidates.For the reasons stated above, the data were repro-cessed, selecting the global minimum of the χ in theOFX algorithm for most ionization pulses, rather thanthe maximum-amplitude method. Because of the largenumber of traces in the raw data consistent with noise,traces with pulses corresponding to charge energies be-low a detector-dependent threshold (0.94 keV ee [38] onaverage across the 14 detectors used in this work) werestill processed with the maximum-amplitude method inorder to save processing time. This has no effect on theresults presented here because this energy is below ourlowest analysis threshold. IV. DATA SELECTION AND EFFICIENCIES
Once the data set was calibrated and position cor-rected, we created several selection criteria (or “cuts”) toproduce the cleanest signal-region sample possible. Sincethe rejection of nonsignal events is typically not perfect,a signal event retention efficiency was computed for eachcut. The cuts and efficiencies are covered in this sec-tion, starting with the efficiency for the trigger–a selec-tion criterion that is made in hardware before the eventsare recorded. The result is a signal-region sample corre-sponding to a well-known exposure.A trigger efficiency for each detector was determined asa function of phonon energy and converted to a functionof recoil energy using our measured NR ionization yieldfrom
Cf data. A given detector’s trigger efficiencywas calculated using all WIMP-search events in whichanother detector caused an experimental trigger. Whenthe global trigger initiates a readout of all the detectors,other delayed instances where phonon pulses were in ex-cess of the corresponding discriminator thresholds wererecorded into a logical buffer. Regardless of the con-tent of these trigger records, for each event and detectorthe optimal-filtering techniques described in Sec. II wereused to reconstruct the total phonon energy ( E p ) fromEq. (1). The trigger efficiency curve is then the ratio oftwo recoil-energy spectra: the spectrum of events with aphonon trigger in both the detector in question and an-other detector divided by the spectrum of events with aphonon trigger in another detector. As shown in Fig. 4–combined across the 14 Ge detectors considered here–thetrigger efficiency is ∼ ∼ recoil energy [keV] e ff i c i e n cy χ FIG. 4. (Color online) Detection efficiency combined acrossall detectors as a function of energy for different classes ofcuts. The curves show the total efficiency as more cuts areadded beyond the hardware trigger (pink triple-dot-dashed):event-level data-quality cuts (green solid), ionization-basedfiducial-volume cut (magenta dotted), ionization-yield cut(maroon dot-dashed), and 5d- χ (see Sec. V C) phonon tim-ing cut (blue dashed). The falloff of the data-quality cutefficiency below ∼
20 keV is due to the interplay between theionization threshold and the requirement that the ionizationyield be 3 σ below the ER band (see main text). The verticalblack dashed line is our standard 10 keV analysis threshold. in bins of recoil energy. Well-motivated functional formswere fit to the bin-wise estimates, with the best-fit re-sults shown in Fig. 4. Most cuts have little effect onthe acceptance of nuclear recoils. The ionization-basedfiducial-volume cut and the phonon timing cuts causethe greatest loss of signal acceptance. Of the three tim-ing cuts described in Sec. V, use of the “5d- χ ” timingcut results in the highest final acceptance above 10 keV,about 50% at ∼
35 keV. The 5d- χ timing-cut efficiency isshown in Fig. 4 as an example while the other timing-cutefficiencies are detailed in Sec. V E. A. Live-time cuts
Live time was removed during periods with disabledreadout channels, poor detector neutralization (charac-terized by increased levels of charge trapping in the crys-tal bulk [52]), decreased resolution, improper experimen-tal configurations, and trigger anomalies that consistof isolated bursts of events or incorrect phonon triggerthreshold. The Soudan Underground Laboratory alsohouses a neutrino detector, MINOS, to measure proper-ties of the “NuMI” neutrino beam originating from FermiNational Accelerator Laboratory [55]. While it is veryunlikely for CDMS to observe any beam-induced events,the loss of live time incurred by removing all time pe-riods coincident with a NuMI beam spill is negligible.The NuMI neutrino-beam cut is implemented as a re-striction on the live time, removing events within 60 µ sof the arrival of the neutrino beam’s 10 µ s spills. Afterapplication of these cuts, 612 kg days of total exposureremain. We define this as our “raw exposure.” Out ofa total of 36% loss of live time (64% data-taking effi-ciency) the largest contributors were bad environmen-tal configurations ( ∼ ∼ ∼ ∼ ∼ B. Quality cuts
The event-level class of quality cuts consists of sev-eral components. A “glitch” cut removes events thathave phonon pulses resembling electronic noise in thephonon readout chain. These events are typically cor-related across multiple detectors and are characterizedby phonon pulses with fall times shorter than the timescale expected for phonon dissipation. Such events areless likely to appear in the ionization readout chain; con-sequently, they are effectively removed with a cut thatrejects events in which the phonon pulse multiplicity issignificantly larger than the ionization pulse multiplic-ity [56]. The muon-veto cut uses the 2 in. thick scintil-lator panels that surround the CDMS II experiment andthe following two rejection criteria: an event was removedif 1) a ZIP event is associated with a muon-like energydeposition (0.58 V or ∼ µs before and 20 µ s after theZIP trigger time, or 2) there is any veto activity abovethe scintillator-panel hardware threshold ( ∼ µ s before a ZIP trigger. Both criteria areused for the muon-veto cut because it makes the cutstronger and has only about a 2% efficiency loss acrossall recoil energies. Finally, because WIMPs will not in-teract in more than one detector within a given event, wedefined a “single-scatter” cut to select events involving asingle ZIP detector as follows: 1) the phonon signal mustbe greater than six standard deviations above the meanelectronic-noise level in the detector under consideration;and 2) the signal in all other detectors must be withinfour standard deviations of the means of their respec-tive noise distributions. Ionization channels were usedfor multiple-scatter rejection for those detectors with de-graded phonon channel performance. Thus, Ge and Sidetectors that are not part of the WIMP-search expo-sure were still live with respect to identification and re-jection of multiple-scatter events. The single-scatter cutefficiency for each detector was estimated from the frac-tion of randomly triggered events (i.e., electronic noise)that satisfy criterion 2). The combined efficiency of the glitch, muon-veto and singles cuts varies by detector andover time, and ranges from 96% to 98%. Another aspectof the singles cut is its use to select certain samples fromcalibration data. For example, multiple-scatter eventsthat have their secondary scatter in the detector above orbelow the triggering detector are said to be “face tagged.”For these events the primary recoil is biased toward thedirection facing the multiples tag. When using the mul-tiples in background estimations, separating the eventsby face helps decrease the systematic uncertainties (seeSec. VII A). Sometimes events with multiples toward thephonon readout side are called “phonon-side” and thosewith multiples toward the charge readout side are called“charge-side.”Some quality cuts depend explicitly on the kinematicvariables. One of these is an ionization threshold cut,which requires events to have reconstructed ionizationsignals greater than 4.5 standard deviations above themean noise level. This cut removes nearly all eventswith zero charge collected (those from very near the sidewalls), resulting in < χ value of the OFX fit to beless than an energy-dependent threshold. This cut se-lects events with high-quality charge energy estimators,suppressing those with excess electronic noise or the oc-currence of multiple pulses within a single event trace(referred to as “pile-up”). Its efficiency was measured asa function of charge energy using Ba calibration dataand translated into recoil energy using the average NRionization yield measured from
Cf calibration. Above ∼
20 keV ee [38] the efficiency is constant and > Cf cal-ibration data, including a small correction ( ∼
13% max-imum across the energy range) accounting for residualleakage of ERs into the NR signal region in these data.Another ∼
5% correction was applied across the whole en-ergy range to account for detector self-shielding to neu-trons and multiple scattering; this correction was basedon a Monte Carlo simulation of neutron scattering.
C. Yield cut
An energy-dependent cut on ionization yield is used asthe primary method for discriminating NRs from ERs.An NR “band” was derived for each detector by fittingthe distribution of NR yields in
Cf calibration with aGaussian hypothesis in bins of recoil energy. The collec-tions of best-fit Gaussian means and widths were then fitwith energy-dependent functional forms, giving smoothparametrizations versus recoil energy of the average yieldand the yield resolution for NRs. The functional formfor the means was inspired by the Lindhard theory [57]( y = a · E br , where a and b are fitted parameters), whereasthe widths are fit to a power law below a fitted energythreshold and a constant above [58]. ER bands were sim-ilarly constructed from Ba calibration. The primaryionization-yield cut requires that events be located withinthe ± σ width of the NR band. By construction, theselection efficiency is ∼
95% and roughly constant withenergy (see Fig. 4). Variations in detector response overthe course of the WIMP search caused slight variationsin each detector’s NR and ER bands. Consequently, thebands depend on both detector and time. The time vari-ation is represented by widths of the band lines shown inFig. 5, where they are plotted with the
Cf data fromwhich the NR bands are derived.To prevent ERs from entering the signal region at lowrecoil energy where the ER and NR bands overlap (seeFig. 5), the yield-based discrimination is refined to in-clude a requirement that candidate events lie at least 3 σ below the mean of the ER band (see Sec. IV D). Thiscondition and the ionization threshold have the greatestimpact on the overall NR detection efficiency at low en-ergies. The efficiencies for these two cuts are 100% forrecoil energies greater than ∼ σ σ recoil energy [keV] i on i z a t i on y i e l d σ σ FIG. 5. (Color online) ER and NR bands defined in theplane of yield versus recoil energy for a representative detector(T1Z2). The ER band mean ± σ curves (red hatched), NRband mean ± σ curves (blue fine cross-hatched), ionizationthreshold 3.0 σ and 4.5 σ curves (green coarse cross-hatched),and the 10 keV threshold line (black dashed) are shown super-imposed onto Cf data (points). The widths of the bandsrepresent the time-period variation of the fits that are usedto define them.
D. Data regions
To assess the effectiveness of the timing cuts definedin Sec. V, several standard data regions in the nontim-ing variables were defined. This terminology will be usedthroughout the rest of this work to describe data sam-ples used for timing-cut tuning and sideband-style back-ground estimations. All of the event samples include thedata-quality and fiducial-volume restrictions. Further,they include the requirement that events have chargeenergies greater than 4.5 standard deviations above themean electronic-noise level.1.
Nuclear-recoil single scatter (NRSS).
Eventsthat are below the ER band mean -3 σ line, withinthe NR band mean ± σ lines, and are single scat-ters. The portion of the NRSS events that pass thetiming cut and are in the range of 10–100 keV makeup the WIMP signal region. The 5d- χ timing cutuses a slightly modified signal region in that thewidth of the NR band is optimized. We use theterm NRSS there as well, leaving the precise widthof the NR band to be determined from context.2. Nuclear-recoil multiple scatter (NRMS).
Events that are below the ER band mean -3 σ line,within the NR band mean ± σ lines, and are notsingle scatters.3. Wide-band (WB).
Events that are below the ERband mean -5 σ line, and above the NR band mean+2 σ line.4. Wide-band multiple scatter (WBMS).
Eventsthat are below the ER band mean -5 σ line, abovethe NR band mean +2 σ line, and are not singlescatters. These events are a good representation ofsurface events. V. PHONON TIMING DISCRIMINATION
ZIP detectors have an excellent ability to discriminateNRs from ERs if the energy depositions occur away fromthe surfaces of the detector. But background surfaceevents can populate the NR band, despite being ER innature, since they can have reduced ionization yield. Weremove surface events from our WIMP candidate sampleusing the timing characteristics of the phonon signals. Asmall number of surface events, however, can survive intothe signal region. We call these leakage events. The def-inition of the timing cut affects the surface-event leakageand the total exposure (through the NR acceptance), so acrucial element in the timing-cut construction is tuningfor optimal WIMP-detection sensitivity. Three timing-cut constructions are reviewed in this section.Having three independent methods for surface-eventrejection gives a handle on the systematic uncertaintyof the leakage estimates. The WIMP limits for these0three realizations of surface-event rejection and sensi-tivity maximization are presented in Sec. VIII. Table Igives a summary of the different timing parameters usedfor each timing cut construction: the “classic,” “neural-network,” and “5d- χ ” analyses. The choice of param-eters used in particular timing-cut constructions is ex-plained in the following sections. Quantity Description Analysis τ VF phonon rise time NN, 5d- χ t del VF phonon delay NN, 5d- χ ˜ τ CF phonon rise time Classic, 5d- χ ˜ t del CF phonon delay Classic, 5d- χ ˜ τ CF phonon 40–70% rise time 5d- χ ˜ w CF phonon pulse width NN P phonon 50–70 kHz power NNTABLE I. Brief description of each of the phonon timingquantities used in this work. The abbreviation “VF” indi-cates the use of variable-frequency filtering prior to the appli-cation of the RT-FT-walk algorithm. “CF” indicates a con-stant (50 kHz) filter prior to that algorithm (see Sec. III A).Under “Analysis,” “NN” refers to the neural-network timinganalysis. Each of the timing cuts was optimized to produce thebest expected sensitivity to WIMPs given the expectedleakage. Since we do not know the WIMP mass, a “tar-get” value is chosen for each analysis (60 GeV/ c formost analyses in this work) and the expected sensitiv-ity is maximized given that WIMP mass. The spectrum-averaged exposure (SAE) is a way to quantify the amountof the raw exposure (MT = detector Mass × live Time)that is utilized toward the WIMP search over the analy-sis energy range, given a WIMP recoil spectrum for mass m χ of f ( E r ; m χ ). The SAE is computed as follows:SAE( m χ , E l , E h ) = MT (cid:82) E h E l dE r (cid:15) ( E r ) f ( E r ; m χ ) (cid:82) E h E l dE r f ( E r ; m χ ) , (4)where E r is the recoil energy, E l ( E h ) is the lower (up-per) signal-region energy limit, and (cid:15) ( E r ) is the cumu-lative signal acceptance efficiency for all cuts at a recoilenergy E r . Note that the SAE is equal to the raw expo-sure only when the analysis efficiency is unity over theentire energy range and that only SAEs with the sameWIMP-mass assumptions and signal-region energy rangeare strictly comparable. Sometimes the right-hand sideof Eq. (4) with MT divided out is called the spectrum-averaged efficiency. The SAE is computed on a detector-by-detector basis and then summed in the final analysis.Before looking at the events in the final signal region(see Sec. VI on “unblinding”) we used the expected leak-age and the SAE of each timing cut optimization to calcu-late the expected sensitivity. The expected sensitivity iscomputed by using the expected leakage and calculatingthe upper limit of counts at the 90% C.L. This is normal-ized by the SAE to produce the lowest WIMP rate theexperiment is sensitive to. For the 10 keV threshold anal-ysis the 5d- χ timing cut had the best sensitivity, while our lower-threshold analysis showed the classic method tohave the best sensitivity. The main results of this workare therefore the limits of the 5d- χ analysis for WIMPmasses above 11.3 GeV/ c , and the limits of the classicmethod for WIMP masses below 11.3 GeV/ c . A. Classic timing analysis
The phonon timing cut strategy that was used inthe original analysis of these, as well as earlier CDMSdata [33, 36, 59] was also used in our reanalysis and pro-vides a point of comparison between the two.Our “classic” timing parameter is defined as the sum oftwo quantities: the delay ( t del ) and the 10–40% rise time( τ ), both derived from the most energetic phonon signalamong the four sensors (see Sec. III). The sum is approx-imately the optimal combination of these two variables,as can be seen in Fig. 6. A timing cut is defined as a setof detector-dependent thresholds on the distributions ofthis parameter, below which all events are rejected. Thethresholds were determined by an optimization schemethat approximately maximizes the sensitivity to a WIMPwith a mass of 60 GeV/ c . Practically, this was accom-plished by maximizing the WIMP-search exposure (asmeasured with Cf NRs) while keeping the total leak-age approximately equal to a “target” leakage of ∼ c WIMP sensitivity while keeping the total expected leak-age well under one event [60].The expected surface-event leakage was estimated fromrepresentative
Ba calibration and sidebands in theWIMP-search data that are insensitive to WIMPs. Thesurface-event background estimates for the three timing-cut strategies are described in detail in Sec. VII A. Al-though the quality cuts remove most unusual events,data-reconstruction artifacts occasionally result in eventswith extreme kinematic quantities (“outliers”). To pre-vent such outliers from skewing the timing-cut optimiza-tion, a consistency cut rejects events for which τ + t del isgreater than 32 µ s or τ - t del falls outside the 0.5% and99.5% quantiles of the combined Ba and
Cf datasets. Figure 6 shows the classic timing cut in the de-lay versus rise-time plane, and Fig. 7 shows the cut inthe yield versus timing-parameter plane for an exampledetector [61]. Applying this timing cut results in a to-tal SAE summed over detectors of 220 kg days between10 and 100 keV for a 60 GeV/ c WIMP and an expectedsurface-event leakage of 0.64 (+0.17 -0.15) stat events, cal-culated after unblinding. The post-unblinding calcula-tion is generally more accurate because it makes use ofthe previously sequestered (see Sec. VI) NR single scat-ters that failed the timing cut (see Sec. VII).1 s] µ rise time [ s ] µ d e l ay [ FIG. 6. (Color online) Classic timing cut in the delay versusrise-time plane for a representative detector (T1Z5).
CfNRSS and NRMS events (blue open circles) and
Ba WBMSevents (red crosses) are shown. Accepted events lie within the“consistency” region (black dashed lines) and to the upperright of the discrimination cut (black solid line). Consideringthe
Ba WBMS events (red crosses), 2381 events fail thetiming cut and four events pass. s] µ timing parameter [ i on i z a t i on y i e l d FIG. 7. (Color online) Classic timing cut in the ionization-yield versus timing-parameter plane for a representative de-tector (T1Z5).
Cf NRSS and NRMS events (blue open cir-cles),
Ba WBMS events (red crosses), and
Ba events inthe ER ± σ band (black filled circles) are shown. Acceptedevents lie to the right of the timing parameter line (blackdashed line) and in the NR band (black solid lines). Consid-ering the Ba WBMS events (red crosses), 2381 events failthe timing cut and four events pass.
B. Neural-network timing analysis
A neural-network technique was used to develop a tim-ing cut using four timing parameters: the previouslydefined phonon delay t del (i) and rise time τ (ii); thephonon pulse width ˜ w (iii), defined as the time differencebetween the 80% points on the rising and falling edgesof the largest-amplitude phonon pulse; and the spectralpower of the largest phonon pulse P (iv), integratedbetween 50 and 70 kHz [62]. These parameters were cho-sen because they showed the most promising discrimina-tion in their one-dimensional distributions.These four variables were fed into a principal compo-nent analysis [63]. Principal component analysis is a sta-tistical method for determining a unitary transformationthat takes N possibly correlated input vectors and re-turns N output vectors that are linear combinations ofthe input vectors ( N = 4 in this case). The output vec-tors are ordered by their statistical variance, so that theoutput vector with the i th-highest variance is called the i th principal component. Since the input vectors canhave different characteristic scales and are dimensional,the input vectors were normalized to zero mean and unityvariance so that the ordering of statistical variances of theoutput vectors is meaningful in an absolute sense.Neural-network computational complexity scalespoorly with the number of input parameters. Therefore,only the first two principal components (i.e., thosewith the highest variance) were selected as inputs forthe neural network. Given input parameters withsimilar intrinsic resolution and physical relevance, thehigh-variance combinations will be those that maximallyseparate the distinct populations (i.e., NRs and ERsurface events). Principal components were selectedseparately for each detector and for several bins oftime spanning the WIMP-search data set. The latter isnecessary to capture changes in detector performancecaused by variations in operating conditions. Theprincipal component rotation showed that all fourtiming parameters contribute significantly to the firstand second principal components in most cases. Thisgenerally indicates that the use of the extra parameters(as compared to the classic analysis) is beneficial eventhough in the end the sensitivity change is not dramatic(see Sec. VIII A).The neural network that was used is a multilayer per-ceptron with one hidden layer, 30 neurons, and a logisticsigmoid activation function. The NETLAB package [64]for MATLAB was used to perform this analysis. Train-ing samples for surface events and NRs were selectedfrom the Ba and
Cf calibration data, respectively.Events from the NR training sample were assigned a tar-get output value of 1, while surface events were assigneda target value of 0. Data were sorted into two bins ofrecoil energy, above and below 30 keV, in order to takethe energy dependence of the timing parameters into ac-count. Finer energy binning was not possible because ofsmall statistics in the training samples.2Separate neural networks were trained for each com-bination of detector, time period, and energy bin us-ing the standard back-propagation-of-errors training al-gorithm [64]. Once the neural networks were trained,they assigned a numerical value to each event in therange 0–1. An output value close to 1 (0) correspondsto events of NR-like (surface-event-like) character. Thedistributions of this output parameter for the calibrationdata, as seen in Fig. 8, were then used to set a thresholdfor each detector, time period, and energy bin.The thresholds were set such that the WIMP-searchexposure is maximized using a target surface-event leak-age of 0.5 events. As in the classic analysis, thistarget leakage approximately maximized the 60 GeV/ c WIMP sensitivity and the same optimization procedurewas used. Following optimization, the total SAE for a60 GeV/ c WIMP is about 216 kg days, with an expectedleakage of 0.87 (+0.24 -0.21) stat events. Similar to theleakage estimate for the classic timing cut, this expectedleakage was estimated after unblinding for better accu-racy (see Sec. VII).
C. 5d- χ timing analysis The 5d- χ surface-event rejection [65] was imple-mented by differentiating events based on a goodness offit to two event-type hypotheses [66, 67]: χ for NR, and χ for surface events. Five timing quantities were usedto form each χ value. Three of the quantities are mea-sures of rise time for the largest-amplitude phonon chan-nel: two based on the 10–40% rise time ( τ and ˜ τ ); andone based on the 40–70% rise time (˜ τ ). The remain-ing two quantities correspond to the delay of the phononpulse relative to the prompt charge pulse – one computedusing a variable and one a constant-frequency RT-FTwalk ( t del and ˜ t del ). The inclusion of different measuresof the same physical quantities (delay and rise time here)increases the robustness of the χ value and is more ef-fective at identifying outliers. Therefore, the delay andrise-time parameters with good one-dimensional discrim-ination and the least redundancy (correlation) were cho-sen.Event samples of neutrons and surface events takenfrom calibration data were used to constrain the timing-quantity distributions for each event type. For eachdetector, calibration data were separated into neutron,charge-side surface-event, and phonon-side surface-eventsamples. The energy-dependent means, µ ( E r ; α )–a vec-tor of the timing-quantity distribution means for eachevent type α –were then fit to the empirically motivatedfunctional form µ ( E r ; α ) = a ( α ) + a ( α ) E r + a ( α ) (cid:112) E r , (5)where the a i ( α ) are free parameters for each timing quan-tity (the vector indices) and for particle type α . The √ E r term is observed to improve the fit. The covariance ma- first principal component − − − − − sec ond p r i n c i p a l c o m pon e n t − − − − − transformed neural-network output − − − − c oun t s FIG. 8. (Color online) Neural-network training distributionsof a representative detector (T1Z5), for the low-energy neural-network bin ( <
30 keV). (Top panel) Contours of constantneural-network output (grayscale solid curves, see text) in theplane of the first two principal components.
Cf NRSS andNRMS events (blue open circles) and
Ba WBMS events(red crosses) are also shown. (Bottom panel) Distributionsof the “transformed” neural-network output: the distributionof NRs from the
Cf data mentioned above (blue solid);and the surface-event distribution from the
Ba data men-tioned above (red dashed). The cut threshold (vertical blackdashed line) with the WIMP-search data passing (green as-terisk) and failing (black × ’s) the timing cut are also plottedalong the top of the plot. The transformed neural networkis used to make the distribution separation more visible; it iscomputed by using the inverse of the neuron response func-tion: φ : ( −∞ , ∞ ) → (0 ,
1) with φ ( x ) = 1 / (1 + e − x ) trix σ ( E r ; α ) was similarly fit using σ ( E r ; α ) = b ( α ) + b ( α ) E r , (6)where the b i ( α ) are matrices of the free parameters foreach pair of timing quantities and particle type α . The3functional form of the variance was motivated by notingthat in a simple model of a pulse with a linear rise butconstant rise time, the probability for a noise fluctuationbefore the pulse rises above the noise is inversely propor-tional to the slope, and therefore is proportional to theinverse of the amplitude (energy).The χ α was then formed for every event according tothe following formula: χ α ( E r ) = ( ξ − µ ) T · ( σ ) − · ( ξ − µ ) . (7)Here ξ is the vector embedding the five timing variablesfor each event and the dependence on E r and α has beenleft implicit on the right-hand side.The surface-event goodness-of-fit variable χ was con-structed for each event by the definition χ ≡ min( χ p , χ q ) , (8)where χ q and χ p are the charge-side surface-event andphonon-side surface-event goodness-of-fit variables re-spectively.Two restrictions were set in the plane of χ versus χ that together complete the definition of the 5d- χ surface-event rejection cut (see Fig. 9). First the potential WIMPevents were required to have χ ≤ c i , where c i is avalue for the i th detector, set by requiring that 90%of the calibration-data neutrons pass the cut. To dis-tinguish potential signal events from surface events, itwas required that χ − χ ≥ η i ( e ) where the index i indicates the detector number and the parameter e in-dicates the event’s energy bin (10–20 keV, 20–30 keV, or30–100 keV).We parametrize–using adjustable parameters t and b –the ionization-yield ( y ) restriction as µ nr − bσ nr ≤ y ≤ µ nr + tσ nr , (9)where µ nr ( E r ) and σ nr ( E r ) are the energy-dependentmean and standard deviation of the ionization yield forNRs, as found using neutron calibration data. The pa-rameters t and b were required to be the same for alldetectors and energy bins. The values η i ( e ), t , and b were determined by a simultaneous optimization for alldetectors and energy bins that maximizes the total SAEfor 60 GeV/ c WIMPs.Optimization of the 5d- χ timing cut and the NR banddefinition was based on requiring the best overall ex-pected sensitivity. For a given timing cut (the set { η i } )and NR band definition, the expected sensitivity is con-structed by dividing the 90% Poisson upper limit on theexpected leakage by the SAE. For an individual detectorwe denote the leakage as L i ( η i ) and the SAE as S i ( η i ).These are smooth functions computed by fitting leakageand exposure evaluations at discrete η i using Ba cali-bration data.Because the cut is defined such that decreasing η i willloosen the restriction, both the leakage and the SAEfor each detector are monotonically decreasing–leading χ B χ recoil energy [keV]
10 20 30 40 50 60 70 80 90 100 N χ - B χ FIG. 9. (Color online) Distributions of calibration data inthe 5d- χ parameters for a representative detector (T1Z2).(Top panel) Cf events in the NR ± σ band (blue opencircles) are shown with the consistency cut that retains 90%of neutrons (vertical black dashed line) in the χ versus χ plane. (Bottom panel) Events passing the consistency cut–NRs from Cf data (blue open circles) and surface eventsfrom
Ba data (red crosses)–in the χ − χ versus recoil-energy plane. The energy-dependent 5d- χ timing cut is alsoshown (black solid line) and events above the line pass the cut.Note the cut is tighter at low energies and looser at higherenergies. to the condition that the slopes d S i /d L i for the opti-mum cut are equal. If not, a unit increase in the leakageof a detector with larger slope could be offset by a unitdecrease in a detector with a smaller slope with a netincrease in exposure, improving the sensitivity. The tim-ing cut optimization was done with respect to this slope,which parametrizes the { η i } uniquely.For the optimum set of timing parameters { η i } the NRband definition is selected by choosing the values of t and b that optimize the overall expected sensitivity.The yield and timing cuts that optimize the expected4sensitivity for a 60 GeV/ c WIMP produced an asym-metric NR band cut with b = − . t = 1 .
8. Thesensitivity optimization gives a total expected leakage of0.5 events and a total SAE of 250 kg days given a WIMPmass of 60 GeV/ c . The estimated leakage after unblind-ing is 1.19 (+0.23 -0.21) stat events. D. Extended analyses
There has been growing interest in low-mass WIMPsearches because of some intriguing published results [39,40] and the suggestion that the baryon asymmetry is re-flected in the dark matter sector [32]. While many pre-vious WIMP searches–guided by the SUSY neutralinoparameter space–paid much attention to ∼
100 GeV/ c WIMP masses, it is interesting in light of these new re-sults to examine data with techniques optimized for muchlower WIMP masses ∼
10 GeV/ c . For this reason thetiming-cut constructions considered so far were extendedwith lower thresholds (down to 5 keV for some detectors),in order to improve sensitivity to low-mass WIMPs.All of the timing cuts presented so far were optimizedfor a 60 GeV/ c WIMP mass, and the correspondinganalyses were restricted to a recoil-energy threshold of10 keV. For events above this threshold the best estimatesof surface-event background leakage are about one event.While the surface-event leakage is worse for the extended-threshold analyses, for low WIMP masses sensitivity im-proves considerably because of the steeply rising WIMPspectrum. Simple extensions to the three timing analyseswere accomplished by lowering the thresholds, confirm-ing the cut efficiencies below 10 keV, and reevaluating thesurface-event leakage estimates. The analysis region forthe extended analyses is approximately 5–15 keV (thresh-old differs by detector; see below), since recoils of lightWIMPs ( (cid:46)
10 GeV/ c ) with energies greater than 15 keVare very rare.For the extended analyses there are several smallchanges to the event selections that help to maximizethe sensitivity to WIMP masses below 10 GeV/ c .A careful study of the leakage induced by lowering thecharge threshold showed that for light WIMPs with sim-ilar cross sections to the CDMS II silicon result [39] anet gain is obtained by lowering the charge threshold to3.0 standard deviations above the mean noise value. Wetherefore modified the charge threshold to this value forthe extended analyses. The recoil-energy threshold wasthen effectively set by the condition that signal eventshave an ionization yield below the mean -3 σ ER bandline. Of course, signal events were still required to bewithin the mean ± σ NR band (with the usual slightmodification for the 5d- χ analyses). Finally, the signalregions were taken to range from threshold up to 15 keVsince this is the region that is expected to add signifi-cantly to the low-WIMP-mass exposure. These changesare implicit where we use the abbreviation “NRSS” inthe context of the extended analyses. Most of the timing cuts (except the 5d- χ , see below) did not have their deci-sion boundaries reoptimized for lower WIMP masses; thethresholds were simply extended and the changes abovewere made.The 5d- χ method, in addition to lowering the thresh-olds, was partially reoptimized for a 8 GeV/ c WIMPmass. The definitions of the χ variables ( χ , χ ) re-mained the same. Additionally, the means and covari-ances were taken to have the same functional depen-dences given in Eqs. (5) and (6). The optimized yieldcut was an asymmetric cut, with parameters t = 1 . b = − .
9, the same as the regular 5d- χ analysis.The reoptimized sensitivity for the 8 GeV/ c was not sig-nificantly better than the standard 5d- χ analysis, andnot as good as the sensitivity for the classic extended-threshold analysis. Therefore, this reoptimization wasnot carried any further and the 60 GeV/ c optimized ver-sion is used. E. Timing efficiencies and summary
Energy-dependent efficiencies for the timing analyseswere computed using neutrons from the
Cf calibra-tion data. Figure 10 presents the combined efficienciesof all cuts (quality, fiducial volume, yield, and timing)for each timing analysis. Above a 10 keV threshold, the5d- χ analysis has the best sensitivity and the highestSAE for a 60 GeV/ c WIMP (see Fig. 10 caption). Theclassic analysis has the second-best SAE, and providescontinuity with the original analysis of this data set [33].The neural-network analysis yields the smallest SAE, buthas the best efficiency just above a 10 keV recoil energythreshold. For the extension of the analysis to below10 keV, the classic timing analysis has the best sensitivity(though not the highest efficiency at low recoil energies;see Fig. 10). As mentioned in Sec. V, the 5d- χ is our“primary” method above 10 keV and the classic is theprimary extended-threshold method.The efficiency function is a necessary ingredient forproducing the limits and any uncertainty on this functionis also present in the final limit. The trigger efficiency un-certainty is ∼
1% across the energy range; this is mostlystatistical uncertainty. Our quality cut efficiency is cal-culated based on baseline noise levels using large eventpopulations and so has negligible uncertainty. Efficien-cies of the other cuts are measured by selecting neutronpopulations in
Cf data and observing the decrease inthe population by the application of the cuts in bins ofrecoil energy. The results are then fit with an empiri-cal functional form with a low-energy falloff similar tothe error function. Fiducial-volume, ionization-yield andphonon-timing cuts each contribute about a 5% statisti-cal uncertainty. This measurement method is, however,prone to error due to the fact that neutrons can havemultiple scatters inside a detector, while WIMPs cannot.In the case of the fiducial-volume cuts we used a MonteCarlo simulation to find a 5% discrepancy for multiple5 recoil energy [keV] e ff i c i e n cy recoil energy [keV] e ff i c i e n cy FIG. 10. (Color online) Total combined efficiencies for all tim-ing analyses. (Top panel) Efficiencies for the 10 keV thresholdanalyses. The SAEs (60 GeV/ c WIMP) are 219.1 kg days forthe classic analysis (black solid), 216.4 kg days for the neural-network analysis (red dot-dashed), and 262.3 kg days for the5d- χ analysis (blue dashed). (Bottom panel) The efficienciesfor the extended-threshold analyses with the addition of the5d- χ analysis optimized for a 8 GeV/ c WIMP mass (ma-genta dotted). All of the analyses have a total exposure (be-fore efficiency reductions) of 612.2 kg days. scatters, and corrected for it. Overall we assign a conser-vative 3% systematic uncertainty on the fiducial-volumeand yield cuts for the effect of multiple scattering. Basedon these estimates, using the averaged sizes of each effi-ciency, we expect the total uncertainty on the efficienciesto be ∼ ∼
10% at the lowest WIMP mass(6.26 GeV/ c ) and drops to ∼
7% by a WIMP mass of7 GeV/ c . VI. UNBLINDING
A blinding technique was used to avoid bias in the set-ting of data selection cuts. The current work used whatcan be referred to as “hidden signal box” analysis [68]which is common in rare-event searches where the signalregion is known a priori . The same data were analyzedpreviously [33], but since all the data were reprocessedwith an upgraded charge reconstruction algorithm, andthe cuts were optimized solely based on calibration dataand distributions outside the newly masked signal region,this is a good approximation of a hidden signal box anal-ysis. Signal events were hidden by removing the single-scatter events in the NR yield band. This technique iseffective for removing bias in the timing-cut preparation,but restricts the information that can be used for estimat-ing surface-event leakage before unblinding has occurred.For this reason the leakage estimates were done beforeand after the unblinding, but consistency between thetwo methods is checked. Including the post-unblindingversion using the NR singles that fail the timing cut madethe final estimate more robust.Upon unblinding, the following number of events passall cuts above 10 keV: zero for the 5d- χ timing analy-sis; two for the classic timing analysis; and one for theneural-network analysis. One of the two candidates fromthe original analysis is a candidate in both the classictiming and the neural-network analyses [33]. The secondcandidate from that analysis, whose poor charge-pulsefitting prompted this reanalysis (see Sec. III), failed alltiming analyses by a substantial margin. Informationabout the two passing candidate events is shown in Ta-ble II. Figure 11 shows the location of these events withinthe signal region of the classic analysis. Detector Recoil energy [keV] Yield AnalysisT1Z5 12.30 0.33 NN, Classic, 2010T2Z3 10.81 0.33 ClassicT3Z4 15.35 0.26 2010TABLE II. Information about the three WIMP candidateevents above 10 keV. Under “Analysis,” “NN” refers to theneural-network timing analysis and “2010” refers to the origi-nal analysis of these data [33]. Note that no candidate eventsare observed for the 5d- χ analysis. Our extended-threshold analyses gave a wide range interms of the number of candidate events. The classicanalysis had six candidates and the neural-network anal-ysis had 16. No events were observed in the extendedversion of the 5d- χ analysis. A description of the can-didate events for our primary extended-threshold anal-ysis (the “classic” cut) is given in Table III. The largenumber of candidate events in the neural-network anal-ysis is attributed to an increased leakage of anomalouslylow-ionization events, those that are normally below theionization threshold in the 10 keV analysis [62]. This in-creased leakage of essentially zero-ionization events wasnot expected and is presumably due to a bias in the train-6 s] µ normalized timing parameter [ − − − − − no r m a li z e d y i e l d − − s] µ normalized timing parameter [ − − − − − no r m a li z e d y i e l d − − FIG. 11. (Color online) Distributions of WIMP-search eventsin the detectors containing the two candidate events for theclassic timing analysis: T1Z5 (top panel) and T2Z3 (bottompanel). The normalized yield is the distance from the averageNR yield in units of the standard deviation of the yield distri-bution; the black horizontal lines on the right side of the plotindicate the ± σ NR band. The normalized timing parameteris the standard timing parameter ( t del + τ ) minus the value ofthe cut boundary. WIMP-search events (red crosses) to theleft of the timing cut (vertical black dashed line) pass all cutsexcept phonon timing and yield; events to the right of thatline pass the timing cut, and those that are also in the NRband (solid black lines) pass the yield cut. The highlightedevent (blue open circle surrounding red cross) in the top panelis candidate 1 and the highlighted event in the bottom panelis candidate 2 (see Table II). All the events shown are in therecoil energy range 10–100 keV. ing set of the neural network. It has negligible contri-bution above the original 4.5 standard deviation chargethreshold. The low number of events in the 5d- χ analy-sis is not unexpected because that cut has an independentenergy bin at 10–20 keV, where the cut is rather stringentbecause of increasing leakage at low energy and the op- timization to exposure at a 60 GeV/ c WIMP mass.
Detector Recoil energies [keV]T1Z5 3.45, 5.73, 12.30T2Z3 10.81T4Z4 7.56T4Z5 7.25TABLE III. Information about the six WIMP candidateevents for the “classic” extended-threshold analysis.
The post-unblinding leakage estimates can be found inTable IV along with the analysis energy ranges, candi-date numbers, exposures and WIMP-mass optimizationassumptions.
VII. ESTIMATED BACKGROUNDSA. Surface electron-recoil background
Surface electron recoils originate from several sources:1) particles emitted from β emitters contaminating thesurfaces of the detector and the material around it (no-tably Pb), 2) photo-electrons emitted from materialneighboring the detector through Compton scattering orthe photoelectric effect, and 3) photons that interact inthe detector within a few microns of the surface. Pho-tons from category 3) can be low-energy xrays or high-energy photons that Compton scatter in the detector.Past studies showed that the dominant contributions arefrom
Pb and photon-induced backgrounds, which con-tribute approximately equally. No other sources werefound to be statistically significant.The expected number of surface events leaking intothe WIMP signal region was calculated using the num-ber of single scatters in the NR band that are rejected bythe timing cut and the surface-event rejection efficiency.The surface-event rejection efficiency was estimated fromthree independent event sets and combined to improveaccuracy. To reduce systematic uncertainties, leakage es-timates were calculated on a detector-by-detector basis(index i below) and where possible the relevant event setwas separated into bins (index j below) of energy andapproximate event position (see “face” bins below). Inevery case the leakage can be expressed as: n = (cid:88) i,j N i s ij m ij M ij , (10)where n is the total expected number of surface eventsleaking into the WIMP signal region. The symbol N i is the number of NRSS events in the NR band rejectedby the timing cut for the i th detector. M ij and m ij arethe number of multiples in (or around; see below) theWIMP signal region failing and passing the timing cut,respectively. The s ij are the fractions of N i in subset j ,which are calculated using the surface-event multiples inthe NR band for the WIMP-search data.7 Method Energy range [keV] Exp. leakage Candidates WIMP mass [GeV/ c ] SAE [kg days]Classic 10–100 0.64 +0 . − . +0 . − . χ +0 . − . > +0 . − . > +0 . − .
16 60 190extended 5d- χ > +0 . − . χ > +0 . − . aa SAE depends on an assumed WIMP mass [see Eq. (4)]; we use the optimization mass in all cases. For this reason SAE is onlycomparable in situations of common optimization mass and signal-region energy range
TABLE IV. Expected leakage and exposure statistics for all of the surface-event rejection methods described in this work. Forthe extended-threshold analyses, we quote 5 keV as an approximate lower limit on the signal region. The actual thresholddepends on the detector and is set by the crossing of the ER band limits and the charge threshold curves (see Sec. V D). Thesymbol > is used to indicate lower limits on the expected leakage. In those situations the event sets that are typically used toestimate the leakage do not have events all the way down to the signal-region threshold (see text). The 14 detectors used in this analysis were split intotwo detector sets according to their positions in thetower: the 12 interior detectors, and two “endcaps,” onthe top or bottom of a stack of six detectors. Equa-tion 10 was applied to each detector set separately be-cause endcaps require additional systematic uncertaintycorrections (discussed below). The M ij and m ij are de-termined independently for three different event sets: 1)NRMS in WIMP-search data, 2) WBMS in WIMP-searchdata, and 3) WBMS in Ba calibration data. For rep-resentative surface-event populations drawn from outsidethe NR band, the index j specifies one of six subsets con-structed from two detector faces (see Sec. IV B on facetagging) and three energy bins (10–20 keV, 20–30 keV,and 30–100 keV). In the case where M ij and m ij weredrawn from the NRMS WIMP-search data, the six eventsubsets were combined for each detector because of thelow statistics, reducing the index set for j to one ele-ment. Each of the three independent estimations of m ij and M ij are used with Eq. (10) to obtain three estimatorsfor n , which were then statistically combined to producethe final leakage estimate.With these estimators, Monte Carlo (MC) simulationsand Bayesian inference were used to estimate the surface-event background and its uncertainties after unblindingthe data [66]. Instead of simulating the posterior distri-bution of n (Eq. (10)), for each event set, MC simulationsare run for the posteriors of the individual Poisson counts m ij and M ij , and the multinomial fractions s ij . Jaynespriors, p ( ρ ) ∝ ρ c with c ≈ −
1, were chosen for the Pois-son distribution. This prior is generally considered an“objective prior,” having the advantage of ensuring in-variance of statistical inference under transformations ofthe Poisson mean [69]. A uniform prior was used for themultinomial distribution. The posterior of n was thencalculated using Eq. (10) with the simulated posteriorsof each component.Systematic uncertainties of different origins were es-timated independently, added in quadrature, and thenincorporated into the posterior of the leakage with a pro- file of the standard normal distribution. Systematic un-certainties from two sources were estimated, for all thedetector and data-set combinations, including the choiceof the prior exponent c and the difference between singlesand multiples.End-cap detectors can have their multiples tagged onlyon one side. This biases the counts of singles upward inthese detectors. Together with the low number of sur-face events in the data, it is challenging to estimate thesurface-event background. As a result, a conservativeestimate where we allow the s ij to be biased higher forendcaps–the detectors with the worst leakage–is used. Asa double check, the component leakages for each detec-tor obtained by the Bayesian approach were compared tothose obtained by the frequentist approach [70]. Thereis good agreement between the two.The final surface-event leakage estimates incorporat-ing both systematic and statistical uncertainties are dis-played in Table IV for all timing cuts presented in thiswork. The analyses with the best expected sensitivitiesgive leakage estimates of 1.19 (+0.23 -0.20) events (5d- χ for the 10 keV threshold analysis) and > σ line and above the NR mean +2 σ line (WBMS) is oneof our typical background samples, but because the ERmean -5 σ and NR band lines cross there is an implicitenergy threshold in this background sample. The signalregion, however, is taken to be the NR band above the 3 σ charge threshold and below the ER mean -3 σ line. Thisleaves a small region (see Fig. 5) that is unaccounted forin the background estimates. The estimate is neverthe-less useful because 1) the implicit background-estimationthreshold can be as low as ∼ B. Cosmogenic background estimate
To estimate the cosmogenic neutron background, acombination of simulation and data was used [71]. Parentmuons that intersect the CDMS II 5-cm-thick scintilla-tor panels are vetoed at nearly 100% efficiency. Our datacontains only a small number of neutrons that accompanyidentified muons, so a high-statistics simulation was usedto estimate how many veto-coincident neutrons produceNRs in the Ge ZIP detectors, compared with those thatare not accompanied by veto activity.Simulation muons generated using the MUSUN [72]program based on slant-path data from Soudan2 an-gular muon flux measurements [73] were used as inputto a GEANT4 MC model of the Soudan CDMS II ex-perimental setup. These muons have a mean energyof ∼
215 GeV and azimuthal and zenith distributionscharacteristic of the overburden of the Soudan site anddepth (2090 m.w.e.). The simulated muons were gen-erated with the appropriate angular distributions andspectra on a five-sided parallelepiped (no floor). Theywere then propagated via GEANT4 through 10 m of rockinto the Soudan hall. All secondaries and the parentmuon(s) were tracked through a complete geometry ofthe CDMS II shielding and detector towers. In the sim-ulations, the CDMS setup is located asymmetrically 2 mfrom one wall of the 8 m east-west cavern dimension toreproduce albedo effects. The statistics correspond to 66live years. A multiple is defined in the simulation outputas an NR in a Ge detector accompanied by energy depo-sition of any sort above 2 keV in any other ZIP detector.In an effort to balance the statistical and systematicuncertainties, three complementary estimators for the to-tal number of unvetoed single-scatter events were con-structed. For the first two estimators, the simulationwas used to calculate a veto ratio (unvetoed to veto-coincident NRs) that was then normalized to the veto-coincident data from the WIMP-search measurements toestablish the background estimate. The second estimatormakes use of the higher-statistics multiple-scatter sam-ple, whereas the first uses only single scatters. The thirdestimator uses only simulated unvetoed single-scatterevents, scaled to the correct experimental live time, andaccounts for the detector efficiencies for each timing anal-ysis.The veto ratio is 0.008 ± ± ± ∼
26% of the vetoedevents should be singles, and the rest multiples. Of thevetoed NR events, five pass the 5d- χ timing cut andfour pass the classic or neural-network timing cuts. Forthe data-driven estimators, we normalize to those passingthe timing cut, rather than introduce additional system-atic uncertainty by applying an average timing-cut effi-ciency to all 14 detectors in the simulation. However, thenumbers obtained either way are consistent, with trade-offs between systematic uncertainty and statistical uncer-tainty in each case. We decided to quote the simulation-driven estimate because it is consistent with the othersand offers the best statistical uncertainty. The system-atic uncertainty is taken as the spread between all theestimators.Based on the 16 unvetoed NR singles over 66 liveyears of simulation data, the 5d- χ timing-cut cosmo-genic neutron background is 0.021 ± stat ± sys events for the 10 keV analysis. The classic and neural-network cuts give 0.019 ± stat ± sys events and0.018 ± stat ± sys events respectively for the10 keV threshold case.The estimates were also made for the extended-threshold analyses by using the energy range 2–20 keVto approximate the extended-threshold energy ranges.Based on 12 unvetoed NR singles in that range over66 live years the 5d- χ extended timing-cut cosmogenicneutron background is 0.009 ± stat ± sys events,the classic is 0.012 ± stat ± sys events, and theneural network is 0.014 ± stat ± sys events. C. Radiogenic background estimate
The radiogenic neutron background was also estimatedusing a GEANT4 simulation of the CDMS II tower con-figuration. Neutrons originating from the Th and U de-cay chains were simulated for each material in the setup(lead, polyethylene, copper). The primary energy spec-tra of the neutrons was generated using the SOURCES4Cpackage [74–76], which computes the neutron spectra dueto spontaneous fission and ( α, n ) reactions of alphas fromthe full decay chains within the matrix material. Theneutrons were propagated through all materials, eventu-ally creating NRs in the Ge and Si detectors. The single-and multiple-scatter rates were tabulated for each detec-tor to create high-statistics files labeled by their contami-nant source. The files were then weighted by the contam-ination level (see below) and normalized to the amount ofmaterial present to determine the final background rate.The source with the highest single-scatter contributionis
U in the copper cryostat enclosures (cans), followedclosely by U and Th contaminants in the lead shielding.The contamination levels were determined by a sep-arate γ simulation, again using GEANT4 and the sameCDMS II geometry. In this simulation, gammas from the Th and
U decay chains, K, and Co were gen-erated from inside the shielding components and tower9structures, and from radon daughters on surfaces. In or-der to reproduce the fiducial-volume cut already appliedto the data, the location of the energy deposition and theelectric field map in the detectors were used to deduce thefraction of charge collected on the inner and outer elec-trodes. The event was cut if it produced a measurablesignal on the outer electrode. A charge threshold wasimplemented in the simulation data by using the exper-imental inner-electrode charge thresholds (in keV ee ) andapplying them to the inner-electrode energy derived fromthe above procedure.To obtain the contamination estimate a χ minimiza-tion was performed to fit the 43 sources considered (seeFig. 12 for a schematic of the setup and Table V for abreakdown of source locations) to the ER data above15 keV ee , avoiding the ∼
10 keV ee Ge activation peaks.Some representative values for this contamination are4 mBq/kg for both
Th and
U in the inner polyethy-lene shield and 6–7 mBq/kg for
Th and
U in thetower 5 copper, the tower with the highest contamina-tion. The comparison between the ER spectrum and thesum of MC results using the best-fit contamination levelsis shown in Fig. 13.
FIG. 12. (Color online) Schematic of the CDMS II shield-ing configuration. The outermost vacuum can surface isshown near the center and radially outward the layers are:inner polyethylene (green), inner lead (light grey), outer lead(darker grey), outer polyethylene (green), and muon-vetoscintillator counters (blue) with attached light guides andphototubes (white and black). Cooling of the detectors isachieved through the “cold stem” from a dilution refrigera-tor (dark blue), while cabling passes through the “electronicsstem” on the other side of the setup.
Because of nearly degenerate fits to the spectrum,there is uncertainty in the final contamination values.The degree of uncertainty was estimated by bracket-ing the largest changes that produce acceptable fits inthe gamma MC and translating these into the resultingchanges in the neutron rate.In the 14 Ge detectors used for the WIMP-search anal- ysis, the expected raw radiogenic neutron single-scatterbackground rate was found to be (1.15 ± × − events/kg/day in the 10–100 keV energy range. Thestatistical uncertainty is negligible because of the largenumber of events simulated while the systematic uncer-tainty listed is due to the uncertainty in the contami-nation values and locations. This singles spectrum wasthen convolved with the final neutron cut efficiencies fromthe three timing analyses (see Fig. 10). The radiogenicspectrum-averaged efficiencies are: 35.3% for the clas-sic, 34.9% for the neural net, and 40.4% for the 5d- χ analysis. For the raw exposure of 612.2 kg days, the finalradiogenic neutron event backgrounds for the three anal-yses are: 0.025 ± stat ± sys events for the clas-sic and neural network, and 0.028 ± stat ± sys events for the 5d- χ analysis.For the extended analyses the estimates were madewith the approximate energy range 2–20 keV. The esti-mates were: 0.0131 ± stat ± sys for the clas-sic, 0.0148 ± stat ± sys for the neural net-work, and 0.0105 ± stat ± sys for the 5d- χ . D. Pb background estimate
Po decays via the α -decay process. This leaves an α particle of 5.4 MeV and a recoiling ( (cid:46)
100 keV)
Pb nu-cleus in the final state.
Po is in the
U decay chainand is part of possible post-
Rn plate-out contamina-tion. It is also the first post- Rn α -decayer after thelong-lived Pb. These facts make this decay likely tobe an important long-lived contamination, and becausea recoiling
Pb nucleus could have produced an energydeposition below ∼
100 keV with an ionization yield con-sistent with being an NR, it makes sense to evaluate thisbackground very seriously.Because of the CDMS II detector geometry (see Fig. 1)much of the surrounding surface of an interior detectoris comprised of another detector.
Po decay events inwhich the decay occurs in one detector and the
Pb re-coil is registered in an adjacent detector cannot contam-inate the signal region of the latter detector because ofthe clear 5.4 MeV α deposition in the former. Thereforethe most important component of the background comesfrom surfaces that are uninstrumented or those that areadjacent to uninstrumented surfaces. Another importantpoint is that since the alphas from this decay give sucha clear signature, and the decay is a two-body decay,an obvious way to estimate the number of unaccompa-nied (single) Pb events is to estimate the number ofunaccompanied α events from the decay. The angulardistribution is isotropic and any surface that can observea single α is approximately equally likely to observe asingle Pb recoil. This estimation of the number ofsingle
Pb events was carried out, and using rough es-timations for the passage fractions of these events whensubjected to the other analysis cuts, we expect approx-imately 0.187 signal-region
Pb events over the whole0
Contamination Outer Lead Inner Lead Inner Polyethylene Cans Innermost Can Surface Outermost Can Surface Towers Co (cid:88) - - (cid:88) (cid:88) (cid:88) (cid:88) K (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) Th (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) U (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) Rn - - - - - (cid:88) -TABLE V. Sources used for the radioactive contamination fitting procedure. The (cid:88) indicates that the respective source was usedin the fitting procedure for the location given. The cans are the nested copper cold stages of the cryostat, where surface sourceswere simulated on both the innermost and outermost surfaces. There are five separate sources for the copper components ofthe five detector towers. No simulated sources were placed in the outer polyethylene or the veto panels because they are behindtoo much shielding to yield a measurable amount of background. ionization energy [keV] - r a t e [ c oun t s ][ ke V k g - d ay ] − − − − −
10 1
FIG. 13. (Color online) Measured event-summed energy across all detectors (black solid line) compared to simulations ofcontaminants from common locations: lead shielding (blue diamonds), inner polyethylene shielding (magenta crosses), coppercans and tower components (cyan filled circles), and outermost can surfaces (brown stars). The sum of all simulated components(red inverted triangles) is also plotted and is in good agreement with the measured data. exposure.Since the
Pb recoil estimates were inferred from α counting, and the Pb recoils can come at differentenergies depending on how deeply the
Po parent isembedded into the originating surface, we do not havea good specification of the energy distribution of suchevents. Therefore, to be conservative we can use the sameestimate for the 10 keV and extended-threshold analy-ses. We expect 0.187 ± stat ± sys events overthe whole exposure; and have assigned a 100% system-atic uncertainty to account for the roughness of this esti-mate. This background estimate is clearly subdominantwith respect to the surface-event background estimates, but is larger than the cosmogenic or radiogenic neutronbackground estimates [62]. VIII. RESULTS
The background estimates for the primary 10 keV andextended-threshold analyses are summarized in Table VI.While the background estimates can be used to interpretthe overall results of the experiments they do not directlymodify the limit curves (see below).Results of direct WIMP-search experiments are usu-ally summarized as upper limits or signal contours in the1
Background 5d- χ (10 keV) Classic (extended)Leakage 1.19 ± > ± ± ± ± ± Pb 0.19 ± ± ± > ± WIMP-nucleon cross section versus WIMP-mass plane.Yellin’s optimum interval method [77] allows derivationof an upper limit on a signal rate in cases with unknownbackground. While the backgrounds in our signal regionare not completely unknown, this is a conservative ap-proach to setting upper limits on a possible signal. Thispresentation also requires assumptions about the WIMPdistribution in the galactic halo, the type of interactionbetween WIMPs and nucleons, and the nuclear form fac-tor for the interaction. The velocity distribution was as-sumed to be Maxwellian and was parametrized by therotational velocity at infinite radius and corrected for thefinite galactic escape velocity [49], which is taken to be544 km/s [78]. A WIMP mass density of 0.3 GeV/ c /cm was used for historical reasons, making the computed lim-its comparable to similar publications [49]. Some recentastrophysical measurements indicate different values [79];a correction for such deviations is a simple multiplica-tive factor and can be easily applied by the reader [80].A most probable WIMP velocity of 220 km/s was usedalong with a mean circular velocity of the Earth withrespect to the Galactic center of 232 km/s. The WIMPinteractions were assumed to be spin independent and theHelm form factor was used [49] for a natural Ge isotopicdistribution.The comparison of the present 10 keV 5d- χ and clas-sic extended results with other published limits and sig-nal contours is shown in Fig. 14. In the figure, our10 keV 5d- χ limit is combined with the CDMS II five-tower exposure acquired before July 2007, resulting ina limit that summarizes the full (and final) CDMS IIhigh-threshold sensitivity. The CDMS II/EDELWEISScombined limit [35] is also shown for comparison. Above ∼
100 GeV/ c WIMP mass the combined limit is compa-rable to our CDMS II combined result owing to the goodefficiency-averaged exposures of both of the experimentsin the relevant energy ranges.
A. Limit cross-checks
To gain insight into the effect of timing-cut figures ofmerit (efficiency, SAE, leakage) on WIMP-search results,we have constructed final limits for all of the timing-cutconstructions in this work. Although the SAE and ex-pected leakage were exclusively used to choose the pri- ] c WIMP mass [GeV/ ] [ c m S I σ -45 -44 -43 -42 -41 -40 -39 FIG. 14. This figure compares the main results fromthis analysis (the current 5d- χ analysis combined with priorCDMS II exposures, black solid and the classic extended-threshold analysis, black dashed) with previously publishedresults (all limits are at 90% C.L.): Darkside-50 [81] (orangetriple-dot-dashed) XENON100 [82] (blue triple-dot-dashed);LUX [83] (red dot-dashed); SuperCDMS low-threshold [84](green dotted); CDMS II/EDELWEISS combined [35] (pur-ple long dashed); CRESST-II [85] (magenta long-dot-dashed);CDMS II Si [39] (90% and 68% C.L. contours, blue, with thebest-fit point marked with a black dot); DAMA/LIBRA [43,86] (3 σ region, light red) and CoGeNT [40] (90% C.L.,brown). mary timing cuts, comparing all of the limits in this waygives cross-checks on how other parameters (efficiency,signal-region events) affect the reach of the experiment.The top part of Fig. 15 shows the limits derived usingthe optimum interval method for the different timing cutswith 10 keV thresholds. Relative to the original publica-tion, the ionization-based fiducial-volume and phonon-timing cuts have improved efficiencies in the analysisreported here, leading to the improved exposure andmore stringent limits. In terms of the overall spectrum-averaged detection efficiency for a 60 GeV/ c WIMP, theclassic timing-cut strategy (see Sec. V A) shows a 12% im-provement over the previously published version, abouthalf of which can be attributed to a reoptimization of theionization-based fiducial volume following the data repro-cessing. A similar improvement was seen for the neural-network timing-cut analysis (described in Sec. V B). Thelargest improvement of 29% in overall SAE efficiency wasachieved for the 5d- χ analysis (described in Sec. V C),owing primarily to an increased timing-cut efficiency inthe 15–90 keV energy range (see Sec. V E).For the 10 keV-threshold analyses, the 5d- χ sets themost stringent limit at a 60 GeV/ c WIMP mass, whilethe neural-network timing cut results in stronger limits2at and below 10 GeV/ c , an important region for furtherstudy [39]. The 5d- χ set weaker limits for low-massWIMPs because the cut was set to maximize sensitivityto a WIMP with mass 60 GeV/ c . This combined withthe fact that the 5d- χ method could set a very tight cutat low energies (it used an independent 10–20 keV binwhereas the other analyses were less granular) produceda poorer WIMP efficiency toward low recoil energies de-spite the lower expected background leakage at those en-ergies, but excellent sensitivity for high WIMP masses.The lower part of Fig. 15 shows the extended limitsin the low-WIMP-mass region. Each extended analysisconstrains the 8–10 GeV/ c mass region more stronglythan the higher-threshold analyses, and the classic tim-ing cut produces the strongest limit near the silicon-detector analysis best-fit point of M W = 8.6 GeV/ c and σ SI = 1.9 × − cm [39]. The extended limits are alsocompared to the previous low-threshold CDMS II results,which did not use a timing cut [34]. That analysis has alarger exposure toward lower recoil energies which ac-counts for the stronger limit set below a ∼ c WIMP mass. The classic analysis presented here hasa stronger limit by a factor of approximately 2.7 at aWIMP mass of ∼ c . IX. CONCLUSION
The reprocessed data did not produce significantchanges in the number of signal-region events, indicat-ing that uncertainties applied in the original processingof the CDMS II data set [33] were robust. All three setsof higher-threshold timing cuts produced similar limits,with small differences consistent with their correspond-ing exposure-optimization procedures. For example, the5d- χ analysis has a high efficiency at moderate recoilenergies (30–60 keV), but has a stringent timing cut atlower energies. It is well suited to provide the strongestlimits at high WIMP mass ( >
60 GeV/ c ), but will pro-duce fewer low-energy signal-region events. On the otherhand, the classic analysis at 10 keV threshold shows aslight weakening of the 90% C.L. limit for WIMP massesbelow about 18 GeV/ c , where sharp increases in thelimit curves indicate systematics near threshold. Theneural-network timing cut has been identified as a ro-bust method with the highest signal efficiency at low en-ergies (see Fig. 10) and good sensitivity at lower WIMPmass ( ≤ c ). The classic extended-threshold limitrules out about half of the silicon 68% C.L. region ob-tained in a previous CDMS II publication [39]. Thisindicates that the low-threshold results from the Ge de-tectors are marginally compatible with the Si-detectormeasurements taken during the same data period, understandard assumptions. Comparisons of such direct de-tection results on different nuclei will be a powerful toolfor understanding WIMP dark matter both in terms ofthe fundamental WIMP interactions [87, 88] and possible ] c WIMP mass [GeV/ ] [ c m S I σ -44 -43 -42 -41 -40 silicon 90% C.L. (2013)silicon 68% C.L. (2013)2010 result (this work) χ ] c WIMP mass [GeV/ ] [ c m S I σ -43 -42 -41 -40 -39 CDMS-II LT (2011) (this work) χ LT 5d-LT Classical (this work)LT Neural-Network (this work)
FIG. 15. (Color online) Experimental upper limits (90% confi-dence level) derived from each of the analyses presented in thiswork compared with the originally published [33] (black dot-ted) limits. The CDMS II Si contour is shown with the best-fitpoint marked with a black dot (WIMP mass of 8.6 GeV/ c and WIMP-nucleon cross section of 1.9 × − cm ) [39].(Top panel) The 10 keV threshold analyses. The 5d- χ limit(blue dashed) is the “primary” high-threshold result to bequoted from this work. The neural-network and classic limitsare shown as red dot-dashed and black solid lines respectively.(Bottom panel) The extended-threshold limits, focused on thelower WIMP-mass region. The same color code applies ex-cept that all of the analyses from this work correspond to theextended-threshold versions. The classic limit (black solid)is the “primary” extended-threshold result to be quoted fromthis work. The extended 5d- χ limit shown corresponds to thetiming-cut optimization assuming a 60 GeV/ c WIMP mass.For comparison, the previous CDMS II low-threshold limit isshown [34] (green triple-dot-dashed). backgrounds.3
ACKNOWLEDGMENTS
The CDMS Collaboration gratefully acknowledges thecontributions of numerous engineers and technicians; wewould like to especially thank Dennis Seitz, Jim Beaty,Bruce Hines, Larry Novak, Richard Schmitt and AstridTomada. In addition, we gratefully acknowledge assis-tance from the staff of the Soudan Underground Lab-oratory and the Minnesota Department of Natural Re-sources. This work is supported in part by the Na-tional Science Foundation, the US Department of Energy,NSERC Canada, and by MultiDark (Spanish MINECO).Fermilab is operated by the Fermi Research Alliance,LLC under Contract No. De-AC02-07CH11359. SLAC isoperated under Contract No. DE-AC02-76SF00515 withthe United States Department of Energy.4 [1] F. Zwicky, Helv. Phys. Acta , 110 (1933); Gen. Relativ.Gravit. , 207 (2008).[2] H. W. Babcock, Lick Observatory bulletin ; no. 498; LickObservatory bulletins ; no. 498. , 41 (1939).[3] Y. Sofue and V. Rubin, Annu. Rev. Astron. Astrophys. , 137 (2001).[4] K. G. Begeman et al. , Mon. Not. R. Astron Soc. ,523 (1991).[5] T. S. V. Albada et al. , Philos. Trans. R. Soc. London,Ser. A , 447 (1986).[6] B. J. Maughan et al. , Mon. Not. R. Astron. Soc. ,1193 (2004).[7] A. D. Lewis et al. , Astrophys. J. , 135 (2003).[8] J. Kneib et al. , Astrophys. J. , 643 (1996).[9] J. P. Kneib et al. , Astron. Astrophys. , 27 (1995).[10] E. S. Sheldon et al. , Astrophys. J. , 2217 (2009).[11] E. S. Sheldon et al. , Astrophys. J. , 2232 (2009).[12] D. E. Johnston et al. , (2007), arXiv:0709.1159.[13] J. D. Bekenstein, Phys. Rev. D , 083509 (2004).[14] M. Milgrom, Astrophys. J. , 365 (1983).[15] F. Kahlhoefer et al. , Mon. Not. R. Astron. Soc. , 2865(2014).[16] D. Clowe et al. , Astrophys. J. , L109 (2006).[17] G. W. Angus et al. , Astrophys. J. , L13 (2007).[18] G. W. Angus et al. , Mon. Not. R. Astron. Soc. , 138(2006).[19] M. Zemp, Mod. Phys. Lett. A , 2291 (2009).[20] J. Diemand and B. Moore, Adv. Sci. Lett. , 297 (2011-02-01T00:00:00).[21] M. Boylan-Kolchin et al. , Mon. Not. R. Astron. Soc. ,1150 (2009).[22] A. G. Riess et al. , Astrophys. J. , 665 (2004).[23] S. Perlmutter et al. , Astrophys. J. , 565 (1999).[24] A. G. Riess et al. , Astron. J. , 1009 (1998).[25] R. H. Cyburt et al. , Journal of Cosmology and Astropart.Phys. , 012 (2008).[26] D. J. Eisenstein et al. , Astrophys. J. , 560 (2005).[27] Ade, P. A. R. et al. , Astron. Astrophys. , A1 (2014).[28] G. Hinshaw et al. , Astrophys. J. Suppl. Ser. , 19(2013).[29] M. Kowalski et al. , Astrophys. J. , 749 (2008).[30] B. W. Lee and S. Weinberg, Phys. Rev. Lett. , 165(1977).[31] I. Aitchison, Supersymmetry in Particle Physics: AnElementary Introduction (Cambridge University Press,2007).[32] D. E. Kaplan, M. A. Luty, and K. M. Zurek, Phys. Rev.D , 115016 (2009).[33] Z. Ahmed et al. (CDMS Collaboration), Science ,1619 (2010).[34] Z. Ahmed et al. (CDMS Collaboration), Phys. Rev. Lett. , 131302 (2011).[35] Z. Ahmed et al. , Phys. Rev. D , 011102 (2011).[36] Z. Ahmed et al. (CDMS Collaboration), Phys. Rev. Lett. , 011301 (2009).[37] D. S. Akerib et al. (CDMS Collaboration), Phys. Rev. D , 052009 (2005).[38] Throughout this work the unit keV shall refer to true re-coil energy whereas keV ee are units of electron-equivalentrecoil energy. This means that an ionization yield y = E q /E r = 1 is assumed when the energy is calculated from the measured signal.[39] R. Agnese et al. (CDMS Collaboration), Phys. Rev. Lett. , 251301 (2013).[40] C. E. Aalseth et al. (CoGeNT Collaboration), Phys. Rev.D , 012002 (2013).[41] G. Angloher et al. , Eur. Phys. J. C , 1971 (2012).[42] C. E. Aalseth et al. (CoGeNT Collaboration), Phys. Rev.Lett. , 131301 (2011).[43] R. Bernabei et al. , Eur. Phys. J. C , 39 (2010).[44] J. Hellmig et al. , Nucl. Instrum. Methods Sec. A ,308 (2000).[45] D. Akerib et al. , Nucl. Instrum. Methods Sec. A , 476(2008).[46] B. Neganov and V. Trofimov, Otkryt. Izobret. , 215(1985).[47] P. N. Luke, J. Appl. Phys. , 6858 (1988).[48] R. C. Alig and S. Bloom, Phys. Rev. Lett. , 1522(1975).[49] J. Lewin and P. Smith, Astropart. Phys. , 87 (1996).[50] S. Butterworth, Exper. Wireless and the Wireless Eng. , 536 (1930).[51] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling, Numerical Recipes in C: The Art of Scientific Computing (Cambridge University Press, 1992).[52] S. Golwala,
Exclusion Limits on WIMP-Nucleon ElasticScattering Cross-Section from the Cryogenic Dark MatterSearch , Ph.D. thesis, University of California, Berkeley(2000).[53] V. Mandic et al. , in
AIP Conference Proceedings , Vol.605 (AIP, 2002) pp. 509–512.[54] Z. Ahmed,
A Dark-Matter Search Using the Final CDMSII Dataset and a Novel Detector of Surface Radiocontam-ination , Ph.D. thesis, California Institute of Technology(2012).[55] A. Abramov et al. , Nucl. Instrum. Methods Sec. A ,209 (2002).[56] We reject events that have phonon pulse multiplicities atleast 4 larger than the charge pulse multiplicity. Further,if the phonon pulse multiplicity is 3 we reject events withcharge multiplicity 0.[57] J. Lindhard et al. , Mat. Fys. Medd. Dan. Vid. Selsk. ,10 (1963).[58] Specifically the form σ y ( E r ≤ ¯ E ) = c · E dr and σ y ( E r > ¯ E ) = e was used, where c , d , e , and ¯ E (the energy thresh-old) are fitted parameters.[59] D. S. Akerib et al. (CDMS Collaboration), Phys. Rev.Lett. , 011302 (2006).[60] J. Zhang, A Dark Matter Search Using the Final CDMS-II Data and 100 mm SuperCDMS Germanium DetectorIonization Test , Ph.D. thesis, University of Minnesota(2014).[61] Note the ERs and NRs have less clear “surface discrimi-nation” within the respective populations. The differencebetween ERs and NRs arises because Luke phonons pro-duce a fast-rising component due to their ballistic natureand consistent near-surface production component. ERsgive a greater fraction of Luke phonons, so will generi-cally be faster.[62] T. J. Hofer,
Development of CDMS-II Surface Event Re-jection Techniques and Their Extensions to Lower EnergyThresholds , Ph.D. thesis, University of Minnesota (2014). [63] K. Pearson, Philosophical Magazine Series 6 , 559(1901).[64] I. T. Nabney, NETLAB: algorithms for pattern recogni-tions (Springer, London; New York, 2002).[65] J. M. Kiveni,
A Search for WIMP Dark Matter using anOptimized Chi-square Technique on the Final Data fromthe Cryogenic Dark Matter Search Experiment (CDMSII) , Ph.D. thesis, Syracuse University (2012).[66] J. Filippini,
A Search for WIMP Dark Matter Usingthe First Five-Tower Run of the Cryogenic Dark MatterSearch , Ph.D. thesis, University of California, Berkeley(2008).[67] J. Sander,
Results from the Crogenic Dark Matter SearchUsing a Chi Squared Analysis , Ph.D. thesis, Universityof California, Santa Barbara (2007).[68] J. R. Klein and A. Roodman, Annu. Rev. Nucl. Part. Sci. , 141 (2005).[69] G. Cowan et al. (Particle Data Group), Phys. Rev. D ,010001 (2012).[70] I. Ruchlin and R. W. Schnee, Nucl. Instrum. MethodsSec. A , 336 (2012).[71] R. Agnese et al. (SuperCDMS Collaboration), Inprogress.[72] V. Kudryavtsev, Comp. Phys. Comm. , 339 (2009).[73] S. Kasahara, A Study of Cosmic Ray Composition in theKnee Region using Multiple Muon Events , Ph.D. thesis,University of Minnesota (1997). [74] W. B. Wilson et al. , in
Proceedings of the 12th Bi-ennial Topical Meeting of the American Nuclear Soci-ety/Radiation Protection and Shielding Division, SantaFe, NM, April 14-18, 2002 (2002).[75] W. Wilson et al. , Radiation Protection Dosimetry ,117 (2005).[76] W. Charlton et al. , Nucl. Instrum. Methods Sec. B ,1 (1998).[77] S. Yellin, Phys. Rev. D , 032005 (2002).[78] M. C. Smith et al. , Mon. Not. R. Astron. Soc. , 755(2007).[79] J. Lavalle and S. Magni, Phys. Rev. D , 023510 (2015).[80] R. Bernabei et al. , Phys. Lett. B , 757 (1996).[81] P. Agnes et al. , Phys. Lett. B , 456 (2015).[82] E. Aprile et al. , Phys. Rev. Lett. , 181301 (2012).[83] D. Akerib et al. (LUX Collaboration), Phys. Rev. Lett. , 091303 (2014).[84] R. Agnese et al. (SuperCDMS Collaboration), Phys. Rev.Lett. , 241302 (2014).[85] G. Angloher et al. , Eur. Phys. J. C , 3184 (2014),10.1140/epjc/s10052-014-3184-9.[86] C. Savage, G. Gelmini, P. Gondolo, and K. Freese, J.Cosmol. Astropart. P. , 010 (2009).[87] N. Anand, A. L. Fitzpatrick, and W. C. Haxton, Phys.Rev. C , 065501 (2014).[88] K. Schneck et al. (SuperCDMS Collaboration), Phys.Rev. D91