First Search for Bosonic Superweakly Interacting Massive Particles with Masses up to 1 MeV = c 2 with GERDA

We present the first search for bosonic superweakly interacting massive particles (super-WIMPs) as keV-scale dark matter candidates performed with the GERDA experiment. GERDA is a neutrinoless double- β decay experiment which operates high-purity germanium detectors enriched in 76 Ge in an ultralow background environment at the Laboratori Nazionali del Gran Sasso (LNGS) of INFN in Italy. Searches were performed for pseudoscalar and vector particles in the mass region from 60 keV =c 2 to 1 MeV =c 2 . Further distribution of work must maintain attribution to the author(s) published No evidence for a dark matter signal was observed, and the most stringent constraints on the couplings of super-WIMPs with masses above 120 keV =c 2 have been set. As an example, at a mass of 150 keV =c 2 the most stringent direct limits on the dimensionless couplings of axionlike particles and dark photons to electrons of g ae < 3 × 10 − 12 and α 0 = α < 6 . 5 × 10 − 24 at 90% credible interval, respectively, were obtained.

No evidence for a dark matter signal was observed, and the most stringent constraints on the couplings of super-WIMPs with masses above 120 keV=c 2 have been set. As an example, at a mass of 150 keV=c 2 the most stringent direct limits on the dimensionless couplings of axionlike particles and dark photons to electrons of g ae < 3 × 10 −12 and α 0 =α < 6.5 × 10 −24 at 90% credible interval, respectively, were obtained. DOI: 10.1103/PhysRevLett.125.011801 The evidence for the existence of nonbaryonic dark matter (DM) in our Universe is overwhelming. In particular, recent measurements of temperature fluctuations in the cosmic microwave background radiation yield a 26.4% contribution of DM to the overall energy density in the ΛCDM model [1]. However, all evidence is gravitational in nature, and the composition of this invisible form of matter is not known. Theoretical models for particle DM yield candidates with a wide range of masses and scattering cross sections with standard model particles [2][3][4]. Among these, so-called bosonic superweakly interacting massive particles (super-WIMPs) with masses at the keV-scale and ultraweak couplings to the standard model can be cosmologically viable and produce the required relic abundance [5,6].
Direct DM detection experiments, as well as experiments built to observe neutrinoless double-β (0νββ) decay, can search for pseudoscalar [also known as axionlike particles (ALPs)] and vector (also known as dark photons) super-WIMPs via their absorption in detector materials in processes analogous to the photoelectric effect (also known as axioelectric effect in the case of axions). The ALP and dark photon energy is transferred to an electron, which deposits its energy in the detector. The expected signature is a full absorption peak in the energy spectrum, corresponding to the mass of the particle, given that these DM candidates have very small kinetic energies [7]. For ALPs the coupling to electrons is parametrized via the dimensionless coupling constant g ae [8], while for dark photons a kinetic mixing α 0 with strength κ [9] is introduced in analogy to the electromagnetic fine structure constant α, such that α 0 ¼ ðeκÞ 2 =4π.
Here we describe a search for super-WIMP absorption in the germanium detectors operated by the GERDA Collaboration, extending for the first time the mass region to 1 MeV=c 2 . At masses larger than twice the electron mass, vector particles can decay into ðe þ ; e − Þ pairs and their lifetime would be too short to account for the DM [6].
The primary goal of GERDA is to search for the 0νββ decay of 76 Ge, deploying high-purity germanium (HPGe) detectors enriched up to 87% in 76 Ge. The experiment is located underground at the Laboratori Nazionali del Gran Sasso (LNGS) of INFN, Italy, at a depth of about 3500 m water equivalent. The HPGe detector array is made of 7 enriched coaxial and 30 broad energy germanium (BEGe) diodes, with average masses of 2.2 kg and 667 g, respectively, leading to a higher full absorption efficiency for the larger coaxial detectors. It is operated inside a 64 m 3 liquid argon (LAr) cryostat, which provides cooling and a highpurity, active shield against background radiation. The cryostat is inside a water tank instrumented with photomultiplier tubes (PMTs) to detect Cherenkov light from muons passing through, and thus reduces the muon-induced background to negligible levels. A detailed description of the experiment can be found in Ref. [19], while the most recent 0νββ decay results are presented in Ref. [20]. Because of its ultralow background level [21] and excellent energy resolution [∼3.6 and ∼3.0 keV full width at half maximum (FWHM) for coaxial and BEGe detectors at Q ββ ¼ 2039 keV, respectively], the GERDA experiment is well suited to search for other rare interactions, in particular for peaklike signatures as expected from bosonic super-WIMPs. Here we make the assumption that super-WIMPs constitute all of the DM in our Galaxy, with a local density of 0.3 GeV=cm 3 [22]. The absorption rate for dark photons and ALPs in an Earth-bound detector can be expressed as [5] and respectively. Here g ae and α 0 =α are the dimensionless coupling constants, A is the atomic mass of the absorber, σ pe is the photoelectric cross section on the target material (germanium), and m v and m a are the DM particle masses. The linear versus inverse proportionality of the rate with the particle mass is due to the fact that rates scale as flux times cross section, where the cross section is proportional to m 2 a and α 0 =α in the pseudoscalar and vector boson case, respectively [5].
We perform the search for super-WIMPs in the (200 keV=c 2 -1 MeV=c 2 ) mass range on data collected between December 2015 and April 2018, corresponding to PHYSICAL REVIEW LETTERS 125, 011801 (2020) 011801-2 58.9 kg yr of exposure. The energy threshold of the HPGe detectors was lowered in October 2017 and enabled a search in the additional mass range of (60-200 keV=c 2 ), corresponding to 14.6 kg yr of exposure accumulated until April 2018. The individual exposures for BEGe and coaxial detectors above (below) 200 keV are 30.8 (7.7) and 28.1 (6.9) kg yr, respectively. The lower energy bounds for our analysis were motivated by the energy thresholds of the Ge detectors and the shape of the background spectrum (dominated by 39 Ar decays) and the size of the fit window, as explained in the following.
In GERDA the energy reconstruction of events is performed through digital pulse processing [23]. Events of nonphysical origin such as discharges are rejected by a set of selection criteria based on waveform parameters (i.e., baseline, leading edge, and decay tail). The efficiency of these cuts for accepting signal events was estimated at > 98.7%. Since super-WIMPs would interact only once in a HPGe diode, events tagged in coincidence with the muon or LAr vetos, or observed in more than one germanium detector, were rejected as due to background interactions. We use the same set of cuts as in the GERDA main analysis for the 0νββ decay [20], with the exception of the pulse shape discrimination cut, which had been tailored to the high-energy 0νββ decay search. The muon and LAr veto accept signal events with efficiencies of 99.9% [24] and 97.5% [20], respectively.
The total efficiency to observe a super-WIMP absorption in the HPGe diodes was determined as where the efficiency of the event selection criteria ϵ cuts and the exposure E of each dataset were taken into account. The index i runs over the individual detectors of that dataset, containing N det detectors, E i is the exposure, f av;i the active mass fraction, and ϵ fep;i the efficiency for detection of the full-energy absorption of an electron emitted in the interaction. With the exception of ϵ fep;i all parameters were identical to those in the analysis presented in Ref. [20]. The full-energy absorption efficiency ϵ fep;i accounts for partial energy losses, for example, in a detector's dead layer. This efficiency was estimated for each detector at energies between 60 and 1000 keV with a Monte Carlo simulation of uniformly distributed electrons in the active volume of the detector using the MaGe framework [25]. Table I shows the average full-energy absorption detection efficiencies ϵ fep and the total efficiencies ϵ tot at the lower and upper boundaries of the search region. At 60 keV, the full-energy absorption was estimated as 99.5% for all detectors, while at 1000 keV it is 95.1% and 96.2% on average for BEGe and coaxial detectors, respectively. The energy dependence of the efficiency is caused by the photoabsorption cross section and the different size of the germanium diodes. The events which survived all selection criteria (with total efficiencies between 85.7% and 81.4%; see Table I) are shown in Fig. 1 for the coaxial and BEGe detector datasets.
The expected signal from super-WIMPs has been modeled with a Gaussian peak broadened by the energy resolution of the HPGe detectors. To estimate the potential signals from these particles we performed a binned Bayesian fit (with a 1 keV binning, while the systematic uncertainties on the energy scale are estimated at 0.2 keV) FIG. 1. The energy spectra of the BEGe and coaxial datasets, normalized by exposure. Only events with energies up to 1 MeV were considered in the analysis. The coaxial dataset shows a significantly higher event rate (mainly from 39 Ar decays) at energies below 500 keV due to the larger surface area of the signal read-out electrodes [20]. The dashed lines indicate the positions of the main known background γ lines, also listed in Table II. of the signal and a background model of the data. The fit was performed within a window of 24 keV in width, centered on the energy corresponding to the hypothetical mass of the particle and sliding with 1 keV step to examine each mass value. The total number of counts from signal and background was determined as follows: where the Gaussian function G 0 models the peak signal of super-WIMPs at a fixed energy E 0 , corresponding to their mass. The Gaussian G γ models the background γ lines with energy E γ listed in Table II in case it is found within the sliding fit window. For more than one background γ line, Eq. (4) is modified accordingly to model all the peaks. N 0 and N γ are the counts in the fitted signal and background peaks, respectively. The effective energy resolutions σ 0 and σ γ of the detectors from the combined spectra are fixed to the values obtained from the regularly acquired calibration data, with systematic uncertainties around 0.1 keV [20]. Finally, the polynomial fit function FðEÞ describes the continuous background, and was chosen as a first-and second-order polynomial for energies above and below 120 keV, respectively. The higher-order polynomial at lower energies is motivated by the curvature of the 39 Ar β spectrum; see Fig. 1. At other energies, the spectrum has an approximately linear shape, and thus a first-order polynomial was judged sufficient. The Bayesian fit was performed with the Bayesian analysis toolkit BAT [26] using the Markov chain Monte Carlo technique [27] to compute the marginalized posterior probability density function (PDF) given energy values of the data E, PðR S jEÞ, where R S is the signal rate, i.e., the number of counts normalized by exposure.
The probability for the signal count rate PðR S ; θjE; MÞ, given data E and a model M, is described by Bayes' theorem as PðR S ;θjE;MÞ ¼ PðEjR S ;θ;MÞπðR S ÞπðθÞ R R PðEjR S ;θ;MÞπðR S ÞπðθÞdθdR S : The denominator defines the overall probability of obtaining the observed data given a hypothetical signal. The numerator includes prior probabilities π for the signal count rate R S and for the nuisance parameters θ (e.g., background shape) estimated before performing the fit. For θ, flat priors were adopted, bound generously according to a preliminary fit with the MINUIT algorithm [28]. For the signal count rate R S the uniform (i.e., constant over the defined range) prior probability was constructed to be positive, with the upper bound defined by the total number of events in the signal region plus 10 times the expected Poisson fluctuations. The conditional probability PðEjR S ; θ; MÞ is estimated according to the super-WIMP interaction model given by Eq. (4) and Poisson fluctuations in the data. The reported results were obtained from the combined fit of the BEGe and coaxial datasets. First, the BEGe dataset was fit using a flat prior for the signal count rate R S , and the obtained posterior was used as a prior for the fit of the coaxial dataset. The results of the latter were then employed to evaluate the corresponding coupling constants of the super-WIMPs. The detection of the signal is ruled out when the significance of the best fit value for the count rate is less than 5σ, estimated as half of the 68% quantile of the posterior PDF. Additionally, if a fitted signal is in close proximity (within 5σ of the energy resolution) to a known background γ line, an upper limit was set irrespective of the mode, as uncertainties in the background rate do not allow us to reliably claim an excess signal above γ lines. An example for the fit using the model described by Eq. (4) for two different background functions FðEÞ is shown in Fig. 2.
The obtained posterior PDFs do not show evidence for a signal in the energy range of the analysis. We thus set 90% credible interval (C.I.) upper limits on the signal count rate, corresponding to the 90% quantile of the posterior PDF PðR S jEÞ, accounting for the detection efficiencies according to Eq. (3).
The 90% C.I. limits on the signal rate R S were converted into upper limits on the coupling strengths using Eqs. (1) and (2). The results are presented in Fig. 3. We compare these to direct detection limits from CDEX [17], EDELWEISS-III [16], LUX [12], the Majorana Demonstrator [14], PandaX-II [11], SuperCDMS [15], XENON100 [13], and XMASS [10], as well as to indirect limits from horizontal branch and red giant stars [5]. Above 120 keV=c 2 indirect α 0 =α limits from decays of vectorlike particles into three photons (V → 3γ) are significantly lower (ranging from 10 −12 at masses of 100 keV=c 2 to 10 −16 at 700 keV=c 2 ) than the available direct limits (not shown) [6]. The improvement in sensitivity with respect to other crystal-based experiments is due to the much larger exposure in GERDA and the lower background rate over all of the search region.
The weakening of our upper limits with increasing mass is primarily due to the steep decrease of the photoelectric cross section from about 45 b at 100 keV to 0.085 b at 1 MeV that overrules both the linear and inverse mass dependence in Eqs. (1) and (2). The fluctuations in the upper limit curves are due to background fluctuations, where prominent peaks come from known γ lines, shown in Table II.
To summarize, in this Letter we demonstrated the capability of GERDA to search for other rare events besides the 0νββ decay of 76 Ge. We performed a search for keVscale DM in the form of bosonic super-WIMPs based on data with exposures of 58.9 and 14.6 kg yr in the mass ranges of 200 keV=c 2 -1 MeV=c 2 and 60-200 keV=c 2 , respectively. Upper limits on the coupling strengths g ae and α 0 =α were obtained from a Bayesian fit of a background model and a potential peaklike signal to the measured data. Our limit is compatible with other direct searches in the mass range (60-120 keV=c 2 ) where the strongest limits were obtained by xenon-based DM experiments due to higher exposures and lower background rates in this low-energy region. Our search probes for the first time the mass region up to 1 MeV=c 2 and sets the best direct constraints on the couplings of super-WIMPs over a large mass range from (120 keV=c 2 -1 MeV=c 2 ). As an example, at a mass of 150 keV=c 2 , the most stringent direct limits on the dimensionless couplings of axionlike particles and dark photons to electrons of g ae <3×10 −12 and α 0 =α< 6.5×10 −24 (at 90% C.I.), respectively, were established. The limits are affected by the known background γ lines, listed in Table II, due to higher background rate at these energies.
The sensitivity to new physics is expected to improve in the near future with the upcoming LEGEND-200 experiment. The experimental program aims to decrease the background rate and increase the number of HPGe detectors operated in an upgraded GERDA infrastructure at LNGS [29].
FIG. 2. Best fit (red lines) and 68% uncertainty band (yellow bands) from marginalized posterior PDFs of the model parameters assuming a hypothetical signal at E 0 in the BEGe dataset. Top: Fit of a signal assumed at E 0 ¼ 520 keV; the excess is at a level of 2.6σ. A first-order polynomial is used for the continuous background and a Gaussian for the background γ line due to the decay of 85 Kr. Bottom: Fit of a signal assumed at E 0 ¼ 87 keV, using a second-order polynomial for the continuous background. Only part of the data was acquired with a lower energy threshold, resulting in a lower exposure for data below 200 keV=c 2 and causing the steplike feature around this energy. Results from other experiments (see text) are also shown, together with indirect constraints from anomalous energy losses in horizontal branch (HB) and red giant (RG) stars (we refer to Ref. [6] for details).