Self-interacting Inelastic Dark Matter in the light of XENON1T excess

We propose a self-interacting inelastic dark matter (DM) scenario as a possible origin of the recently reported excess of electron recoil events by the XENON1T experiment. Two quasi-degenerate Majorana fermion DM interact within themselves via a light hidden sector massive gauge boson and with the standard model particles via gauge kinetic mixing. We also consider an additional long-lived singlet scalar which helps in realising correct dark matter relic abundance via a hybrid setup comprising of both freeze-in and freeze-out mechanisms. While being consistent with the required DM phenomenology along with sufficient self-interactions to address the small scale issues of cold dark matter, the model with GeV scale DM can explain the XENON1T excess via inelastic down scattering of heavier DM component into the lighter one. All these requirements leave a very tiny parameter space keeping the model very predictive for near future experiments.


I. INTRODUCTION
There exist convincing number of evidences suggesting the presence of a non-luminous, non-baryonic form of matter in the present universe, popularly known as dark matter (DM).
This form of matter constitute a significant portion of galaxies, clusters as well as the entire universe. Data from cosmology experiments like Planck which measures the cosmic microwave background (CMB) anisotropies very precisely, predict the amount of DM in the present universe to be around 26.8% of the present universe's energy density. In terms of density parameter Ω DM and h = Hubble Parameter/(100 km s −1 Mpc −1 ), the present DM abundance is conventionally reported as [1]: Ω DM h 2 = 0.120 ± 0.001 at 68% CL. Similar evidences exist in galactic and cluster scales as well, collected over a long period of time since 1930s [2][3][4]. It should be noted that the Planck estimate of present DM abundance relies upon the standard model of cosmology or ΛCDM cosmology which has been very successful in overall description of our universe at large scale (≥ O(Mpc)). Here CDM refers to cold dark matter while Λ denotes the cosmological constant or dark energy. CDM, a pressure-less or collision-less fluid acts like a seed for structure formation providing the required gravitational potential well for ordinary matter to collapse and form structures.
Since none of the standard model (SM) particles can be a viable CDM candidate, several beyond standard model (BSM) proposals have been put forwarded out of which the weakly interacting massive particle (WIMP) paradigm is the most widely studied one. In this framework, a WIMP candidate typically having interactions and mass in the electroweak regime, naturally satisfies the correct DM relic abundance, a remarkable coincidence often referred to as the WIMP Miracle [5].
While ΛCDM is in excellent agreement with large scale structure of the universe, yet there exist some discrepancies between its prediction and observations, particularly at small scales. In particular, too-big-to-fail, missing satellite and core-cusp problem are three such well known cases where ΛCDM appears to be in conflict with observations. For recent reviews of these issues and possible solutions, please see [6,7]. One interesting solution to this puzzle was proposed by Spergel and Steinhardt [8] where they considered an alternative to collision-less CDM in terms of self-interacting dark matter (SIDM) 1 . While SIDM solves the problems at small scales, it reproduces the CDM halos at large radii, thus consistent with observations. This is simply due to the fact that self-interacting scattering rate is proportional to DM density. The required self-interaction rate is often quantified as a ratio of cross section to DM mass as σ/m ∼ 1 cm 2 /g ≈ 2 × 10 −24 cm 2 /GeV [10][11][12][13][14][15]. Such selfinteracting cross sections can be naturally realised in models with very light mediator. For such a scenario, self-interactions can be shown to be stronger for smaller DM velocities such that it can have large impact on small scale structures while being consistent with usual CDM predictions at larger scales [10][11][12][13][16][17][18][19]. From particle physics point of view, such self-interactions can be naturally realised in Abelian gauge extensions of the SM. While DM sector can not be completely hidden and there should be some coupling of the mediator with SM particles as well, which can ensure that DM and SM sectors were in thermal equilibrium in the early universe. The same coupling can also be probed at DM direct detection experiments as well [20,21]. Several model building efforts have been made to realise such scenarios. For example, see [22][23][24][25][26][27] and references therein.
DM with light mediators have also received attention very recently after XENON1T collaboration published their latest results in June 2020 where they have reported the observation of an excess of electron recoil events over the background in the recoil energy E r in a range 1-7 keV, peaked around 2.4 keV [28]. While the excess can be explained by solar axions at 3.5σ significance or neutrinos with magnetic moment at 3.2σ significance both these interpretations face stringent stellar cooling bounds. While there is also room for possible tritium backgrounds in the detector, which XENON1T collaboration neither confirm or rule out at this stage, there have been several interesting new physics proposals in the literature.
Thus we noticed that in a class of models, the DM interpretation of XENON1T excess as well as SIDM phenomenology rely on light mediators. This motivates us to propose a common platform to show that the self interaction of DM arising via light mediators in such models can also give rise the observed XENON1T excess. In other words, the proposed scenario provides a unique way of probing the parameter space of SIDM at direct DM search experiments like XENON1T. To be more specific, we consider a dark sector consisting of sub-GeV inelastic DM with keV scale mass splitting and a corresponding massive vector boson Z [35,40]. Unlike earlier works where DM and Z masses are in the same regime so that DM relic is governed by resonant 2 → 2 annihilations, here we consider light mediators (order of magnitude lighter than DM mass) motivated from SIDM point of view. While the self interaction of DM is realised via Z -exchange, the latter can mix with U (1) Y gauge boson to provide a unique portal for detecting the DM at direct search experiments. The scalar field which leads to spontaneous breaking of dark sector gauge symmetry also induces a tiny Majorana mass to a singlet Dirac fermion field leading to an inelastic DM scenario [55,56].
In this setup we first find the DM parameter space consistent with velocity dependent selfinteraction rates explaining the data at the scale of clusters, galaxies and dwarf galaxies. We then confront the SIDM parameter space with the observed XENON1T electron excess while being consistent with other experimental bounds. We show that these two requirements make pure thermal relic DM insufficient to produce the observed relic and therefore we consider a hybrid setup where both freeze-out and freeze-in mechanisms can play non-trivial roles in generating DM relic. As we discuss in the upcoming sections, a long lived scalar singlet has to be invoked whose late decay into DM helps in generating correct DM relic in such a hybrid setup. This paper is organised as follows. In section II we briefly discuss our model followed by the analysis for dark matter self-interaction in section III. In section IV, we discuss production of self-interacting DM from a hybrid of freeze-in and freeze-out formalism. In section V we discuss the possible origin of XENON1T excess in our model. We finally summarise our results and conclude in section VI.

II. THE MODEL
We consider a simple Abelian extension of the SM. Under this U (1) X gauge symmetry, the SM fields do not have any charge while there exists a SM singlet Dirac fermion Ψ with U (1) X charge 1. A SM singlet scalar Φ having U (1) X charge −2 introduced which not only break the new gauge symmetry spontaneously but also splits the Dirac fermion into two pseudo-Dirac components as we discuss below. The Dirac fermion Ψ is identified as the DM field. The relevant part of the DM Lagrangian is where D µ = ∂ µ + ig Z µ and B αβ , Y αβ are the field strength tensors of U (1) X , U (1) Y respectively and is the kinetic mixing between them. The Lagrangian involving singlet scalar can be written as where H is the SM Higgs doublet. The scalar fields which acquire non-zero vacuum expectation value (VEV) can be represented as The singlet scalar VEV gives rise to U (1) X gauge boson mass M Z = 2g u while Higgs doublet gives rise to the usual SM particle masses.
The scalar singlet Φ also breaks U (1) X spontaneously down to a remnant where sin θ ≈ m − /M . The mass splitting between the two mass eigenstates is given by In order to address the XENON1T anomaly, we take ∆m ∼ 2 keV. While we stick to such minimal DM models in this work, such Abelian gauge extensions can be motivated from other phenomena like origin of light neutrino masses as well, as discussed in several works including [40,[57][58][59][60][61][62][63][64][65][66][67].

III. DARK MATTER SELF-INTERACTION
As mentioned before, we have an inelastic DM scenario where two DM components with tiny mass splitting of keV scale can populate the universe. We are also considering light mediator of DM interactions which is motivated from SIDM solution to structure formation problems. See [68][69][70] for earlier studies on self-interacting inelastic DM which were primarily motivated from the requirements of avoiding strong direct detection constraints or to explain some anomalous observations at indirect detection experiments like monochromatic photon lines. While in our model DM candidates can interact among themselves either via scalar or vector mediators, we consider only the vector mediator or Z to be light and hence consider this only to constrain the parameter space from required self-interactions. The relevant Lagrangian for DM interactions with Z can be rewritten as Ignoring the second term which is suppressed by m − /M , we can write down the corresponding potential for two Majorana fermion DM with a light mediator of dark photon type as [68][69][70][71] The two body Schrodinger equation for relative motion is where M ψ is the mass of the dark matter, ignoring the tiny mass splitting ∆m , v is the individual velocity of either of the dark matter particles in the centre of mass frame (half the relative velocity), ∆m is the mass splitting between two DM candidates and Ψ( r) is M ψ α and writing rΨ( r) = ψ(r), the s-wave Schrodinger equation is given by As shown in Feynman diagrams of figure 1, DM candidate of one type can scatter off each other while remaining in the same state, only at one-loop level, due to the off-diagonal nature of DM-mediator couplings. Using these, we constrain the DM parameter space from the required self-interactions at different scales while considering the mass splitting between the two DM candidates to be 2 keV, as favoured from XENON1T excess. The relevant cross sections are given in appendix A 1. For a more general analysis, one may refer to [68].
Using these self-interaction cross sections and using the required σ/m from astrophysical observations at different scales, we constrain the parameter space of the model in terms of DM (ψ 1,2 ) and mediator Z masses. As our study is motivated from explaining the XENON1T excess, we keep the required mass splitting between two DM candidates to be 2 keV. In figure 2, we show the allowed parameter space in DM mass versus Z mass plane which gives rise to the required DM self-interaction cross-section (σ/m) in the range 0.1 − 1 cm 2 /g for clusters (v ∼ 1000 km/s). The corresponding region of parameter space for galaxies (v ∼ 200 km/s)and dwarf galaxies (v ∼ 10 km/s) are shown in figure 3 and figure 4 respectively. It should be noted that for dwarf galaxies, due to smaller DM velocities we do not get sufficient self-interaction cross section (σ/m) from up scattering processes in the entire parameter space considered and hence the corresponding plot is not shown in figure 4. This is due to the fact that, lighter DM, due to low velocities, do not have sufficient kinetic energies to scatter efficiently into heavier DM resulting in a large self-interaction cross section. We will finally compare these regions of parameter space of GeV scale DM mass in the context of XENON1T excess and other phenomenological constraints.
The self-interaction cross section per unit mass of DM as a function of average collision velocity is shown in figure 5 as measured from astrophysical data. The data includes measurements from dwarfs (red), LSBs (blue) and clusters (green) [17,72]. here can explain the astrophysical observation of self-interaction of DM appreciably well.
See [73] for discussions on astrophysical probes of such inelastic DM with a light mediator. the mediator m ψ 1,2 > m Z , DM can have a large annihilation cross section to Z affecting its relic abundance. For example, the thermal averaged cross section for the t-channel process where α x = g 2 /(4π) and for typical gauge coupling and DM mass of our interest we have GeV. This leads to a cross section which is at least two order of magnitudes larger compared to the typical annihilation cross section of thermal DM. This reduces the relic abundance by same order of magnitudes, as seen from figure 7 showing the comoving number density of DM, assuming it to be a purely thermal relic. Before calculating DM relic, we first compare rates of different annihilation processes. Note that for the purpose of numerical analysis, the model has been implemented in LanHEP [74] and CalcHEP [75] and the cross-sections required has been fed into Mathematica [76] from CalcHEP. η † η → H † H, and decays back to DM after the dark sector freezes-out, thus filling the deficit in relic abundance. The relevant Lagrangian for η decay is given by: be larger than the observed DM abundance. As a result the late decay of η → ψ i ψ j can give rise to ample amount of DM. In Eq. 12, we use appropriate Boltzmann equations to get the correct relic density of DM. While we incorporate this additional scalar singlet η only to satisfy DM relic through its late decay, it can serve other purposes as well. One such possibility is to realise cosmic inflation. There have been proposals where a single scalar field is shown to play the role of inflation as well as thermal DM whose relic is generated via usual freeze-out. For example, see [78] and references therein. The same idea can be implemented here as well while noting that the scalar field is not perfectly stable but decays at late epochs into the DM. We however, do not discuss such additional roles the scalar singlet might play in our minimal scenario discussed here.
From figure 8, it is evident that the process DM e → DM e which is responsible for keeping both the dark and visible sector in kinetic equilibrium decouples around x ∼ 0.03, after which the temperature of the dark sector (denoted by T ) evolves independently of the thermal bath (temperature T) until x ∼ 100 when all the dark sector particles becomes nonrelativistic (and hence ceases to contribute to the relativistic degrees of freedom). Between these two epochs, the ratio of the two temperatures can be obtained by conserving the total entropy separately in the two sectors. Considering the kinetic decoupling temperature to be T D , we can relate the temperature of the two sectors as Here g SM * s (T ) is the relativistic entropy degrees of freedom in the standard model which goes into the calculation of relativistic entropy density s(T ) = 2π 2 45 g * s (T )T 3 . Since the above relation (10) is for T < T D , we naturally have g SM * s (T ) < g SM * s (T D ) leading to T < T . This is also understood from the fact that SM bath temperature receives additional entropy contributions from the species which keep getting decoupled gradually. Within the decoupled dark sector itself, the DM particles can transfer their entropies into lighter Z bosons once T falls below DM mass. This corresponds to an enhancement of dark sector temperature for T < m DM by (13/6) 1/3 , a factor close to unity. We have ignored this additional enhancement in the calculations.
Due to different temperatures of dark sector and SM bath after some epoch, we accordingly divide the range of integration for solving the Boltzmann equations as follows: • From the epoch of reaching kinetic equilibrium between DM-SM sectors till x < 0.03 (see figure 8), both the dark and the visible sectors share the same temperature T = T .
• One with 0.03 < x < 100 where the dark sector is decoupled from the thermal bath and its temperature evolves according to (10). Accordingly, one can define a new dimensionless parameter and relate to the usual parameter We can now write down the Boltzmann equations for two DM candidates ψ 1,2 and the scalar singlet η whose late decays into DM is crucial to generate correct DM relic. Unlike DM whose interactions with the SM bath are suppressed due to tiny kinetic mixing, the scalar singlet can be in thermal equilibrium with the SM due to large quartic couplings followed by freeze-out 2 . Thus, we define comoving number densities of these particles as Y ψ 1,2 = n ψ 1,2 /s (T (T )), Y η = n η /s(T ). The relevant coupled Boltzmann equations can then be written as where, ignoring the tiny mass splitting ∆m.
We solve these coupled Boltzmann equations taking into account of different temperatures of DM and SM sectors after kinetic decoupling, as given in (10). The corresponding evolutions of different comoving number densities are shown in figure 9. In figure 9, the dot-dashed dark blue line shows the equilibrium number density of the singlet scalar η with mass m η ∼ 1 TeV, which was initially in thermal equilibrium with the SM bath. As its interaction rates falls below the expansion rate, it freezes out leaving a thermal relic, shown by the green dot-dashed line, assuming it to be stable. The blue dot-dashed line shows the freeze-in production of DM only from the process e + e − → DM DM without considering subsequent annihilation of DM into Z pairs. When we take into account both its production from e + e − → DM DM and subsequent annihilations into Z bosons via DM DM → Z Z its abundance is depicted by the pink line. The sharp contrast is due to the strong DM DM → Z Z annihilation rate which reduces the abundance of DM produced from freeze-in. As the number density of DM increases due to freeze-in production, the annihilation rate into Z pairs also increases leading to the first depletion in the pink line around x = 0.1. Shortly after that, DM production from freeze-in again balances DM annihilation rate leading to a plateau region all the way till x = 1. However, since freeze-in production from thermal bath becomes negligible beyond x = 1, we see further depletion in DM density due to its annihilation into Z pairs leaving an under-abundant relic beyond x = 10. Note that, at this point we have not considered scalar decay contribution to DM.
Since freeze-in production of DM from the thermal bath followed by DM annihilation into Z pairs lead to under-abundant relic density, we now consider the additional contribution from scalar singlet decay. The red dot-dashed line shows the evolution of comoving number density of DM after taking scalar decay contribution into account. The corresponding evolution of the scalar number density is shown by the maroon coloured dot-dashed line.
Clearly, once the number density of the scalar falls due to its decay, the DM number density gets uplifted. Once the decay is complete, DM relic also saturates beyond x ≈ 30. It should be noted that, the scalar decay occurs after DM annihilation to Z pairs freezes out around x = 10 to avoid further depletion. Also, while considering freeze-in production of DM from the thermal bath, we considered the contribution of electron-positrons only, for simplicity.
If we consider all the particles in the thermal bath, we will get more freeze-in production of DM and the final required abundance of DM can be realised by appropriate tuning of scalar decay width without affecting rest of the analysis related to self-interaction and XENON1T excess.
Note that the lines showing the evolution of DM number density in figure 9 considers both the DM components ψ 1,2 . Since their mass splitting is very small ∆m ∼ O(keV) they behave very similarly as far as calculation of relic abundance goes. However, once the net relic is generated, there can be interconversion between two DM components dominantly through Z -mediated t-channel process ψ 2 ψ 2 → ψ 1 ψ 1 . We take this into account and show that the effect of such interconversion with such small mass splitting (∆m = 2×10 −6 GeV) is negligible. This can be seen from figure 10, where the fractional contributions Y DM 1 /Y DM Total and Y DM 2 /Y DM Total for mass splitting ∆m = 2 × 10 −6 GeV are shown. We have also taken into account the Sommerfeld effect induced by the multiple Z boson exchange in the interconversion process [80]. Clearly, such interconversions lead to negligible effects on individual DM relic abundance and hence we consider them to be equally dominant in rest of our analysis.

V. THE XENON1T EXCESS
The direct detection prospects of such self-interacting DM can be addressed through the recently reported excess in the electron recoil events at XENON1T experiment. We assume ψ 2 is heavier than ψ 1 with a small mass splitting ∆m = M 2 − M 1 between the two components. Because of this inelastic nature of these DM candidates and since the mass splitting ∆m is kept fixed at keV scale, it can successfully explain the recently reported XENON1T anomaly [28]. For a fixed incoming velocity v of heavier DM ψ 2 , the differential scattering cross section for the down scattering process ψ 2 e → ψ 1 e (with electrons inside the Xenon atom) can be written as where m e is the electron mass, σ e is the corresponding free electron cross section at fixed momentum transfer q = 1/a 0 with a 0 = 1 αme being the Bohr radius, α = e 2 4π = 1 137 being the fine structure constant, E r is the recoil energy of electron and K(E r , q) is the atomic excitation factor. For our calculations, the atomic excitation factor is adopted from [81].
We assume the DM form factor to be unity.
However, to include velocity dispersion in Eq. (13), we use the following distribution function (obtained after angular integration of a Maxwellian velocity distribution boosted in earth's rest frame) velocity dispersion of DM, can be rewritten as [38,51,81,82] d σv Where v esc is the DM escape velocity in the Milky Way which is of the order v esc ∼ 533 +54 −41 km/s [83]. In inelastic DM scenarios, the minimum DM velocity (v min ) required by the DM to upscatter to the NLSP and register a recoil inside the detector is decided by the kinematics of scattering. However, it is worth mentioning that in the case of an inelastic down scattering of DM with electron, which we consider here, there is no kinematic limit on the minimum velocity of DM as the incoming particle with almost vanishing velocity can still down scatter to the lighter component with the mass splitting between the DM components being transferred to the electron recoil energy, without violating anything kinematically.
The free electron scattering cross-section for the process ψ 2 e → ψ 1 e is given by where α Z = g 2 4π , α = g 2 4π and is the kinetic mixing parameter between Z and Z gauge bosons. For chosen values of DM and mediator masses in our work, this kinetic mixing is required to be ∼ 10 −8 . It should be noted that, for GeV scale DM, σ e is independent of DM mass as the reduced mass of DM-electron is almost equal to electron mass. The limits of integration for the inelastic scattering in Eq. (15) are determined depending on the relative values of recoil energy (E r ) and the mass splitting between the two DM components.
And for E r ≤ ∆m The dependency of atomic excitation factor on the momentum transferred q is shown in figure 12. Here the dominant contribution comes from the bound states with principal quantum number n = 3 as their binding energy is around a few keVs. In the right panel of figure 12, we have shown the plot for the integration of momentum transferred times the atomic excitation factor i .e.K int (E r , q) = q+ q− qdqK(E r , q) as a function of the recoil energy E r for M 1 = 0.3GeV and ∆m = 2keV. The figure shows a peak around E r ∆m since the q − approaches to zero and the momentum transfer maximising this factor is available. It is worth mentioning that such kind of enhancement is a characteristic feature of inelastic scattering.
The differential event rate for the inelastic DM scattering with electrons in Xenon atom, i.e ψ 2 e → ψ 1 e, can be given as: where n T = 4 × 10 27 Ton −1 is the number density of Xenon atoms and n DM is the number density of the dark matter particle.
The detected recoil energy spectrum can be obtained by convolving Eq. (19) with the energy resolution of the XENON1T detector. Incorporating the detector efficiency γ(E), the energy resolution of the detector is given by a Gaussian distribution with an energy dependent width, where γ(E) is reported in figure 2 of [28] and the width σ det is given by with a = 0.3171 and b = 0.0037. Thus the final detected recoil energy spectrum is given by On the other hand, in the bottom panel of figure 13, we have shown the fit considering different velocity dispersion for the DM particle as we have no observational constraints on f (v) apart from numerical simulations. Clearly as we increase the velocity dispersion the peak in the spectrum giving an appreciable fit gets flatten out and no longer explain the XENON1T signal within E r = 2 − 3keV for larger σ v .

VI. SUMMARY AND CONCLUSION
We summarise our key findings in figure 14. We show all the relevant constraints as well as favoured parameter space in the g −M Z plane. In figure 14, all the coloured regions (except the blue one which favoured from XENON1T excess) represent disfavoured regions from different bounds. The green patch represents the region where the DM self-scattering crosssection is not large enough to solve the astrophysical problems discussed in section III. To be more quantitative, the green shaded regions correspond to DM self-scattering cross-section σ/m < 0.1 cm 2 /g. The triangular region on upper left corner of figure 14 is disfavoured from lower bound on lifetime of heavier DM. Since the mass splitting between ψ 1 and ψ 2 is kept at keV scale ∆m = O(keV ), there can be decay modes like ψ 2 → ψ 1 νν mediated by Z − Z mixing. If both the DM components are to be there in the present universe, this lifetime has to be more than the age of the universe, that is τ ψ 2 > τ Univ. . The decay width of this process is Γ(ψ 2 → ψ 1 νν) = g 2 g 2 2 (∆m) 5 This tiny kinetic mixing also prevents DM from reaching chemical equilibrium with the SM requiring its non-thermal or freeze-in production from the SM bath. However, due to large coupling of DM with U (1) X gauge boson Z , they can annihilate strongly into much lighter Z bosons depleting the number density generated from freeze-in. To fill the gap, we introduce another long-lived scalar singlet which freezes out from the thermal bath and decays very late into DM generating the required relic. As seen from the summary plot in figure 14, matter ψ 1 and ψ 2 , which we assume to constitute about 1% of the total relic (see figure   7). Since the relic density is smaller by two orders of magnitude than the observed one, the corresponding DM-electron cross-section σ e (ψ 2 e → ψ 1 e) has to be increased by two orders in order to explain the observed XENON1T excess. This can be achieved by increasing by one order of magnitude, since σ e ∝ 2 . However, increasing by one order of magnitude will not satisfy the lifetime bound on ψ 2 as τ The scattering cross sections can be derived as [68] σ ψ 1 ψ 1 →ψ 1 ψ 1 = π where we have defined, ∆ = 2 v − 2 δ , µ and V 0 are defining parameters for the exponential potential V 0 e −µr , given by, Here r M is chosen from the relation e − φ r M /r M = max( 2 δ /2, 2 φ ). The terms Γ v , Γ ∆ are given by with Γ denoting the gamma function.

DM Velocity distribution function
The distribution function used in Eq. 14 can be obtained as follows. Let − → u and − → v are the velocities of dark matter in the rest frames of galaxy and earth respectively. If − → v E is the velocity of earth with respect to the galactic rest frame then we have − Assuming that the velocity distribution of dark matter with respect to the galactic rest frame is Maxwellian, we can write where N is the normalisation constant and σ v is the velocity dispersion. Assuming spherical symmetry and considering z-axis in the direction of − → v E which subtends an angle θ with − → v , we can write : . Now carrying out the integration for the angular co-ordinates φ and θ, we obtain where we have neglected e