A Search for Light Fermionic Dark Matter Absorption on Electrons in PandaX-4T

We report a search on a sub-MeV fermionic dark matter absorbed by electrons with an outgoing active neutrino using the 0.63 tonne-year exposure collected by PandaX-4T liquid xenon experiment. No significant signals are observed over the expected background. The data are interpreted into limits to the effective couplings between such dark matter and electrons. For axial-vector or vector interactions, our sensitivity is competitive in comparison to existing astrophysical bounds on the decay of such dark matter into photon final states. In particular, we present the first direct detection limits for an axial-vector (vector) interaction which are the strongest in the mass range from 25 to 45 (35 to 50) keV/c$^2$.


INTRODUCTION
The nature of dark matter (DM), which makes up around 85% of the total mass in the Universe [1], is a top scientific mystery. The cold dark matter featured in the standard cosmological model (ΛCDM) is consistent with large scale cosmological observations at different eras of the cosmic evolution. Leading cold dark matter candi-dates include the Weakly Interacting Massive Particles (WIMPs) and axions, which have been searched with tailored experiments for decades [2][3][4][5][6][7][8][9][10][11][12][13]. Warm dark matter, on the other hand, is also a viable dark matter solution, which mitigates the so-called small scale problems of ΛCDM [14][15][16][17][18] with a larger free streaming distance [19][20][21], while maintaining excellent agreement with observations at the large scale. A representative warm dark matter candidate is the keV-scale sterile neutrino [22], which received particular attention recently due to the 3.5 keV unidentified excess in the x-ray spectrum from satellites [23][24][25]. However, comprehensive analyses later have challenged the sterile neutrino interpretations of the excess [26,27].
Nevertheless, a more general light neutral fermionic particle can still be a warm or cold DM candidate, depending on the initial conditions in the thermal history. If such fermionic particles (noted as χ hereafter) couple to neutrinos or charged leptons via an effective interaction, they can be probed experimentally [28][29][30]. In a DM direct detection experiment, the tree-level process of χ absorbed on an electron with an outgoing active neutrino (χe → eν) can be sensitively searched, as the electron obtains a kinetic recoil through the mass of χ [29]. Such a search becomes particularly attractive for axialvector (A) and vector (V) interactions, where the satellite x-ray limits are relatively weaker since χ → νγ can only be produced by two-loop diagrams due to the gauge symmetry and quantum electrodynamics charge conjugation symmetry [30]. In fact, the recent electron recoil excess observed in XENON1T was interpreted as a fermionic DM with a mass of 60 keV/c 2 absorbed on electrons via a vector interaction [29][30][31]. However, such interpretation has some tension with cosmological observations as the invisible decay of χ → 3ν would have injected more radi-ation into the early universe, leaving traces in the Hubble constant and matter power spectrum [30]. For these DM candidates, constraints also arise from the relic density to avoid overproduction through the "freeze-in" mechanism [32]. In this letter, we present the first dedicated search on fermionic DM-electron absorption using recent data from PandaX-4T experiment to further elucidate the situation.
For the effective interaction, the V and A operators can be written as [30] where the standard model (SM) left-handed neutrino is taken, and 1/Λ 2 is the Wilson coefficient with a dimension of inverse mass square. Note that our experimental constraint stays the same if a light sterile neutrino (m ν m χ ) replaces the active neutrino in Eq. 1, and the changes in other astrophysical and cosmological constraints are also negligible.
Given the tiny mass of the active neutrino, the absorption signal has distinguishable peak features in the electron energy spectrum. Neglecting the halo velocity of the DM, we have where m χ is the static mass of χ, E nl (E nl > 0) is the binding energy of the initial electron on the state |nl , |q| is the outgoing neutrino energy and E R is the recoil energy of the electron. If the binding energy is omitted, E R = m 2 χ /2(m χ + m e ). The expected event spectrum of the visible energy E vis summing up electron recoiling and deexcitation photon energy is where N T is the number of electron targets per unit mass, m e is the electron mass, and ρ χ ∼ 0.3 GeV/cm 3 /c 2 is the local DM density [33,34]. K nl is the atomic K factor [35] also known as the ionization form factor [29,36]. |M (q)| 2 is the scattering amplitude with the leading term for O V eνχ and O A eνχ , respectively. We define σ e ≡ m 2 χ /4πΛ 4 is the free-electron-χ total cross section for the vector operator when m χ m e . Two examples of the event spectra with m χ = 50, 130 keV/c 2 are shown in This work utilizes the data from the commissioning run (Run 1) of PandaX-4T [2], which is located in the B2 Hall of China Jinping Underground Laboratory [37,38]. The dual-phase liquid xenon time projection chamber (TPC) contains 3.7 tonne xenon in the sensitive region, viewed by two arrays of three-inch photomultipler tubes (PMTs) from the top and the bottom. A recoiling event generates the prompt scintillation photons (S1) and the delayed electroluminenscence photons (S2). We use the same dataset as in the search for spin- independent WIMP-nucleus scattering [2], with a fiducial volume (FV) 2.67 ± 0.05 tonne xenon and a 0.63 tonneyear exposure. The Run 1 data are divided into 5 subsets with slightly different field configurations and background levels. For the range of the electron energy considered in this Letter (< 30 keV), the relativistic correction to the ionization form factor is less than 10%, posing a minor effect on our search results. The electron-equivalent energy of each event is reconstructed as where 0.0137 keV is the average work function of liquid xenon [39], PDE is the photon detection efficiency, EEE is the electron extraction efficiency, and SEG b is the single electron gain, using S2 b collected by the bottom PMT array. G2 b is conventionally defined as EEE × SEG b for the correction in S2 b . These parameters are predetermined using monoenergetic electronic recoil (ER) peaks at the energy range from 41.5 to 408 keV [2]. Details of the ER signal response model will be discussed elsewhere [40], and we will only outline steps relevant to this Letter. The detector ER response is calibrated by the 220 Rn injection data (see Ref. [41] for the method). The model based on the noble element simulation techniques (NESTv2.2.1) [42,43] is then used to fit the data below 30 keV, shown in Fig. 2. The so-called ionization recombination coefficient r, an intermediate random variable in the simulation, follows a Gaussian distribution with the median (µ recomb ) and the fluctuation (σ recomb ). In the fit, µ recomb and σ recomb are allowed to vary from nominal values from NEST by a quadratic function in E vis , and a scaling factor (p f ), respectively. Other detector parameters, PDE, G2 b , and SEG b , are treated as nuisance parameters. After the fit, the quadratic parameters on µ recomb are linearly cast into three orthogonal variables via a principle component analysis, with the uncertainty dominated by the major parameter p 0 along the "long axis" (30% uncertainty). The fitted model agrees well with the distribution of calibration data. The complete set of fitted detector parameters p * are summarized in Table I.
Because of the peak feature, understanding the detector resolution of ER events is critical for this search. Strong validation comes from the overall agreement between the data and model quantiles along each equienergy line in Fig. 2. Further validation comes from the measured 1σ fractional energy resolution of 83m Kr data (41.5 keV and 6.8%), which agrees with the prediction of our model (7.0%) [44]. Approximately, our model leads to a 1σ resolution in E rec as 0.073+0.173E vis −6.5× where E vis is in unit of keV, and the smeared spectrum should be further corrected by the efficiency function in Fig. 1. In our later analysis, however, the fit is performed in S1 and S2 b using the full model, including uncertainties in p * . Rn calibration data in log 10 (S2 b ) vs. S1, taken under a drift field of 93 V/cm, and the comparison of 10% and 90% quantiles between the data (blue) and the best fit model (green). The average 1σ variation of p * in the quantiles is 9%. The dashed gray lines are the equi-ER-energy lines.
The main backgrounds in the light fermionic DM search follow the same compositions as in the WIMP search [2], summarized in Table. I. The ER background components include tritium, flat ER (including radon, 85 Kr, material background, and solar neutrino), 127 Xe Lshell electron capture (5.2 keV), and 136 Xe two-neutrino beta decay. The accidental, surface, neutron, and coherent solar neutrino-nucleus backgrounds are much less important in this ER analysis. , and best fits for the detector response parameters p * (upper), and signal and individual background components (lower). Parameters p * include PDE, G2 b , SEG b and ionization recombination parameters p 0 and p f (see text for details). Similar to Ref. [2], the common DM signal and tritium background for each set are floating in the fit. The best fit values of the number of events have been corrected for their efficiencies. The change of the tritium rate between set 4 and 5 is attributed to the removal effects from the hot purifier [45]. The variation in radon rate is primarily introduced from the krypton distillation tower. The candidate signals are the same 1058 events as in Ref. [2]. To carry out the search, the data are fitted with a profile likelihood ratio analysis [46]. The likelihood function for a given signal value µ, evaluated in S1 and S2 b , is constructed as where L n is the likelihood function for dataset n, and the common Gaussian penalty terms Gs are related to uncertainties of the background and response parameters, with corresponding nuisance parameters (1σ uncertainties) δ b and δ p * (σ b and σ p * ), respectively (Table. I). The measured number of events for each set N n meas is compared to the expected event number N n fit with a Poisson distribution. In each dataset, N n fit is a sum of fitted signal (N n µ n µ ) and background event numbers (N b n µ ), with corresponding probability density distribution function (PDF) P n b and P n µ in S1 and S2 b . The test static, q µ , is computed as the difference of −2 log L tot between the test signal value µ and the best fit [46]. A cross-check of the fit performed in one-dimensional reconstructed energy is made, and the results are consistent.
In the fit, signal and background PDFs need to be up-dated for each set of response parameters p * , which is computational expensive. The overall efficiencies n b,µ are also functions of p * for different event components. To include the uncertainties in p * , we utilize a Monte Carlo reweighting technique [47][48][49][50] that a pool with millions of simulation events generated and saved ahead of time together with intermediate random variables drawn from p * . When p * change to a different set of values, a new weight for the saved event k is computed as the product of probability density ratios of all related intermediate variables. For example, if parameters (p 0 , p f ) are updated to (p 0 , p f ), the weight factor for the ionization recombination coefficient r k for the event k is computed as  Table. I (see legend). The gray histogram is the sum of the best fit background. The dotted pink histogram is the best fit of XENON1T's low energy ER excess (m χ = 59 keV/c 2 and Λ = 0.98 TeV) from Ref. [30]. Bottom: the relative deviation of the data in the background-only model fit compared to the statistical uncertainty (green: ±1σ, yellow: ±2σ) in each 1 keV energy bin.
The search of DM was carried out for m χ larger than 10 keV/c 2 and smaller than 180 keV/c 2 , well covered by the Run 1 data release [2]. The global best fit is found at m χ = 130 keV/c 2 , with all key components projected to E rec (Eq. 5) in Fig. 3, together with the best fit signal to XENON1T's ER excess at m χ = 60 keV/c 2 (Ref. [30]) for comparison. The computed χ 2 using statistical uncertainties is 27.7 for 30 bins, indicating excellent agreement between the data and fit. The best fit detector and background nuisance parameters are all consistent with 0 within 2σ (Table I). For the best fit signal strength, the local significance at m χ = 130 keV/c 2 is 1.7σ relative to zero, based on the p value of q 0 (data tested with background-only hypothesis) evaluated using pesudodata generated with background-only simulations. The significance reduces to 0.6σ with the look-elsewhere effect considered [34,51]. Note that the signal peak shape is rather smeared due to atomic effects and detector resolution. For example, at m χ = 130 keV/c 2 , the full width at half maximum is 5.2 keV. The slight and gentle excess between 9 to 14 keV is picked up by the fit, but much less so for narrower features (or fluctuations) in the spectrum, for example, at 4 keV. The ∼3σ data deficit between 20 and 21!keV was investigated and concluded to be most likely a background downward fluctuation.
With no significant excess identified, the 90% confidence level (CL) upper limits on the χe → eν for the axial-vector and vector operators are set with the standard profile likelihood ratio analysis procedure, with background downward fluctuation power-constrained to −1σ [34,52]. The limit on axial-vector interaction is three-fold stronger than the vector interaction (Eq. 4). The behavior of the sensitivity as a function of m χ (blue line in Fig. 4) is driven by several factors. According to Eqs. 3 and 4, for a constant σ e v χ , the scattering rate between DM and free electrons scales with 1/m 3 χ . In addition, both the detector threshold and background level contribute. For m χ less than 20 keV/c 2 , the sensitivity is worsened due to the threshold, whereas in the region where m χ > 90 keV/c 2 , the sensitivity flattens due to the reduction of background rate for E rec 6 keV (Fig. 3). Note that for m χ > 100 keV/c 2 , the ionization form factor of the nonrelativistic atomic response contains non-negligible uncertainty, so we plot the limits in dashed lines [53]. The local upward fluctuation at m χ = 130 keV/c 2 , and the downward structures at m χ = 90 and 160 keV/c 2 , correspond to local data fluctuations around 13, 7, and 20 keV in Fig. 3, respectively. For comparison, the astrophysical constraint from the leading visible decay channel (χ → γγν) for O A eνχ from Ref. [30] is overlaid on Fig. 4 (blue line) combining results from Insight-HXMT [54], NuSTAR/M31 [55] and INTEGRAL/08 [56], together with the cosmological constraint (dashed magenta) and the freeze-in overproduction constraint (orange) [30]. Our dark matter direct detection data have produced leading limits in the mass region 35 to 55 (25 to 45) keV/c 2 for the vector (axialvector) operator. The 1σ amd 2σ contours obtained by fitting XENON1T data from Ref. [30] are also overlaid.
Our limit is touching the center region of XENON1T's best fit. The 90% CL exclusion limits (red lines) and ±1 and 2σ sensitivity (green and yellow band) on σ e v χ of fermionic DM absorption on electrons with the PandaX-4T commissioning data for the vector (upper) and axial-vector (lower) operators. Upper limits from the leading visible decays from x-ray satellites (cyan) whereas the leading process χ → γγγν for the vector operator is not shown which is above the parameter space selected, and cosmological constraints from leading invisible decays (dashed magenta) and DM overproduction (gray), and the best fit together with 1σ and 2σ contours associated with with 0.65 tonne·year XENON1T data (black) [30] are overlaid.
To summarize, we present the first sensitive experiment search on the absorption signals of fermionic dark matter on electron targets with an outgoing active neutrino with PandaX-4T. No significant dark matter signals are identified from the data. Our data present the strongest constraints on the fermionic DM in the mass range 35 to 55 (25 to 45) keV/c 2 for the vector opera-tor (axial-vector operator), in comparison to constraints from x-ray satellites and large scale observations, illustrating the strong physics potential of PandaX-4T. The best fit to the XENON1T excess with the same model lies marginally within our constraint, which calls for further investigation with more data [4,5,8,57].
We thank Tien-Tien Yu for the discussion on the ionization form factor, and Kenny C. Y. Ng and Brandon M. Roach for the discussion on the keV sterile neutrino. This project is supported in part by a grant from the Ministry of Science and Technology of China (No. 2016YFA0400301), grants from National Science Foundation of China (Nos. 12090061, 12090064, 12005131, 11905128, 11925502, 11775141), and by Office of Science and Technology, Shanghai Municipal Government (grant No. 18JC1410200). This project is also funded by China Postdoctoral Science Foundation (No. 2021M702148). We thank supports from Double First Class Plan of the Shanghai Jiao Tong University. We also thank the sponsorship from the Chinese Academy of Sciences Center for Excellence in Particle Physics (CCEPP), Hongwen Foundation in Hong Kong, and Tencent Foundation in China. Finally, we thank the CJPL administration and the Yalong River Hydropower Development Company Ltd. for indispensable logistical support and other help.
N ote added.--Recently, the XENONnT Collaboration released a new dataset [58], excluding XENON1T's excess. The new data are also consistent with our constraints.