Direct detection of freeze-in inelastic dark matter

We show that the current sensitivities of direct detection experiments have already reached the interesting parameter space of freeze-in dark matter models if the dark sector is in the inelastic dark matter framework and the excited dark matter state is cosmologically stable. Using results recently presented by the XENON1T experiment, we present constraints on these models. We also show that these models can explain the reported excess in the electron recoil signals if the mass gap between the ground state and the excited state is at keV scale.


Introduction
The particle nature of dark matter (DM) is one of the most prominent mysteries. Till now, all the evidence of the existence of DM is from gravitational effects. The relic energy density of DM in today's universe is measured to be about one-quarter of the total energy density. A successful DM model must be able to provide a mechanism to understand this number. The freeze-in scenario of DM production provides such a mechanism [1]. In this scenario, the DM particles live in the dark sector, very weakly connecting to the standard model (SM) sector through a portal. The portal is usually assumed to be a new vector kinetically mixed with the photon field or a new scalar very weakly coupled to the SM Higgs field. It is assumed that after inflation, only the SM sector is reheated, and through the portal, the energy in the SM sector leaked into the dark sector. In this scenario the observed relic density of DM can be nicely produced. However, the direct detection channel in freeze-in models is also proportional to the portal and, therefore, strongly suppressed. Inelastic DM models were first introduced to explain the excess observed in the DAMA/LIBRA experiment [2,3], since an enhanced annual modulation can be generated due to the extra cost of the kinetic energy in the up-scattering process. The down-scattering process in inelastic DM models is usually ignored since the population of the excited state is usually exponentially suppressed. In this work, we consider freeze-in inelastic DM models. We show that due to the possibility of the large down scattering rate, the sensitivity of the XENON1T experiment [4] has already achieved the interesting parameter space of inelastic freeze-in models. In Ref. [4], an excess around 1∼5 keV in electron recoil events is also reported, which cannot be accounted for by known backgrounds. Since the report of this excess, there have been active investigations trying to under it with new physics models . In this work we show that this excess can be explained in the framework of inelastic freeze-in models.
Vector portal freeze-in models In this paper, we consider two typical models, one is with a complex scalar dark matter and the other with a Dirac spinor dark matter. In both the models the dark portal is assumed to be a vector field V , we also call it dark photon in the following discussion. In both the models we assume the U(1) symmetry is broken softly by an explicit mass splitting. In the scalar case the masses of the real and imaginary parts are split, and in the Dirac spinor case the spinor is split into two Majorana spinors. The DM part of Lagrangian can be written as In L sc , χ is a complex scalar, χ 1 and χ 2 real scalars, whereas in L sp , χ is a Dirac spinor and χ 1 and χ 2 twocomponent Weyl spinors. Then the Lagrangian for the dark photon part is In the freeze-in process χ 1 and χ 2 must be produced in pair. As long as ∆m ≡ m 2 − m 1 is smaller than m V and 2m e , the dominant channels for the excited state χ 2 to decay to χ 1 are the three-photon channel and the neutrino channel. The former is suppressed by a factor of κ 2 (∆m/m e ) 8 and the latter is suppressed by 4 . Therefore in this case the lifetime of χ 2 can be much longer than the age of the universe. For the purpose of direct detection we assume ∆m Freezing-in the DM To calculate the relic density of DM we need to solve the Boltzmann equation where n D is the DM number density, H is the Hubble parameter, and Γ f i is the production rate of DM per volume. When the temperature T SM > 1 MeV, the universe is filled with relativistic plasma, therefore a Γ f i can be estimated as κ 2 e 2 D α em T 4 SM . The time interval at certain T SM can be estimated as H −1 ∼ m pl /T 2 SM , where m pl ≈ 1.22 × 10 19 GeV, is the Planck mass. Then the arXiv:2006.15672v1 [hep-ph] 28 Jun 2020 produced number density per entropy can be written as y D ≡ n D /s ∼ κ 2 e 2 D α em × (m pl /T SM ). Therefore in this scenario, DM is mainly produced at low temperate. The relation between y D and T SM stops when T SM hits either m V , m χ or m e . For m e < m V < 2m χ , which will be motivated later, the freeze-in process stops at T SM ∼ m χ . Consequently the dependence of today's relic density Ω D ∝ y D m D on m D is canceled. As an order of magnitude estimation, we can roughly get where m p is the proton mass. As a result, to get the observed relic abundance, the product of the dark coupling e D and kinetic mixing κ is fixed to κe D ∼ 10 −13 ∼ 10 −12 .
The numerical results of κe D required to produce the observed relic abundance for different choices of m V and m D are shown in Fig. 1. The direct detection rate is also proportional to the factor κ 2 e 2 D α em , making it difficult to search for the freeze-in model. It is also proportional to the number density of the DM, and therefore in favor of low mass DM as long as the energy deposit can surpass the thresholds of the experiments. Therefore in this work, we focus on the region where both m χ and m V are around MeV scale.
In this regime there are two main contributions to Γ f i , one is through e + e − annihilation, the other is through plasmon decay. The details of these two processes can be found in Ref. [41] (see also [42]). In our case since we require that m V < 2m χ the freeze-in cannot go through on-shell V the contribution from plasmon decay is always subdominant (The similar phenomenon is also found in the freeze-in process of on-shell dark photon [43]).
where f (x, y, s) = 1 2πy with a = e x/2 and b = y between between Γ sc f i and Γ sp f i is due to the difference of the spin structure and the factor of (1−4m 2 D ) 3/2 is due to the p-wave nature of the decay of virtual V into scalars.
Here and the following we use m D ≈ m 1 ≈ m 2 in the calculation when ∆m can be neglected.
Theoretical rates and recoil spectrum When χ 2 particles fly into the XENON detector, the ionization process through the down-scattering of χ 2 can happen. To calculate the ionization rate, we assume the scattering electron is approximated by a plane wave, and the xenon atom is isolated and described by the Roothaan-Hartree-Fock ground state wave functions [44]. It follows that the velocity averaged differential ionization cross section times velocity for electrons in the (n.l) shell can be written as (7) where E r and k are the kinetic energy and momentum of the outgoing electron, q is the momentum-transfer, and v is the velocity of the incoming χ 2 . According to the principles of quantum mechanics, scattering states and bound states from the same Hamiltonian must be orthogonal to each other. Therefore, in the domain that | q · x| 1 the plane wave approximation overestimates the cross section. To avoid this spurious contribution we subtract the bound state component from the outgoing wave function, that The form factors calculated in this way agree reasonably well with the ones used in [45]. In the case of downscattering, where E B is the absolute value of the binding energy. If E B ∆, the electrons indeed can be treated as free particles, and we can have an estimation of σv without going through the complicated form factors that where µ = m e m D /(m e + m D ). The differential cross section contributed from each energy level d σ nl v /d ln E r as a function of E r is shown in Fig. 2, where one can see that for each individual contribution the spectrum is bump-like and centered at the region m D ∆m/(m e + m D ). The width of each spectrum can be estimated as E r m e /m χ . This can be understand from the form factor in Eq. (7).
The differential ionization rate in the detector can be written as  and a signal spectrum with mD = 1.2 MeV, ∆m = 5 keV and a floating f2. The blue curve shows a signal with mD = 1.7 MeV, ∆m = 24 keV that corresponds to 95% CL limit shown in Fig.(4). .
We perform statistical analysis for both parts assuming the same test statistic as in Ref. [4] but without systematic uncertainties, which are neglectable comparing to the statistical ones. For the high energy sideband, we use the XENON1T data to constrain some parameter space. To ease numerical analysis, we assume that the data is consistent with a floating background component and the likelihood function is maximized by signal plus background events being equal to the data events. This assumption enables us to evaluate the test statistic without doing a fit and preserves its positivity provided there is no large signal-like excess. Using the asymptotic distribution of the test statistic in the large sample limit, we use the experiment spectrum in a range E R = ∆mm D m D +me ± 2∆E R to set observed 95% CL limit on Δm=8 keV, f2=0.5 Δm=16 keV, f2=0.5 Δm=24 keV, f2=0.5 Δm=8 keV, f2=0.2 Δm=16 keV, f2=0.2 Δm=24 keV, f2=0.2 Δm=8 keV, f2=0.05 Δm=16 keV, f2=0.05 Δm=24 keV, f2=0.05 the signal strength. The ∆E R is the detector resolution extracted from Ref. [46]. As shown in Fig.(4), the excluded region on the lower-left is dependent on the size of ∆m and f 2 . The constraints are stronger in regions of larger ∆m and f 2 . With different choices of ∆m and f 2 the parameter regions excluded by the XENON1T data are shown in Fig. 4. For the benchmark point in Fig.(3), the best fit signal with f 2 = 0.063 is found to be 3.7σ favored over the background-only hypothesis, which decreases to 2.7σ if the shape of tritium contribution is included as an unconstraint component. To find the parameter region consistent with this excess, we consider spectra whose E R ≈ m D ∆m/(m D + m e ) ≈ 2.7 keV and adjust f 2 such that the yield equals to that of the best fit benchmark. Results of the parameter scan are shown in Fig.(5) as contours of m V . We found that the model in this work can provide viable solutions to the XENON1T excess for dark photon masses between one and two MeV. The contour plot with different choice of m V from 2m e to about 2 MeV to fit the excess is also shown in Fig. 5, one can see that to get enough number events f 2 must be larger than about 0.02.
De-excitation in the early universe Although χ 2 is cosmologically stable, the χ 2 χ 2 → χ 1 χ 1 in the early universe can de-excite χ 2 into χ 1 . During freeze-in the temperature or "average kinetic" T D is comparable to T SM and is roughly equal to m χ . As ∆m m χ , the occupation of χ 1 and χ 2 are almost equal to each other. During expansion of the universe, T D evolves nonrelativistically. As a result we can parameterize T D as solving the Boltzmann equation. The Boltzmann equation for n 1 at T D ∼ ∆m m D ∼ m V can be written as where the collision coefficient C can be written as for both the scalar and spinor models introduced in Eq. (1). Here K 1 is the modified Bessel function with index 1. The decoupling of the de-excitation happens during T D ∼ ∆ T SM , so the universe is still at the radiation dominated era, and therefore the Hubble expansion rate is determined by the energy density of the SM sector. Defining f 2 = n 2 /n D , and x = ∆m/T D we have where and H 0 is today's Hubble parameter. The dependence of today's f 2 on A is shown in Fig. 6, where one can see that f today 2 ≈ 0.5 for A 1, and f 2 ≈ (2A) −1 for A 1.
MeV scale dark photon can be copiously produced during the supernova explosion, and therefore reduce the amount of the released neutrinos, which will conflict with the observations. For MeV scale dark photon, the constraint on the kinetic mixing can be written as [45] κ < 2.5 × 10 −9 (1 MeV/m V ) 2 .
For fixed m V and m D and by requiring all the DM particles observed today are produced from this mechanism, it can be translated to a lower bound on e D . Take the scalar model as an example. Eq. (16) means roughly e D > 10 −3 . Then from Eqs. (14) and (15) one can get that in this case which is about one order of magnitude smaller than needed to fit the excess (see Fig. 5). Similarly, for the spinor DM case, in the minimal model (1), f 2 is also about one order of magnitude needed. However, one can easily solve this problem by introducing more interactions in the dark sector. In the scalar DM scenario, a fourpoint interaction, which we omit in the Lagrangian (1), L λ = −λ(|χ| 2 ) 2 /2 cannot be forbidden by any symmetries. For the sake of vacuum stability, without inducing other interactions in the potential, λ must be positive. Then the combined square of the absolute value of the matrix element for the down scattering in the NR limit becomes One can see that these two contributions always have different signs. To reach the f 2 value required to fit the excess, only about 60% cancellation is needed. In the spinor DM case, one can introduce a scalar-carried force in the dark sector with the Lagrangian L h = −yhχχ, which equals to − 1 2 yh(χ 1 χ 1 − χ 2 χ 2 ) + h.c.. The down scattering process conducted by h is s-channel; its sign depends on the mass of h. Therefore, it is always possible for this new contribution to cancel part of the contribution from the exchange of V . In both the scalar and spinor DM scenarios, the newly introduced interactions in the dark sector contribute neither to the freeze-in process nor the direct detection.
Summary Freeze-in models with small couplings to the SM field is known to be challenging to search and constrain. However, in models that the DM is composed of a two-state system and with the excited state still populated in the universe, the direct detection signal can be enhanced in both the recoil energy and the rate. We show that the XENON1T experiment's sensitivity has already reached the exciting parameter space of these models. We also present a parameter region where the model predictions can explain the excess in the XENON1T ionization result.