Improved limits on solar axions and bosonic dark matter from the CDEX-1B experiment using profile likelihood ratio method

We present improved constraints on couplings of solar axions and more generic bosonic dark matter particles using 737.1 kg-days of data from the CDEX-1B experiment. The CDEX-1B experiment, located at the China Jinping Underground Laboratory, primarily aims at the direct detection of WIMPs using a p-type point-contact germanium detector. We develop the profile likelihood ratio method for analysis of data in the presence of backgrounds. The background modeling is compatible with the data and no excess signals are observed. An energy threshold of 160 eV was achieved. This significantly improve the sensitivity for the bosonic dark matter below 0.8 keV. Limits are also placed on the coupling $g_{Ae}<2.26 \times 10^{-11}$ from Compton, bremsstrahlung, atomic-recombination and de-excitation channels and $g^{eff}_{AN} \times g_{Ae}<4.14 \times 10^{-17}$ from a $^{57}$Fe M1 transition at 90\% confidence level. All the constrains improve over our previous results.


I. INTRODUCTION
For the CP problem of strong interactions, the Peccei-Quinn (PQ) mechanism [1] is still the most compelling solution in which a new kind of U(1) symmetry would be spontaneously broken at large energy scale f A . After this original solution to the CP conservation in QCD, a new Nambu-Glodstone boson called axion is proposed later by Weinberg [2] and Wilczek [3] through the PQ symmetry. Axions are pseudoscalar particles with properties closely related to those of neutral pions and their mass m A is fixed by the scale f A of the PQ symmetry breaking, m A ≈ 6 eV ( 10 6 GeV/f A ). The range of scale f A can not be restricted by theory but the order of the electroweak scale has been excluded by experiments. At a higher symmetry-breaking energy scale, 'invisible' ax-ion models such as hadronic model KSVZ (Kim-Shifman-Vainstein-Zakharov) [4,5] and non-hadronic model DFSZ (Dine-Fischler-Srednicki-Zhitnitskii) [6,7] are still allowed. Another interest in this paper is more general bosonic dark matter (BDM) like pseudoscalar axion-like particles (ALPs) and vector bosons, which also have couplings to electrons.
A new 1−kg scale detector with lower energy threshold was fabricated for the CDEX-1B experiment. A data size with exposure of 737.1 kg−days is adopted for this analysis. This same data set has been used in earlier WIMP search analysis with integrated rate [34] and annual modulation [30] signatures, and with the presence of Midgal effects [29]. We describe in this article the statistical model with profile likelihood ratio method applied to this CDEX-1B data, and report the results on the searches of solar axions and BDM.

A. CDEX-1B setup and overview
The CDEX-1B experiment adopts one 939g singleelement PPCGe crystal with dead layer of 0.88 ± 0.12 mm [35]. Outside of the PPCGe detector is the passive shielding system and the detailed information is described in Ref. [34]. A well-shaped cylindrical NaI(Tl) crystal surrounding the PPCGe detector is used as the anti-Compton detector. The coincidence events both in germanium and NaI(Tl) crystals denoted as AC + are discarded to depress the γ background.
The schematic diagram of electronics and data acquisition (DAQ) system is shown in Ref. [34]. Four identical energy-related signals were out of the p + point-contact electrode after a pulsed-reset feedback preamplifier. Two of them were distributed into the shaping amplifiers at 6 µs (SA 6µs ) and 12 µs (SA 12µs ) shaping time for low energy region (0-12 keV). The output of SA 6µs provided the system trigger of the DAQ. The other two outputs were fed into timing amplifiers (TA) which provide the accurate time information. One with high gain (TA 1 ) is limited to medium energy region (0-20 keV), and the other one with low gain (TA 2 ) for high energy can reach 1.3 MeV. The energy resolution of TA 1 output is similar to the SA 6µs,12µs . As a result, the spectrum below 12 keV is from SA 6µs and above 12 keV is from TA 1 in our analysis. The energy resolution (σ) at 10.37 keV is about 75 eV in this analysis.

Solar Axions
The sun is a potential source of axions and in this article we concentrate on two different mechanisms.
The first is the 14.4 keV monochromatic axions from the M1 transition of the 57 Fe in the sun, i.e. 57 Fe * → 57 Fe+A, due to the stability and the large abundance of 57 Fe in the sun.
The Lagrangian coupling axions to nucleons is [12]: where ψ N is the nucleon isospin doublet, φ A is the axion field, and τ 3 is Pauli matrix. g 0 AN and g 3 AN are the modeldependent isoscalar and isovector axion-nucleon coupling constants [36,37]. Introducing g ef f AN ≡ (−1.19g 0 AN +g 3 AN ) as the effective nuclear coupling adapted to the case of 57 Fe, the corresponding axion flux is given by [12,38]: where κ A and κ γ are the monenta of the outgoing axion and photon respectively. Given the axion-nucleon couplings g 0 AN and g 3 AN for specific models such as DFSZ and KSVZ, the axion flux can be evaluated.
The second from the Compton-like scattering(C), axion-bremsstrahlung(B), atomic-recombination(R) and atomic-deexcitation(D) processes, and the corresponding effective Lagrangian is given by [12]: where g Ae is the dimensionless axion-electron coupling constant. Its flux depends on the g 2 Ae : where the units of fluxes are cm −2 s −1 keV −1 and axion energy E A is in unit of keV. For the RD process, the tabulated spectrum in Ref. [39] is used. As discussed in Ref. [39], the flux is valid for light axion; hence, we consider only the axion mass below 1 keV.
The axion-electron coupling is depended on models. In the DFSZ model, the coupling is proportional to cos 2 β, where tanβ is the ratio of the two Higgs vacuum expectation values. In the KSVZ model, it depends on E/N , the ratio of electromagnetic to color anomalies. E/N = 0 and cos 2 β = 1 are used in this analysis [12].

Bosonic Dark Matter
The main cosmological interest in bosonic particles such as pseudoscalar ALPs and vector BDM arises from their possible role as the dominant component of dark matter, the nature of which is still unknown. The absorption via ionization or excitaiton of an electron in target atom makes BDM experimentally interesting and PPCGe detectors have advantages to study BDM due to their excellent energy resolution, sub−keV threshold and low radioactivity background. Decay modes of these two kind of dark matter are different. For pseudoscalar ALPs, the two−photon decay is the dominant channel, while this channel is strictly forbidden for vector bosons, whereas decay into three photons is allowed at loop level [13].
Assuming that these bosonic particles constitute all of the galactic dark matter, we get the total average flux of dark matter axions on Earth: where ρ DM ∼ 0.3 GeV/cm 3 is the dark matter halo density [40], m A is the axion mass, v A is the mean axion velocity distribution with respect to the Earth and β is the ratio of the axion velocity to the speed of light for cold dark matter. This flux is independent of any axion coupling.

C. Particle interactions in CEDX-1B
The axion detection channel studied in this paper is the axio-electric effect illustrated in Eq. (1). The axioelectric cross-section as described in Ref. [41][42][43] is given by: where σ pe (m A ) is the photo-electric cross-section for germanium in the unit of barns/atom, m A is the mass of axion, α is the fine structure constant, m e is the electron mass and β is the ratio of the axion velocity to the speed of light. The expected axion rates of CBRD process and 57 Fe are displayed in the Fig. 1.
In the situation of pseudoscalar ALPs in cold dark matter model (β ≈ 10 −3 ), the coupling to electrons of which is same with the case of QCD axions. For the vector BDM, the absorption cross section σ abs can be written as: where m v is the mass of the vector BDM, α and α are the fine structure constant and its vector boson equivalent, respectively. Using the parameter mentioned above, the interaction rate in the direct detection experiment can be written as: for pseudoscalar ALPs and for vector BDM, where A is mass number of germanium. The expected rates of these two kinds of particles are shown in Fig. 2.  The expected event rate of pseudoscalar ALPs (red solid line) and vector BDM (blue solid line) at different masses. The red dashed line is the maximum event rate of ALPs Gaussian distributions versus their masses, while the blue dashed line is corresponding to vector BDM. The couplings used here are gAe = 2 × 10 −11 and α /α = 5 × 10 −25 . The widths of these peaks are determined by the energy resolution.

A. Data Selection
As discussed in earlier analysis [34], the background spectrum is derived by the following steps: 1. Stability check, removing the time periods of calibration or other testing experiments.
2. Anti-Compton veto, disgarding the the events in coincidence with the anti-Compton detector and retains the anti-coincidence events.
3. Basic cuts, removing the electronic noise through getting rid of the abnormal pulses and spurious signals.
4. Bulk and surface event selection, rejecting the surface events by pulse shape analysis using their characteristic slower rise time.
The detailed information of these selection efficiencies and trigger efficiency are discussed in Ref. [34]. An improved Ratio Method, which base on the bulk/surface rise time distribution PDFs (probability density functions), is developed to reject the surface event and derive the signal-retaining and background-leakage efficiencies [44]. This method has been proved correctly above 160 eV. So in this analysis, 160 eV is selected as the physics analysis threshold, at which the combined efficiency including those from the trigger, basic cuts and anti-Compton veto is 17%.

B. Background and Understanding
With an exposure of 737.1 kg-days, the bulk spectrum from 160 eV up to 20 keV after data selection and efficiency correction is displayed in Fig. 3(a). The background consists of several K shell X-rays and their corresponding L shell X-rays from the cosmogenic isotopes and a continuous background with a smooth, slightly increasing profile as the energy decreases [34]. Considering the low muon flux mentioned above, the contribution from muon can be neglected. The continuous background below 20 keV was expected to probably originate from the the 238 U, 232 Th and 40 K in the materials in the vicinity of the PPCGe detector, radon gas penetrating through shielding and cosmogenic 3 H in the crystal. A detailed modeling of the continuous background is beyond this work and will be studied in our future work.
However, axion analysis is not sensitive to the accurate background assumption because the signatures of axion are significantly different from the continuous background. As can be seen from Fig. 2, the signal signatures of 57 Fe and BDM are monochromatic and of Gaussian distribution with widths determined by the energy resolution. As to the continuous CBRD solar axion, a saw-tooth-like profile arises between 0.9 keV and 1.6 keV considering the axion mass below 1 keV. So in the following fitting procedure, the background model can be described by a continuous background plus the peaks from K/L-shell X-rays. Benefiting from the low threshold and better energy resolution of CDEX-1B, the L-shell X-ray peaks at low energy region can be clearly distinguished. So in the background model, the amplitude of the K-shell X-ray peaks and the corresponding L-shell X-ray peaks are limited by each other using the K/L-shell X-ray ratios mentioned in Ref. [45,46]. In the ultra low energy region around the threshold, M-shell X-rays are also taken into consideration in the background model.
Besides the bulk events, surface events also enter the later analysis as a kind of background and the surface event distribution in energy is displayed in Fig. 3(b).

C. Axion Analysis Method
A profile likelihood analysis, as described in Ref. [47], is used to get the constraints and the test statistics is: where L is the likelihood function. Quantity µ is a parameter corresponding to the strength of signals and θ denotes all of the nuisance parameters. The quantitŷ θ denotes the value of θ that maximizes L for the specified µ, while the denominator is the maximized likelihood funciton, i.e. µ and θ are their ML estimators. To obtain the 90% C.L. bounds on the signal strengths µ, we have to calculate the probability distribution functions (p.d.f) of the test statistics q µ . Let f (q µ |µ) be the p.d.f of the test statistic q µ under the signal strength hypothesis µ, and f (q µ |0) be the p.d.f of q µ under the background-only hypothesis. Since downward fluctuations of background might lead to a much stringent exclusion result, we used the CL s method [48] to get rid of this effect. The 90% up limits µ up are defined as: where F is the cumulative distribution function of the test statistic.
The specific likelihood function we used in this analysis can be written as follows and the parameter of interest becomes the axion number N A instead of µ mentioned above: where g b,j (E i ) and g s,j (E i ) are the bulk and surface event p.d.f of rise time at the certain energy bin i, while σ b,j (E i ) and σ s,j (E i ) are their corresponding errors. Parameter f b represents the normalized p.d.f. of the background including the continuous component and K/L/M shell X-ray peaks. Parameter f s represents the normalized p.d.f. of the efficiency-corrected surface spectrum S r in Fig. 3(b), derived from fitting S r by a smooth curve, while f A is the p.d.f of signal. All these three distribution functions are modified by the data selected efficiency parameters e. The parameter n ij is the measured event number in the energy spectrum bin i and the rise time spectrum bin j, while the detailed information of the 2-dimension grid is displayed in the Fig. 4(a). N bulk and N surf are the corrected number of events in bulk region and surface region after bulk/surface discrimination. Take the system uncertainties caused by efficiency calculation and bulk/surface discrimination into consideration, we used parameters e, t b and t s to construct the additional terms L 2 and L 3 .
where e 1 , e 2 are parameters in an error function used to describe the efficiency line of Basic cut, while e 3 , e 4 are used in another error function to describe the trigger efficiency. The center values µ ei are the best-fit results used in Ref. [34] and V ij is the covariance matrix between center values.
Parameters t b and t s describe the entire systematic deviation from bulk and surface ratio after selection, while t b,s = ±1 corresponds to ±1σ .

A. 14.4 keV Solar Axion
The signal of solar axions produced in the 57 Fe magnetic transition on the spectrum is a monochromatic Gaussian peak around 14.4 keV with width determined by resolution, which is about 84 eV (σ) under this situation. The fitting range is limited to 14.06 keV to 14.76 keV, about ±4σ, and a polynomial function is used to described the background in this range. The 90% C.L result is shown in Fig. 5 and the rate of this kind of axion is found to be less than 0.029 counts·kg −1 ·d −1 . For a low-mass axion at 0 keV, this result translates to a 90% C.L. constraint on the coupling:

B. CBRD
For CBRD solar axions, the fitting range is from 0.8 keV to 2.0 keV, because there is a saw-tooth-like profile arising in this energy range which is different from the continuous background. Using the analysis procedure mentioned above, we get the constraints on g Ae :  [18], EDELWEISS-II [12], EDELWEISS-III [20], Majorana Demonstrator [14] and PandaX-II [16]. The yellow band represents the 1σ expected sensitivity. The CDEX-1B 90% C.L. on CBRD solar axions (solid red line), together with astrophysical bounds [8,11] and other direct search experiments [12, 15-17, 20, 23]. The yellow band represent the 1σ expected sensitivity.
better results of BDM comparing with CDEX-1. The 90% C.L. limits on g Ae of ALPs and α /α of vector BDM are displayed in Fig. 9 and Fig. 10 respectively. Due to the much lower energy threshold, we can extend the first point of exclusion line down to the 185 eV.  9. The CDEX-1B 90% C.L. upper limit on coupling of pseudoscalar ALPs as a function of mALP, together with the constraints set by CDEX-1 [18] and other experiments [9, 10, 12, 14-17, 20, 21]. The yellow band represents the 1σ expected sensitivity.

V. SUMMARY
Based on the data size of 737.1 kg-days from CDEX-1B experiment, we set constraints on the couplings of solar axions and BDM using profile likelihood ratio method. The results of solar axions are better than CDEX-1, while the constraints on ALPs and vector BDM improve over the previous results below the mass of 0.8 keV due to the low energy threshold and excellent energy resolution. The 90% C.L. upper limit on the coupling of vector BDM from CDEX-1B together with the result of EDELWEISS-III [20], Majorana Demonstrator [14], XENON100 [19], TEXONO [21] and astrophysical bounds from [13]. The yellow band represents the 1σ expected sensitivity.
The impacts of systematic uncertainties have been taken into consideration for all of analysis presented here. The major source of systematic uncertainties includes the background assumption, energy resolution, signal selection efficiency and bulk/surface discrimination. The uncertainties on selection efficiency and bulk/surface discrimination are taken into account in the profile likelihood function via the nuisance parameter e and t. The uncertainties of the background assumption are evaluated by using different background model. The variation of background models causes the change of constraints less than 5% for CBRD axion, less than 20% for BDM, and less than 8% for 57 Fe solar axion. Varying the energy resolution by ±10%, the change of result is less than 17% for 57 Fe solar axion, less than 9% for BDM and negligible for CBRD axion.