Explorations of pseudo-Dirac dark matter having keV splittings and interacting via transition electric and magnetic dipole moments

We study a minimal model of pseudo-Dirac dark matter, interacting through transition electric and magnetic dipole moments. Motivated by the fact that xenon experiments can detect electrons down to $\sim$\,keV recoil energies, we consider $\mathcal{O}$(keV) splittings between the mass eigenstates. We study the production of this dark matter candidate via the freeze-in mechanism. We discuss the direct detection signatures of the model arising from the down-scattering of the heavier state, that are produced in Solar upscattering, finding observable signatures at the current and near-future xenon based direct detection experiments. We also study complementary constraints on the model from fixed target experiments, lepton colliders, supernovae cooling and cosmology. We show that the latest XENONnT results rule out parts of the parameter space for this well motivated and minimal dark matter candidate. Next generation xenon experiments can either discover or further constrain how strongly inelastic dark matter can interact via the dipole moment operators.

The scattering of typical WIMP-like DM, with velocity distributions following the Standard Halo Model, against electrons lead to recoil energies of µ e,DM v 2 DM ∼ O(eV) 1 , where µ e,DM is the electron-DM reduced mass and v DM is the incoming DM velocity. The O(keV) recoil energies that XENON1T probes currently can broadly arise from various DM models like boosted DM [45,53,97,122,144] that acquire higher velocities than the typical DM halo velocities, absorption of O(keV) particles [146], or through down-scattering of inelastic DM with two nearly degenerate DM states with mass splittings of O(keV) [9,23,51,111,120,125]. The coincidence of similarity in the recoil energies probed by XENONnT and the temperature at the Sun's core (∼ 1.1 keV) motivates the study of inelastic DM which is upscattered in the Sun. Specifically, if the heavier particle is not cosmologically stable, and the lighter particle constitutes the DM, they can be excited to the heavier particle by upscattering against Solar electrons. It can subsequently down-scatter via electron scatterings at XENONnT, depositing an energy approximately equal to the mass splitting. 1 See ref. [33] for a discussion on the recoil energies in DM scattering against electrons inside a nucleus. For comparison, a DM particle with a form factor of FDM (q) = 1 is shown to be unable to explain the excess without being in conflict with XENON1T S2-only analysis. While DM with form factors FDM (q) ∝ q, q 2 mediated by heavy particles are shown to be give recoil energies of O(keV).
With this in mind, we study the most minimal model of pseudo-Dirac DM with two Majorana fermions. The lowest dimension interactions with SM allowed, in the absence of any additional dark sector particles, are the transition electric and magnetic dipole moments. Dipolar DM [32,130] with scattering processes mediated by the photon, gives enhanced scattering rates in the small velocity limit applicable in terrestrial direct detection (DD) experiments [107]. We therefore study: • The freeze-in (FI) production [106,131] of pseudo-Dirac DM with mass splittings of O(keV).
• The DD of this DM via up-scattering in the Sun followed by down-scattering in electron recoil events at the XENONnT experiment (as well as projections for future runs of XENONnT and DARWIN). We study the DD rates for generic parts of the parameter space without the relic density constraint from FI production.
• Complementary bounds on generic parts of the parameter space. These include constraints from fixed target experiments, lepton colliders, SN1987A and N eff constraints.
Inelastic DM with interaction via transition dipolar moments have previously been studied with focus on different parameter spaces [130,137,148]. Refs. [46,95,114] have also studied pseudo-Dirac type DM in a dark photon model. Dipolar DM motivated by lepton sector minimal flavor violation was recently studied in [67].
In section II, we discuss our model and the production of the DM via FI mechanism. In section III, we discuss the signatures at DD experiments from electron and nuclear scatterings. In section IV, we discuss complementary bounds on the model from various existing experiments and observations. Finally, we present the results in section V and conclude in section VI.

II. MODEL FRAMEWORK AND FREEZE-IN PRODUCTION
We consider a dark sector consisting of two Majorana fermions χ 1 and χ 2 that form a pseudo-Dirac state with mass splitting δ ≡ m χ 2 − m χ 1 ∼ O(keV). The lowest dimension interaction operators allowed in the absence of any additional new particles are through the DM dipole moments where the interaction is mediated by photons. Since the only dipolar type interactions that are allowed for production from plasmon decay.
Majorana fermions are of transition type, we get the following Lagrangian for DM interaction [130] where µ χ and d χ are the transition magnetic dipole moment (MDM) and electric dipole moment (EDM), respectively.
We consider the production of DM by the FI mechanism which is operative when the DM coupling to SM is small enough for the DM to have never entered thermal equilibrium with the SM bath. With an assumption of negligible DM number density at the earliest epoch, it is produced via annihilation and decay of SM particles; some relevant Feynman diagrams for dipolar DM are shown in fig. 1. As the temperature of SM bath falls, the DM production ceases as follows. For m DM > m SM , as the bath temperature falls below the DM mass, the DM production becomes kinematically suppressed.
While for m SM > m DM , as the bath temperature falls below the SM particle mass, it freezes out and its annihilation rate is suppressed. Thus the DM number density per comoving volume becomes a constant, and the DM is said to have frozen into the current relic abundance.
We solve the Boltzmann equation to calculate the DM relic density. For convenience, we define the number per co-moving volume as Yield, Y ≡ n/s, with n representing the number density and s being the entropy density. The Boltzmann equation then gives [80,106,131] dY dT where T is the temperature of the SM bath, M P l is the Planck mass, g s * and g * are the temperature dependent relativistic degrees of freedom contributing to entropy and energy densities, respectively. R(T ) is the rate for DM production and is a sum of rates of production from 2 → 2 annihilation, R 2→2 , and from the decay of plasmon, R γ * , We discuss these two modes of production in the following subsections.
The total DM relic density is where we have assumed only the lighter particle χ 1 with mass m χ 1 to be stable on cosmological timescales. Here, ρ crit is the critical density and T 0 is the temperature of current epoch, with the values for the physical constants taken from ref. [24].

A. Annihilation production
The rate of production of particles χ 1 , χ 2 from the annihilation of SM fermions, as shown in fig. 1 (a), is [28,80,106] where N f c are the color degrees of freedom of SM fermion f with mass m f . The total rate is a sum over the SM fermions f . K 1 is the modified Bessel function of the second kind, s min = Max (m χ 1 + m χ 2 ) 2 , 4m 2 f , and the 3-momentum of χ 1 in the center-of-mass (CM) frame is The squared-amplitude |M| 2 is summed over initial and final spins [130] dΩ ⋆ where we give the expressions upto leading order in δ ≡ m χ 2 − m χ 1 . In the limit of vanishing SM fermion and DM masses, we find that the production rate scales with temperature as R(T ) ∝ T 6 .
From eq. (2), we see that this leads to UV (UltraViolet) FI with maximum yield produced at the   largest temperature of the SM bath, equal to the reheating temperature for models with instantaneous reheating [52,60,85,100].
Of phenomenological interest is the parameter space with large coupling which can lead to large scattering rate at DD experiments. We note that for m DM ≥ T RH , there is kinematic suppression in DM production, requiring large couplings to reproduce the observed relic density. For the light DM we consider in this work, m DM ≤ 1 GeV, we choose low reheating temperatures (T RH = 5 MeV and 10 MeV) from allowed 2 values. We verify our calculations using the expressions given above, with those obtained from micrOMEGAs 5.0 [24], for fermionic channels of production.
We discuss the results in section V.

B. Plasmon decay
The SM bath consists of charged particles that couple to the photon. An electromagnetic wave propagating through this bath behaves qualitatively differently from that in vacuum. Coherent vibrations of the electromagnetic field and the density of charged particles results in a spin-1 particle with 1 longitudinal and 2 transverse polarizations, propagating at a speed less than the speed of light in vacuum. Its dispersion relation depends on the properties of the SM plasma and is therefore known as the 'plasmon'. We follow the discussion in refs. [38,140] for the plasmon properties and decay rates.
The significance of plasmon decay leading to production of DM in FI scenarios was noted in ref. [83] (see also [50,62,108]). The massive plasmon can decay to DM in cases where the DM couples to the photon. This is true in the dipolar interaction case we study here. For DM produced via FI mechanism, the final abundance is accumulated over time and summed over the various production channels of SM particles annihilating and decaying to DM. Therefore, the production from plasmon decay can play a significant role in FI produced DM (as opposed to freeze-out production).
2 A robust lower bound on reheating temperatures of TRH > 4 MeV was derived by combining cosmic microwave background, large scale structure and light element abundances data in ref. [110].
To calculate the DM abundance resulting from plasmon decay, we must begin with the modified dispersion relations which depend on the temperature T and the net density of the charged particles in the SM plasma. Since ref. [38] was addressing the energy loss from the plasma in a supernova with typical energies of O(10) MeV, only electrons and positrons were considered as charged particles of relevance that modify the dispersion relations. But for the case of a UV FI with large enough T RH , all charged particles with masses much smaller than the reheating temperature would add to the plasmon effect (see appendix A for further details). The rate of DM production from plasmon decay is where we assume that only the lighter particle χ 1 survives. The transverse and longitudinal polarizations are g T = 2 and g L = 1, respectively. The distribution for the plasmon is given by The decay width of a plasmon Γ pol with four momentum k = (ω, ⃗ k) in the medium frame, with a definite polarization 'pol' is [140]  with |M| 2 γ * pol →χ 1 χ 2 being the squared amplitude for plasmon decay to DM, expressions for which are given in eq. (A11). For dipolar DM we find that the rate of DM production from plasmon decay is also maximised at the largest temperatures of the SM bath. We call this a UV plasmon production.
We show in fig. 2, the scaled DM yield as a function of temperature, from plasmon production for dipolar interaction (UV plasmon) and millicharged interaction [83,84] (Infrared, IR, plasmon). Here, we have considered a DM of mass 1 keV for both the cases. The UV plasmon production can be seen to be maximum at the reheating temperature, taken to be 1 GeV in this figure. As T RH is changed, the IR produced yield is unchanged. While the UV plasmon line would shift towards right (higher temperatures) as T RH is increased .
Although we have only shown the yield for EDM interacting DM in fig. 2 for the UV plasmon, the same is also true for MDM interaction. Note that we show the UV plasmon production for T RH = 1 GeV to clearly show the difference in IR vs UV production; the maximum production of DM happens at the largest temperatures even for smaller values of T RH . This is in contrast with IR plasmon production where the maximum production takes place at the lowest temperatures so that the plasmon production mechanism begins to dominate for small DM masses, m DM ≲ 400 keV [83]. This can be understood by noting that the 2 → 2 annihilation process becomes ineffective when electrons freeze out, while the plasmon production continues to be effective.
While for UV plasmon production (with small T RH ), the plasmon decay production does not take over the annihilation production as we decrease the DM mass, since both the processes maximize at the largest temperatures of the SM bath. For the small reheating temperatures we consider, mindful of the sub-GeV DM masses that are of interest to the DD discussed in the following, the plasmon production is always subdominant and can be safely ignored. For much larger reheating temperatures T RH ≳ 100 GeV though, the plasmon production can dominate and must be taken into account.

III. DIRECT DETECTION
We discuss the DD of inelastic DM with interactions via transition EDM and MDM operators. As mentioned above, an O(keV) splitting is of interest for production of the excited state in the Sun [23], as well as for detection at electron scattering experiments. Additionally, mediation by the photon leads to an enhanced scattering cross section at low velocities [108].
For couplings of interest that reproduce the observed relic density, the heavier state χ 2 is not stable on cosmological scales and χ 1 makes up the full DM density. We assume that this is also true for other parts of the parameter space. For the given model, a scattering process can either proceed through the inelastic scattering process of χ 1 e → χ 2 e or through a loop-suppressed elastic scattering process, χ 1 e → χ 1 e. The up-scattering process requires larger DM velocities than those allowed in the Galactic halo, while the elastic scattering rate is loop-suppressed. The only avenue 3 for DD is then from the up-scattering of χ 1 against electrons in the Sun to χ 2 [7,8,23]. The χ 2 that have outgoing velocities large enough to overcome the gravitational potential well of the Sun, can then reach the Earth (as depicted in fig. 3) and down-scatter at DD experiments via χ 2 e → χ 1 e. In the following, we follow the discussion in refs. [7,23,87] to find the event rates at DD experiments.

A. Up-scattering from the Sun
The maximum velocity of DM falling into the Sun is (v esc where v esc ⊙ is the Solar escape velocity. We take the value of the Solar escape velocity to be the number averaged 4 escape velocity, over the Solar core [87] (since most of the scatterings are expected to happen inside the core), giving v esc ⊙ ≃ 1307 km/s [87]. DM particle's most probable Galactic velocity 5 is given by v halo DM ≃ 220 km/s. This is much smaller than the most probable electron velocity in the Solar core v e ⊙ (core) = 2T ⊙ m e ≃ 0.066, for T ⊙ = 1.1 keV, 3 Cosmic ray upscattering rates will be suppressed since the constituent particles have relativistic speeds and the dipolar scattering rates are inversely proportional to the DM-target relative velocity (see eq.(29)). 4 ⟨vesc⟩ = rcore 0 dr vesc(r)n(r)/ rcore 0 dr n(r) where n(r) is the electron number density [87] and rcore ≃ 0.2R⊙. 5 Using a halo distribution averaged velocity leads to an overall factor of ≃ 0.89 in the upscattered flux.
where T ⊙ is the temperature at the Sun's core. Hence, the DM particles can be understood to be at rest, with Solar electrons scattering against them. The steady state DM number density in the Sun is given by Here, the gravitational focusing effect, that enhances the area at spatial infinity, is given by the factor , and the factor (v halo DM /v esc ) accounts for the decrease in the number density of DM owing to its larger velocity near the Sun (from conservation of flux). The velocity distribution of electrons in the Sun is taken to be Maxwell Boltzmann (MB) where m e and v e are the electron mass and velocity, respectively. The differential flux of particles χ 2 generated with recoil energy K χ 2 is given as where n e = 2 × 10 25 cm −3 and V ⊙ = 2.2 × 10 31 cm 3 are the Solar mean electron number density and volume, respectively [19,23]. The factor of 1/4π(1AU) 2 is the scaling in the flux on account of traveling from the Sun to the Earth, with the distance travelled being 1 AU (astronomical unit). The reduced mass is given by Here, we have plugged in the explicit expressions for the temperature averaged differential cross section for up-scattering, given by where v min (K χ 2 ) is the minimum velocity an electron must have in order to up-scatter a χ 1 into a χ 2 , with recoil energy K χ 2 , given as and in the limit δ << m e and m χ 1 . Here,σ e is the DM-electron reference cross section [91] with the 3-momentum transfer fixed at q = (αm e ) (appropriate for atomic processes), The q-dependence of the matrix elements are encoded in the DM form-factor F DM (q) [91]. |M χ 1 e (q)| 2 is the squared-amplitude for DMelectron scattering, averaged (summed) over initial (final) spin states.
The MDM contribution consists of one part that is independent of relative velocity v and with a form factor F DM (q) = 1 similar to that for contact interaction, and another part proportional to v 2 and with a form factor same as that of EDM.
Note that the χ 2 particles that reach the Earth are the ones that overcome the gravitational potential Therefore, the flux on the Earth is attenuated from the produced flux as The Heaviside theta function ensures that the condition in eq. (20) is satisfied. We also shift the flux to account for the reduction in χ 2 velocity in overcoming the Sun's gravitational potential. In the following, we denote the flux on the Earth by dΦ/dK χ 2 and drop the explicit notation. There will be a further suppression in the flux of χ 2 from its decay in the time taken to travel from the Sun to the Earth. In fig. 4 we show the flux per unit energy from eq. (14) at production (blue-dashed), the flux that escapes the Sun's gravitational well from eq. (21) (blue-solid) and the attenuated flux at the Earth after taking into account the χ 2 decay in travelling from the Sun to the Earth (black-solid).
Note that the flux that overcomes the gravitational potential has a lower cut-off, from the Heaviside theta function in eq. (21).
We note that the plasmon does not play any role in the up-scattering or production of DM in the Sun. The plasmon frequency in the non-relativistic limit [38] gives a plasmon mass order 1 smaller than the Sun's average temperature for the Solar electron number density. Therefore, the plasmon effect is subdominant in scattering. In addition, there can be no plasmon-sourced production as we consider DM masses much larger than the Solar temperature.
B. Down-scattering on the Earth

Electron scattering:
The χ 2 flux received on the Earth is depleted if the decay lifetime of χ 2 is comparable to the time taken by it to travel to the Earth. We take this into account, and calculate the down-scattering event rate in electron scattering experiments. With the DM flux per unit energy dΦ/dK χ 2 given in eq. (21), we can write the electron recoil energy spectrum per detector mass per unit time as [23,33,89] where ∆E e is the energy transferred to the electron 6 , E nl is the ionization energy of the n, l orbital of given atom, ϵ(∆E e ) is the signal efficiency given in fig. 2 of [13]. The maximum energy transferred to the electron for a given incoming DM kinetic energy K χ 2 is given by ∆E max DM-electron scattering in xenon, we get the number of targets per tonne as n T = 2×4.2×10 27 /tonne 7 .
The exponential factor accounts for the depletion in χ 2 flux in travelling from the Sun to the Earth and t(K χ 2 ) = 1AU × m χ 2 /2K χ 2 is the time taken by χ 2 to travel to the Earth. The decay width of The form factor for ionization of an electron in n, l orbital with a total of ∆E e energy transferred to the electron, is given by |f nl→∆Ee−E nl (q) 2 . We use QEDark [89] to extract these ionization form factors.
The limits of integration q ± are obtained from energy conservation in the down-scattering process and given as 6 The energy transferred to the electron is a sum of the energy of the outgoing electron at asymptotically large distances from the nucleus, ER, and the ionization energy of the shell it originated from, E nl , i.e., ∆Ee ≡ ER + E nl . The energies are assumed to be emitted almost simultaneously and the collection of the energies of the electrons and photons emitted at the de-excitation and the ionization together is assumed to be equal to ∆Ee, that is we assume that the ionization energy is released completely [116]. 7 We take Z eff = 2 since the electrons in different orbitals are accounted for by summing over ionization form factors for all the accessible orbitals where θ is the angle between the momentum of χ 2 and the transferred momentum q and We note that the factor from integration over transferred momentum q in eq. (23) leads to an enhancement in the F DM (q) = 1 part of the MDM scattering rate but a small suppression for the remaining part of MDM as well as the EDM scattering rate.
We find the limits from latest results from XENONnT experiment [14] with a total exposure of 1.16 tonne × year and background rate of 15.8 ± 1.3 events/(tonne × year × keV). We also find the projected limits from future runs of XENONnT and DARWIN experiments, with projected exposures of 20 tonne × year and 200 tonne × year, respectively. We use the same background rate for future projections, as that of the latest XENONnT run, to get the most conservative limits. We also use the same efficiency as of XENON1T, which gives a conservative bound, since the efficiency can only be expected to improve. We discuss the results in section V.

Scattering from Migdal effect:
In addition to the direct ionization of an electron, a DM particle scattering off of a nucleus can also lead to subleading electronic energy deposition into detectors via nuclear scattering. In a typical nuclear recoil the electron cloud is assumed to follow the recoiling nucleus instantaneously. But if the effect of the sudden acceleration of the nucleus, with the electron cloud still in its original position, is taken into account, it is known to deposit electronic energy via ionization/excitation of the recoiling atom (the Migdal effect) or the emission of a Bremsstrahlung photon [12,27,116,133].
In the Migdal approximation, the whole electron cloud is assumed to recoil with the same velocity with respect to the nucleus, with no change in its shape. The rate for a nuclear scattering with recoil energy E NR accompanied with a Migdal electron recoil with energy E e from n, l orbital, leading to a deposition of total energy ∆E e ≃ E e + E nl , is [26,74,116] where, E det = E e + E nl + LE NR . The Lindhard quenching factor is denoted by L and is the fraction of the nuclear recoil energy observed in the electron channel. Its value is well approximated to L ≃ 0.15 [26,89]. The ionization form factor is given by Z(E det − LE NR − E nl ) 2 and its values for different orbitals are given in fig. 4 of ref. [116]. The nuclear differential scattering rate is: The nuclear differential scattering cross sections, approximated to the elastic case, are where m N and Z are the mass and atomic number of the nucleus, respectively, and E NR is the nuclear recoil energy.
The most sensitive low-energy analysis comes from the S2-only data set from the XENON1T experiment [11]. The S2-only differential rate is where the minimum kinetic energy of the incoming DM χ 2 to downscatter to χ 1 , with a nuclear recoil energy of E NR along with an electronic deposition of energy E det via Migdal effect, is given by where, v Mig,inel We find the total number of event at XENON1T by integrating eq. (30) over the range E det = 0.19−3.8 keV ee , with an exposure of 22 tonne-day [11]. A total of 61 events were observed at XENON1T over this exposure, while the expected number of background events was 23.4. This gives an upper limit of 49 events expected from DM at 90% confidence, and can be used as an upper limit to derive constraints on DM interactions with SM. For our model of dipolar DM though we find less than 1 events over the full parameter space of interest, so the scattering rate from Migdal effect is not large enough to derive any bounds on Solar upscattered dipolar DM.
We note that this is expected from the discussion in [92] for long range interactions, like EDM F DM (q) = 1/q 2 , with a bias towards low q values where the Migdal effect rates are always smaller than electron ionization rates by a factor of ∼ Z 2 m 2 e /m 2 N .

IV. COMPLEMENTARY BOUNDS ON SUB-GEV DARK MATTER
In this section, we discuss the constraints on sub-GeV DM with O(keV) splittings from existing laboratory experiments and astrophysical sources.

A. Fixed target experiments
At lepton-fixed target experiments, an electron beam of fixed energy is dumped against an active target, comprised of some heavy nucleus, that is either a part of the detector itself or a separate target.
Dark sector particles can be pair produced from the electrons scattering off of nucleons, giving rise to missing-energy final states [10,101]. These experiments employ stringent selection criteria making it possible to conduct an essentially background free search for such missing-energy signals [3,31,93,103,147]. In particular, we use results from the NA64 experiment at CERN SPS [20] to derive constraints on our model via the production of DM particles from an off-shell photon, where γ vir is the off-shell photon that subsequently produces the dark sector particles. The leading processes are shown in fig. 5. Note that we do not consider the processes where the virtual photon originates from the nucleus since these processes are suppressed by a factor of (Zm e /m N ) 2 for coherent photon emission. We follow the discussion in ref. [63] in the following. Assuming that the fermion pair is produced within the first radiation length gives the target length to be L target = X 0 = 0.56 cm for Pb. The search region in the ECAL is limited by the energy threshold for detection of electron on one side and the requirement for missing energy to be larger than half the beam energy. This gives the selection criteria for the energy of the outgoing electron E 4 as [21,123] 0.3 GeV < E 4 < 50 GeV.
The polar angular coverage for the outgoing electron is 0.0 < θ 4 < 0.23 radians.
The number of signal events with these geometric and angular cuts are [63] where, N EOT = 2.84 × 10 11 are the number of electrons incident upon the target [21]. The target material density and nuclear mass are given by ρ target and m N , respectively. The detector efficiency of NA64 is known to depend only marginally on the energy and is taken to be constant, ϵ ef f (E 4 ) ≃ 0.5 (averaging over the total signal efficiencies of the various runs [21]). The integration limits of E 4 and θ 4 are from eqs. (34) and (35), respectively. The double differential cross section for production of the processes shown in figs. 5a and 5b is given by dσ prod /dE 4 d cos θ 4 with the full expression given in eq. (C3).
Since the beam energy is much larger than the O(keV) splitting we consider, and the signal being observed is of missing energy in final state, the constraints on our model do not differ from those of the elastic DM cases discussed in ref. [63]. We follow the discussions in [63] and re-derive the constraints from their run as mentioned in ref. [21]. The details of the calculation are given in Appendix C. We derive constraints by demanding that N sig < N 90% = 2.3 where the latter corresponds to the 90% C.L. for the number of signals events given zero observed events. The resulting constraints are shown in figs. 7 and 11.
This type of search for missing-energy in the final state gives stronger bounds for a feebly interacting dark particle, as compared to experiments where the dark particle is detected via its scattering off electrons/nuclei in the main detector (for example in the mQ experiment at SLAC [139]), since processes of the latter kind are suppressed by further powers of the small-valued DM-SM effective coupling. We, therefore, do not study these latter processes.
Constraints are also derived from proton fixed target experiments, and as shown in ref. The FSR processes are suppressed by the small DM-photon couplings, so we only consider the ISR process, as shown in fig. 6, for deriving constraints.
We follow the discussion in ref. [62] for constraints on DM dipole moments from pair production, adapting them for the inelastic case. For DM mass splittings much smaller than the CM energy of the χ 1 , χ 2 system, δ ∼ O(keV) << √ s χ 1 χ 2 , these constraints for inelastic DM are equal to the ones for elastic DM.

BABAR
We use the data from the search for mono-photon events in decays of Υ(3S) at BABAR detector at the PEP-II asymmetric-energy e + e − collider at the Stanford Linear Accelerator Center (SLAC), with 28 fb −1 of data collected at a CM energy [16] √ s ≈ m Υ(3S) ≈ 10.355 GeV.
The single-photon events were chosen based one two trigger criteria, High-E region: 3.2 < E γ < 5.5 GeV, − 0.31 < cos θ γ < 0.6 radian, full data, Low-E region: 2.2 < E γ < 3.7 GeV, − 0.46 < cos θ γ < 0.46 radian, 19 fb −1 of data, where θ γ is the CM polar angle, and E γ is the photon energy in the Υ rest frame. For each region, the number of signal events is given by FIG. 6: Processes for production of DM along with initial state radiation, at e + e − colliders, leading to missing energy plus mono-photon signatures.
Here, ϵ eff is the total efficiency, L is the integrated luminosity, √ s is the CM energy of the e + e − system, √ s χ 1 χ 2 ≡ (1 − x γ )s is the CM energy of the χ 1 χ 2 system with The photon makes an angle of θ γ with respect to the beam in the CM frame. Following ref. [90] we apply a non-geometric cut ϵ eff of 30% and 55% in the high-E and low-E regions, respectively.
The ISR production cross section is approximated by dressing the cross section for DM pair production (without ISR), with an angle-dependent radiator function [62,134] as: where the cross-section for e + e − → χ 1 χ 2 at the energy scale reduced by photon emission is = α(s + 2m 2 e ) 6s 3 and, the radiator function taking into account all soft and collinear corrections upto O(m 2 e /s) is [134] H α (x γ , θ γ , s) = α π To constrain the DM couplings, we require that the expected number of events be smaller than the observed number of events in each bin i at 90% C.L., such that N obs . We show these bounds in figs. 7 and 11.

LEP
We consider high energy colliders like the Large Electron Positron (LEP) collider. The model considered here can give rise to events with one photon and missing energy by the same process as at BABAR, see section IV B 1. The only difference lies in the high CM energies 189 GeV ≤ √ s ≤ 209 GeV and high luminosities (a total of 619 pb −1 of data for the single-and multi-photon + missing energy final states) that the LEP operated at [2]. The cross sections for these processes having been found to be in agreement with the SM expectation from e + e − → ννγ(γ) give constraints on any BSM physics model that can also lead to the same final state. But the high CM energies lead to constraints that are independent of the DM mass, for light DM (m χ < GeV here). This was studied in ref. [98] and bounds obtained from the mono photon channels give where µ B ≡ e/2m e is the Bohr magneton. We note that at these large energies of production, the small splitting of O(keV) has negligible effect and the constraints for such inelastic DM are the same as those for the elastic case. The upper bound is shown in figs. 7 and 11.
Note that the LHC probes heavier masses and smaller couplings ∼ O(0.01) GeV −1 [22], and is not of relevance in this work.
since positrons are thermally supported in SN. Constraints on models are derived in two limiting cases leading to two-sided bounds on DM couplings for each mass as follows: 1. Weak coupling: This is applicable in the limit of small interaction strengths of DM with SM such that for any smaller strengths there would be too little production of DM to cause any significant change to the SN cooling rate. In this limit, any DM that is produced escapes the SN with almost 100% probability, such that it is possible to derive lower bounds on the DM effective coupling, by constraining only the production rates. This is given by the "Raffelt criterion" [140] which says that any "exotic" cooling will not change the neutrino signal significantly, as long as the emissivity obeys the conditionĖ < 10 19 erg/g/s.
This condition is easily converted into a condition on energy emitted per unit time per unit volume by noting that the density for t ≃ 1s is nearly constant (see fig. 5 of ref. [42]), giving a conversion between d 3 r and dM . We take ρ ≈ 3 × 10 14 g/cm 3 . The emissivity (energy emitted per unit volume per unit time) is defined as [75] where f i are the Fermi-Dirac distribution functions We ignore the final state Pauli blocking for small DM number densities and simplify the expression using where we use energy conservation E 1 + E 2 = E 3 + E 4 and F ≡ s(s − 4m 2 e )/2 [102]. The cross-section of production for the process in eq. (43) with squared-amplitude summed over initial and final spins 10 , is given in the limit of vanishing electron mass m 2 e < s ≈ T 2 avg ≃ (30 MeV) 2 as 10 σ ′ = g1g2σ where g1, g2 are the spin degrees of freedom of incoming SM fermions, and σ is the usual cross-section defined with the squared-amplitude averaged(summed) over initial(final) spins.
Also simplifying the remaining part of eq. (45) as [80,102] with E ± ≡ E 1 ± E 2 , we rewrite eq. (45) in the limit of vanishing electron mass as , computed at radius r 0 = 10 km, where the emmissivity can be seen to be maximum. We use the temperature and chemical potential radial profiles as given in ref. [129]. The limits of E − integration are E max,min 2. Large coupling: In the opposite limit of large DM couplings, the cooling process is dictated by the probability of escape or mean free path (MFP) of DM. The relatively larger density of SN results in a trapping of DM particles produced inside the SN, giving an upper bound on the DM effective coupling with SM. We adapt this bound from ref. [63], shown in figs. 7 and 11.

V. RESULTS
We summarise the results for sub-GeV pseudo-Dirac DM with mass states χ 1 and χ 2 having mass difference O(keV) for the two interactions: 1. Transition EDM interaction: • As discussed in section II A, the FI production is UV sensitive. We show the relic density contours for two reheating temperatures, 5 MeV and 10 MeV, in fig. 7. For higher values of reheating temperatures, the contours would shift to lower d χ values and the upward bend would occur at higher m χ 1 , going outside the range shown here. The T RH = 5 MeV contour can be seen to cut off at m χ 1 ≃ 100 MeV as the observed relic density cannot be reproduced via FI production for larger masses. We show the relic contours for two values of mass splittings, δ = 1 keV (dashed) and δ = 10 keV (dotted). These coincide at small masses for a given reheating temperature, since δ << m χ 1 . But they begin to diverge for larger masses, m χ 1 ≳ T RH , as the Boltzmann suppression leads to an exponential that is more sensitive to the difference of the two masses.
We can see that the bounds from SN1987A are applicable on parts of these contours.
• The total number of events at various xenon based DD experiments are shown in the color palette in fig. 7. The points shown correspond to a mass splitting of δ = 1 keV. We find that masses less than 12 MeV and 4 × 10 −6 GeV −1 ≲ d χ ≲ 10 −5 GeV −1 are ruled out by XENONnT. These results are to be understood to be correct upto O(10%) since we ignore astrophysical uncertainties (solar parameters and DM halo distribution) in probing orders of magnitude of DM masses.
• We note that there is a competition between the decay rate of χ 2 and its down-scattering cross section at DD experiments. This is because increasing d χ gives larger scattering cross sections as seen from eq. (19), but also increases the decay width, Γ χ 2 ∝ d 2 χ , so that the flux received on the Earth decreases. Therefore, we see from fig. 7 that starting from the smallest values of d χ , the rate initially increases with increasing d χ , maximizing at some value of d χ and then falls quickly with further increase in d χ .
• This competition also leads to a maximum DM mass that can be probed at DD experiments. The XENONnT experiment currently probes masses m χ 1 ≲12 MeV, while in the future XENONnT can probe m χ 1 ≲ 18 MeV and DARWIN can probe m χ 1 ≲ 23 MeV. We note that a large part of the points probed by XENONnT lie in the parameter space ruled out by N eff bounds [50] for this model. These arise because of thermalization of the dark sector for large enough couplings, assuming standard cosmology 11 .
• Larger parts of the parameter space that aren't ruled out by N eff bounds will be probed by future runs of XENONnT and DARWIN. Further, we show the event shapes for some benchmark points in figs. 8, 9 and 10, showing in red the signal+background rates, and the blue band corresponds to background rates with Poissonian uncertainties (± √ N ). We see that the signals may show up as an excess in the lowest recoil energy bins, and can also be distinguished from the background by the shape of the spectra. shaded region is ruled out by N eff constraints in standard cosmology [50]. Also shown are the contours that lead to observed relic density for T RH = 5 MeV (brown) and 10 MeV (black). The dashed and dotted lines for each correspond to mass splittings of δ = 1 keV and 10 keV, respectively.
The points with color palette show the total number of events for mass splitting δ = 1 keV for various xenon experiments, as mentioned below each figure.
We can see the decrease in event rates when δ is increased from 1 keV in fig. 8 to 1.5 kev in fig. 10. This is because as δ increases, the decay width of χ 2 increases, Γ χ 2 ∝ δ 3 . In addition, the up-scattering rates in the Sun get suppressed since the only electrons with enough energy to cause the up-scattering are those from the high velocity tail of MB distribution, given an average temperature of T ⊙ ≃ 1.1 keV. Together, these prohibit the detection of large mass splittings (δ ≳ 1.2 keV at XENONnT current run, δ ≳ 1.7 keV at XENONnT future projection and δ ≳ 1.8 keV at DARWIN) via electron scattering.
• The minimum d χ values that can be probed by DD experiments are such that the parameter space that leads to the production of observed relic density by FI production, is not within the reach of the DD experiment. While smaller reheating temperatures could lead to an overlap between the two, its value is bounded from below by T RH > 4 MeV [110].

Transition MDM interaction:
Here, we only discuss the features that are distinct from the EDM case, with all other features being the same. shaded region is ruled out by N eff constraints in standard cosmology [50] for this model. These arise because of thermalization of the dark sector for large enough couplings, assuming standard cosmology. Also shown are the contours that lead to observed relic density for T RH = 5 MeV (brown) and 10 MeV (black). The dashed and dotted lines for each correspond to mass splittings of δ = 1 keV and 10 keV, respectively.
• We note from eqs. (17) and (19) that the MDM differential cross section dσ/dE R ∝ {v 2 /E R , 1} (corresponding to the two form factors F DM = 1, 1/q, see eq. (19)) which is suppressed with respect to the corresponding EDM dσ/dE R ∝ 1/v 2 E R by factors of Here we use the fact that q 2 ≃ 2m e E R and the suppression comes from v and E R being small. Therefore, the scattering rates for MDM interaction are highly suppressed and do not lead to significant event rates at the xenon DD experiments, thus are not shown in fig. 11.

VI. CONCLUSIONS AND OUTLOOK
We have studied a model of inelastic DM interacting with SM via transition electric and magnetic dipole moments. We first address the production of DM by the FI mechanism taking into account both 2 → 2 annihilation production and production from decay of plasmon. We find that both these processes are UV dominant with the rate of production (or dY /dT ) maximised at the largest temperatures (near T RH ). We observe that the plasmon production is insignificant for the parameter space we are interested in here, although it can become the dominant source of production for much larger reheating temperatures.
The heavier states are not stable on cosmological scales and their flux is produced by Solar upscattering of the lighter, stable χ 1 particles that are assumed to make up the entirety of DM. We study the constraints on this model from DD experiments where we can observe the down-scattering of the heavier mass state in electron recoil events. We find that DM with masses less than 12 MeV and transition EDM 4 × 10 −6 GeV −1 ≲ d χ ≲ 10 −5 GeV −1 are ruled out by XENONnT [14], already probing parameter space not ruled out by any other constraints.
In addition, we find that future results from DD experiments, XENONnT and DARWIN, by virtue of larger exposures and lower backgrounds, can lead to the discovery of pseudo-Dirac DM with EDM interaction and mass splittings less than 2 keV. Notably, this parameter space is not probed by any current experiment. The reach of DD experiments will further improve by lowering of detection thresholds (as suggested by the S2 only analysis from the XENON1T experiment [11]). We also show complementary constraints on the transition dipolar DM model from current fixed target experiments, e + e − colliders and information from SN cooling. Projections from Belle-II show that additional parameter space will be probed in the future [63].
We thus study the most minimal model of inelastic DM with electric and magnetic dipolar couplings.
With a focus on xenon based direct detection experiments we find that future DD experiments have a great potential to discover this minimal model. Our work provides further motivation for an in-depth exploration of low-energy electron recoil events in xenon based DD experiments.

VII. ACKNOWLEDGEMENTS
We thank Itay Bloch, Jae Hyeok Chang, Xiaoyong Chu, Hyun Min Lee, Hongwan Liu, Tarak Nath Maity, and Mukul Sholapurkar for correspondence and helpful discussions. RL acknowledges financial support from the Infosys foundation (Bangalore), institute start-up funds, and Department of Science and Technology (Govt. of India) for the grant SRG/2022/001125.

Appendix A: Plasmon Production
The plasmon frequency is a good measure of the magnitude of medium effects and is given as [38,140]: where the sum is over contribution from SM fermions ψ. The velocity of SM particle ψ is given by v = p/E and the Fermi-Dirac distributions corresponding for SM fermions and anti-fermions are given by f ψ and fψ, respectively. To arrive at the second line, we assume that the chemical potentials are zero and that each antiparticleψ stops contributing at temperature Tψ, which we take to be Max[2mψ, Λ QCD ]. The first mode frequency of the plasma is given by: using which we can define a velocity v * = ω 1 /ω p . If we consider ψ ∈ {e}, then v * can be understood as the typical electron velocity. With these definitions, the general dispersion relations for the transverse and longitudinal polarizations are given by the following approximate expressions [38] that are correct to order α. Here, k max is the maximum wavenumber upto which longitudinally polarized plasmons can be populated The in-medium couplings of the photon to the SM particles are modified by vertex renormalization constants Z T,L given by [140] Since the dispersion relations of transverse and longitudinal polarizations of the thermal photons are distinct, we separate the two polarizations in the follwoing calculation. The decay width of a plasmon with four momentum k = (ω, ⃗ k) in the medium frame, and a definite polarization is [140]: where the squared amplitude, summed over incoming and outgoing spin states is The first term for each case in eq. (A11) is integrated in the rest frame of the plasmon as shown in eqs. (B5-B9), giving: with C 1 as given in eq. (B10). The integration for the last three terms are done using Lenard's formula [126] updated for the massive, inelastic case Lenard's formula for our case is [126]: Multiplying both sides by g µν and contracting, we get where k 2 = s. Here, the integration on RHS has been carried out in the rest frame of plasmon: where the ⃗ p χ 1 integral just sets ⃗ p χ 1 = ⃗ p χ 2 . Then we do the ⃗ p χ 2 integral in the rest frame of plasmon, defining all quantities in this frame with a * and redefining p * We simplify the delta function as where, and, Here, λ is the Källén function defined as λ(a, b, c) = a 2 + b 2 + c 2 − 2ab − 2ac − 2bc. Plugging this into eq. (B5) we get where, Here, we have used s = k 2 = ω 2 * (in CM frame, ⃗ k * = 0 =⇒ k 2 * = ω 2 * ). Subsequently, multiplying both sides of eq. (B1) by k µ k ν and contracting, we get and substituting from eq. (B9), we get: Putting together eqs. (B4) and (B12), we get for the inelastic Lenard's formula, eq.(B1), with C 1 given in eq. (B10)

Appendix C: Fixed Target: NA64 Events
We review the discussion in ref. [63] and begin by noting that in full generality, a 4-body phase space has 12 degrees of freedom. In particular, though, there are redundancies from invariance in rotation around the beam line, and from imposition of energy-momentum conservation, leaving us with 7 independent degrees of freedom. With this knowledge, the 4-body phase space can be written as where the kinematic quantities are, Here, q 2 = p 2 − p 4 , q 1 = p 1 − p 3 and the remaining momenta are as shown in fig. 5. The azimuthal angle of p 3 in the frame where ⃗ p 3 + ⃗ p χ 1 + ⃗ p χ 2 = 0 is given by ϕ R3χ 1 χ 2 3 and the solid angle between the DM particles is Ω Rχ 1 χ 2 χ . The Källén function denoted by λ is defined as The Jacobian of the transformation from E 4 , cosθ 4 to s 3χ 1 χ 2 , q 2 2 , required to connect between eq. (36) and eq. (C1), is given by The double differential cross section for the production corresponding to fig. 5 is where E 2 = E 0 is the incoming electron (beam) energy with its velocity given as The squared amplitude for the process shown in fig. 5 is with χ µν the DM emission piece of the amplitude, Here, the interaction operators are given by Γ µ (q) = −d χ σ µν q ν γ 5 for EDM and Γ µ (q) = iµ χ σ µν q ν for MDM. The term corresponding to electron scattering, with averaging over initial and sum over final spins is with, L ρσ,µν a = 1 2 Tr ( / p 4 + m e )γ µ ( / p 4 + / q + m e )γ ρ ( / p 2 + m e )γ σ ( / p 4 + / q + m e )γ ν The hadronic tensor describing the response of the nuclear target is Assuming Pb to be a scalar target gives W 1 (q 2 ) = 0, 1+t/d(A) , t = −q 2 1 , a(Z) = 111 Z 1/3 me and d(A) = 0.164 GeV 2 A −2/3 . The mass number and atomic number of the target nucleus are given by A and Z, respectively. We have neglected the magnetic form factor for Z >> 1.

3
where G n is the Gram determinant of dimension n and ∆ n is the Cayley determinant (symmetric Gram determinant) [44,109].
And finally, since the NA64 experiment isn't sensitive to the angular distribution of the outgoing DM particles, we can integrate over the solid angle between the DM particles, Ω Rχ 1 χ 2 χ in eq. (C1).
Appendix D: Derivation of the electron recoil rate formula for upscattered DM We follow the discussion in Appendix A of Essig et al [89] to re-derive the electron scattering rate for a general DM flux.
The integrated flux changes as follows in going from the Standard Halo Model DM distribution (g χ (v)) to a generic incoming DM flux per unit energy (dΦ/dK χ 2 ) Then, starting from eq. A12 of [89], we rewrite the cross section for a χ 2 to excite an electron from level 1 to level 2 of an atom as The rate of excitation events, for a given transition and given target electrons, is found by multiplying eq. (D2) by the incoming χ 2 flux. Using eq. (D1) we get this to be The rates for ionization of electrons bound in isolated atoms can be calculated with the simplifying assumptions of a spherically symmetric atomic potential and filled shells. The ionized electron can be treated as being in one of a continuum of positive-energy bound states, approximated to free particle states at asymptotically large radii. The ionization rate for an atom is found by taking eq. (D2), summing over occupied electron shells, and integrating over all possible final states. For ionization, with the final states being a continuum, the phase space is [89] ionized electron phase space= Here, l ′ , m ′ are the angular quantum numbers of the ionized electron final state and k ′ is its momentum at asymptotically large distances from the nucleus, with energy E R = k ′2 /2m e . Plugging this in, the ionization rate is given as with |f ion (k ′ , q)| 2 ≡ 2k ′3 (2π) 3 occ. states l ′ m ′ |f i→k ′ l ′ m ′ (q)| 2 .
We can write the electron recoil energy spectrum per detector mass per unit time as dR ion d∆E e = n Tσ e 64µ 2 χ 2 ,e n,l 1 ∆E e − E nl dK χ 2 dΦ dK χ 2 m χ 2 K χ 2 dq q|F DM (q)| 2 |f ion (k ′ , q)| 2 , where the total energy transferred to the electron is ∆E e = E R + E nl and n T is the number of targets per tonne. To explicitly show the order of integration, including the factor for depletion in numbers due to decay of χ 2 in travelling from the Sun to the Earth, and the energy dependent detector efficiency (ϵ(∆E e )), we get: dq q|F DM (q)| 2 |f nl→∆Ee−E nl (q)| 2 .

(D10)
Appendix E: Effect of accounting for detector resolution in recoil spectra We integrate over the theoretical differential event rate as given in eq. (23) to get the total number of events (shown by the color palette in fig. 7). For comparison with the experimental data, the theoretical spectra can be further smeared using a Gaussian distribution with energy-dependent width [13,125] where σ(E) is the recoil energy dependent energy resolution of the detector. In figures 13 and 14, we show the smeared spectra finding that in each case the signal+background rates still shows an excess over the background rates. We use the detector resolution from XENON1T with σ(E) = a. √ E + b.E and a = 0.310 ± 0.004 √ keV, b = 0.0037 ± 0.003 [13].  (23)).