Resonant Pseudo-Dirac Dark Matter as a Sub-GeV Thermal Target

Dark matter (DM) could be a pseudo-Dirac thermal relic with a small mass splitting that is coupled off-diagonally to a kinetically mixed dark photon. This model, particularly in the sub-GeV mass range, is a key benchmark for accelerator searches and direct detection experiments. Typically, the presence of even a tiny fraction of pseudo-Dirac DM in the excited state around the time of recombination would be excluded by DM annihilation bounds from the cosmic microwave background (CMB); thus, viable thermal histories must typically feature an exponential suppression of the excited state. We revisit assumptions about the thermal history in the resonant regime, where the dark photon mass is slightly more than twice the DM mass (to within $\sim10\%$), leading to an $s$-channel resonance in the annihilation cross section. This resonance substantially reduces the couplings required for achieving the observed relic abundance, implying that in much of the parameter space, the DM kinetically decouples from the Standard Model well before the final DM relic abundance is achieved. We find that the excited state is not thermally depopulated in this regime. In spite of this, we find that the presence of the excited state does $\textit{not}$ violate CMB bounds, even for arbitrarily small mass splittings. The present-day abundance of the excited state opens up the possibility of signatures that are usually not relevant for pseudo-Dirac DM, including indirect detection, direct detection, and self-interacting DM signatures.


I. INTRODUCTION
The origin of dark matter (DM) remains an elusive mystery.If the DM thermalizes with the Standard Model (SM) plasma in the early Universe, then thermal freezeout provides a compelling explanation for the observed abundance of DM.In particular, thermal freeze-out is relatively insensitive to the initial conditions of the early Universe and the relevant couplings can be probed in a number of ways using direct detection, indirect detection, and collider observables [1].If the DM is lighter than the ∼GeV scale, then the Lee-Weinberg bound [2] implies that the mediator for DM-SM interactions cannot be a SM force carrier, which opens up the possibility of a "dark sector," with auxiliary forces and matter fields beyond just DM.A simple, technically natural example of a new mediator is a dark photon that kinetically mixes with the SM photon [3][4][5][6][7][8][9][10][11][12].Given the null detection of weak-scale DM thus far (see e.g.Refs.[13,14]), lighter dark sectors are of increasing interest to the community (see e.g.Refs.[15][16][17][18][19]) and there are a range of new and proposed experimental methodologies that will be sensitive to these DM candidates [20][21][22][23][24][25][26].
The standard thermal freeze-out mechanism for sub-GeV DM is subject to strong bounds from cosmic microwave background (CMB) anisotropies assuming annihilation through an s-wave process to visible SM particles.Even after freeze-out, DM can still annihilate at a sub-Hubble rate and inject considerable energy into the SM plasma near the time of recombination.Energy injection in the form of visible particles would observably modify the properties of the plasma even if the DM annihilations are extremely rare, since a ∼part-per-billion fraction of the DM annihilating would be enough energy injection to ionize all the atoms in the Universe.Considering the effects on CMB anisotropies as measured by Planck, current bounds on DM annihilation rule out swave thermal freeze-out of DM below ∼ 10 GeV, with the exact value depending on the SM final state [27,28].
A well-studied way to bring s-wave freeze-out into consistency with CMB constraints is by introducing a small mass splitting between non-degenerate DM states [29][30][31][32][33][34][35].There is no symmetry that prevents Dirac fermions from splitting into two Majorana mass states, and this is easily realized in models where DM is charged under some new dark gauge symmetry at high energies which is broken at low energies [36,37].The dark matter multiplet, consisting of χ 1 and χ 2 , acquires a mass splitting δ = m χ2 − m χ1 , which can be naturally small if the dark symmetry is approximate (for example, the small neutron-proton mass splitting is protected by an approximate isospin symmetry).In this work, we do not specify the origin of the mass splitting and treat it phenomenologically.However, we note that small mass splittings are generally easy to accommodate from a modelbuilding perspective in situations with a small overall mass scale and small couplings [29,[38][39][40][41].In this model, the couplings of DM with the dark photon are purely off-diagonal, i.e. the only vertex with the dark photon couples χ 1 and χ 2 .For this reason, annihilation rates to SM final states can be significantly reduced because the leading-order tree-level annihilation process requires a large χ 2 abundance, which may be thermally depleted like e −δ/T at temperatures T ≲ δ [41,42].Given the ∼eV-scale temperatures that are relevant for recombination, mass splittings with δ ≳ 1 eV can be compatible with CMB constraints in parts of the parameter space (with larger mass splittings being unconstrained for a wider range of couplings and DM masses).
Alternatively, if the DM annihilation occurs close to a pole in the cross-section (for instance when the mediator is close to twice the DM mass), then the relevant couplings to achieve the observed relic abundance can be lowered substantially [43,44].If this pole is relevant for setting the DM abundance at early times but not during the recombination epoch, then the CMB bounds are relaxed because of the lower off-resonance annihilation rate to SM final states.The CMB bounds were carefully studied in Ref. [45] for the case of Dirac DM interacting with a dark photon with m A ′ ≈ 2m χ .Meanwhile, the resonant regime for inelastic pseudo-Dirac DM has been studied primarily in the context of its signature at colliders given the modification to the predicted couplings [38,[46][47][48].The cosmology of resonant pseudo-Dirac DM, on the other hand, has yet to be studied in detail.In particular, the substantial effects of early kinetic decoupling of the DM have been overlooked so far.
In this work, we perform a comprehensive study of the cosmology of pseudo-Dirac DM in the resonant regime.We find that even in the mildly resonant regime with (m A ′ −2m χ )/m A ′ ∼ 10%, pseudo-Dirac DM can have arbitrarily low mass splittings without violating limits from the CMB.Moreover, we find that in most of the parameter space the excited state has a high relic fraction.This provides a strong contrast to most pseudo-Dirac thermal histories which feature an exponential suppression of the excited state due to thermal depletion.Accordingly, the cosmology and astrophysics of this DM candidate are quite different from the usual pseudo-Dirac parameter space, as are the direct and indirect DM detection signatures.
The rest of this article is organized as follows.In Section II, we review the model and the relevant processes that affect the cosmology of this DM candidate in the early Universe.In particular, we solve the Boltzmann equations for the density and temperature evolution of the DM states χ 1 and χ 2 .In Section III, we consider cosmological and astrophysical signatures including Big Bang Nucleosynthesis, the CMB, self-interacting DM (SIDM), and indirect detection.In Section IV we discuss prospects for detecting this DM candidate using terrestrial experimental methods.Discussion and concluding remarks follow in Section V

A. Pseudo-Dirac DM Parameter space
We consider a light (m χ ≲ 10 GeV) pseudo-Dirac DM model with its relic abundance set by annihilation to SM final states via a dark photon mediator.We focus on this mass range primarily because Dirac DM with m χ ≲ 10 GeV is excluded by the CMB for s-wave freeze-out to visible final states [28].In this model, the interaction terms are where χ 2,1 are the excited and ground states, respectively, that couple with interaction strength g χ to a vector mediator A ′ that kinetically mixes with the SM photon with mixing parameter κ.We focus on parameter space where the mass splitting is much smaller than the DM mass, δ = m χ2 −m χ1 ≪ m χ1 .In the following discussion, we denote the average mass of the two states as m χ .We are furthermore interested in the resonant regime where m A ′ ≈ m χ2 + m χ1 = 2m χ1 + δ = 2m χ .We parameterize the proximity to resonance with the parameter where s 0 = (m χ1 + m χ2 ) 2 .In this work, we consider ϵ R ∈ [0.001, 0.1].The lower limit is motivated by photodisassociation bounds coming from BBN as discussed in Sec.III A. On the other hand, as ϵ R > 0.1, we approach the non-resonant limit.We remain agnostic about the mechanism responsible for generating the mass splitting as well as the resonance, however the parameters we consider are self-consistent even in a minimal UV setup.
For instance, we could consider a complex dark Higgs with a dark charge of 2 interacting with a Dirac fermion with a mass m D [36,37,49].The breaking of the dark U (1) symmetry through the vev of the dark Higgs, v D , then results in both the mass term for the dark photon as well as a Majorana mass term for the Dirac fermion which generates the mass splitting.If the dark Higgs is heavy with v D ≫ m D then it will not participate in the dynamics that determine the DM thermal history.With this hierarchy in mind, having the dark photon mass near its resonant value pushes g χ ∼ m D /v D ≪ 1, which we show below is consistent with setting the observed DM relic abundance in the resonant regime.Finally, the mass splitting is determined by the Yukawa coupling of the Dirac fermion y χ , with δ ∼ y χ v D .
Though the dark photon is the mediator of this model, it can be resonantly produced on shell via inverse decays.The dark photon can subsequently decay invisibly, or visibly to SM final states, where R(m A ′ ) is the empirically determined branching ratio of σ(e + e − → hadrons)/σ(e + e − → µ + µ − ) [50,51] at centre-of-mass energy The total on-shell decay width is Γ A ′ ≡ Γ DM + Γ SM ≡ Γ DM + Γ e + e − /B e , where B e is the branching ratio of the dark photon to electrons.T dec [GeV] The SM temperature as a function of the dark photon mass at which DM-SM scattering (dashed), DM-DM scattering (dotted) and DM-DM annihilation (dot-dashed) decouple for ϵR = 0.001 (light green) and ϵR = 0.1 (dark green), with gχ = 0.01.The solid lines show when annihilation becomes resonant, with T ∼ ϵRm A ′ .For each parameter point in this plot κ is chosen so as to obtain the observed DM relic abundance ignoring the effects of early kinetic decoupling.The shaded areas between dotted lines correspond to varying δ between 1 − 100 eV, from bottom to top.In some parts of the parameter space, resonant depletion of the DM abundance occurs much later than other decoupling processes, indicating that the early decoupling can influence the subsequent relic DM abundance.

B. Relic Density Assuming Thermal Equilibrium
In our parameter region of interest, DM obtains its relic abundance when s-channel processes like χ 1 χ 2 → SM SM become inactive, usually well after chemical freeze-out for ϵ R ≪ 1.Using the formalism developed in Ref. [47], the corresponding thermally averaged crosssection can be written as where x = m χ1 /T , ϵ = (s − s 0 )/s 0 is a dimensionless measure of the kinetic energy and with, The thermally averaged cross-section can be further reduced to semi-analytic forms in the non-relativistic (ϵ ≪ 1) and resonant (ϵ ≈ ϵ R ) limits [45].
In order to calculate the total DM relic density, one needs to solve the Boltzmann Equation for χ 1 χ 2 → SM SM , where Y tot = Y χ1 + Y χ2 is the total comoving density for DM with Y χ1,χ2 = n χ1,χ2 /s, s is the entropy density, and the effective thermally averaged cross-section is [43,47] ⟨σv⟩ eff = 2(1 + δ/m χ ) 3/2 e −xδ/mχ 1 (1 + (1 + δ/m χ1 ) 3/2 e −xδ/mχ 1 ) 2 ⟨σv⟩.(10) The total DM relic density is then given by where g * correspond to the effective relativistic degrees of freedom in the early universe.It can be seen from Eqs. ( 6) and ( 7) that the cross-section is resonantly enhanced when χ 1 and χ 2 have enough energy to produce the A ′ on shell, leading to very efficient annihilation at a temperature, T ∼ ϵ R m χ1 .Due to the resonant enhancement, as ϵ R → 0 even very tiny couplings are able to efficiently deplete the DM in the early Universe to obtain the observed DM relic density.

C. Early Kinetic Decoupling
A crucial caveat to the relic density calculation detailed above is that it assumes that the DM and SM remain in kinetic equilibrium while annihilations are active (including after chemical freeze-out) through scattering processes, primarily off of electrons χ 1 e ↔ χ 2 e.Since this is a t-channel process, it does not benefit from the same resonant enhancement as the s-channel annihilations.Therefore, when the couplings between the two sectors are small, the assumption of kinetic equilibrium may no longer hold and DM can kinetically decouple well before DM annihilation χ 1 χ 2 → SM SM hits its resonance (compare e.g. the dashed and solid lines in Fig. 1).As a result, to accurately calculate the DM relic density, one needs to solve a coupled system of differential equations tracking the evolution of both the DM number density as well as the dark sector temperature, y(x) ≡ m χ1 T DM s −2/3 [52][53][54], Dashed lines indicate CMB annihilation constraints on those couplings, while solid lines are consistent with the CMB.In this parameter space, the smaller of the two couplings (corresponding to whichever decay channel is the bottleneck) determines the relic abundance, with the larger coupling being irrelevant.Small deviations from this behaviour occur at large values of mχ and ϵR where early kinetic decoupling has a significant effect on the relic abundance.These trends contrast with the parameter space for standard thermal freeze-out where the product of couplings is the most relevant for setting the relic abundance.For comparison, we show the couplings for thermal freeze-out for m A ′ = 3mχ (off resonance) as dot-dashed lines.
where H is the normalised Hubble rate as defined in Ref. [53] and the subscript "neq" denotes that the corresponding thermal average is over the DM phase distribution at a temperature T DM distinct from the SM temperature T , i.e., assuming DM is not necessarily in kinetic equilibrium with the SM.Additionally, ⟨σv⟩ 2 is a temperature-weighted analog of the usual thermally averaged annihilation cross-section ⟨σv⟩, as defined in Ref. [53], and γ(T ) is the DM-SM momentum transfer rate which is a measure of DM-SM elastic scattering.
Terms involving scattering and annihilation can both keep the DM temperature coupled to the SM.These Boltzmann equations have been extensively studied for the case of elastically decoupling relics [52][53][54][55].For inelastic DM models, the corresponding Boltzmann equations may have an additional dependence on the masssplitting δ, which would appear in various cross sections, and there could also in principle be separate thermal evolution of the ground and excited state species.However, in our parameter region of interest, δ ≪ m χ and δ is also much smaller than the decoupling temperatures for all relevant processes.Therefore, we explicitly find that the thermal history of this model behaves, to a very good approximation, as a strictly Dirac model during the temperatures relevant for freeze-out.Therefore, to calculate the relic pseudo-Dirac DM density, we use the publicly available Boltzmann solver DRAKE [54] modified for a resonant Dirac DM model.

D. Parameter Space
The couplings required to reproduce the observed DM abundance are shown in Fig. 2 for different values of m χ and ϵ R .Note that δ is generally much smaller than all the temperatures that are relevant for setting the relic abundance, and therefore it is irrelevant in determining the couplings.As expected, for a given DM mass, smaller values of ϵ R (indicated by thinner lines) result in a larger resonant enhancement in the annihilation crosssection, and corresponds to smaller couplings reproducing the relic density, thereby shifting the lines downwards and to the left.For a fixed value of m χ and ϵ R , the shape of the curve in the g χ − κ plane can be explained by considering the thermally averaged cross-section in the limit ϵ R ≪ 1.In this case, as was pointed out in Ref. [45], the slower decay (Γ SM vs. Γ DM ) is the bottleneck in terms of determining the final DM abundance, This implies that for g χ ≪ κ (κ ≪ g χ ), the relic density becomes independent of κ (g χ ) resulting in the asymptotic behavior seen in Fig. 2.
In the limit ϵ R → 0, we find that the relic density obtained using the coupled system of Boltzmann Eqs. ( 12)-( 13) differs from the standard Boltzmann treatment (i.e.assuming identical SM and DM temperatures) by at most a factor of ∼ 2. For ϵ R ∼ 1, one would naively expect the difference to be even smaller since we are not only further off-resonance but are also pushed toward larger couplings where the expectation is that kinetic equilibrium should be maintained more easily.However, we find that the difference between the two treatments can be as large as an order of magnitude in the relic density in the region ϵ R ∼ 0.1.This can be attributed to the deviation of the DM temperature from the SM temperature and was earlier discussed in the context of scalar singlet DM in Ref. [56].In particular, for ϵ R ≪ 0.1, the deviation is very small and positive, T DM ≳ T whereas for ϵ R ∼ 0.1 the deviation is large and negative, T DM ≪ T .Although one can numerically estimate the sign and magnitude of this deviation by studying the interplay of the various terms on the right hand side of Eq. ( 13), they can also be understood qualitatively by considering the underlying DM phase space during and after chemical freeze-out.
Under the assumption that DM has already kinetically decoupled, the final DM density can be assumed to be proportional to the annihilation cross-section averaged over the DM temperature, ⟨σv⟩ neq (in analogy with Eq. ( 11)) In general, as ϵ R → 0, the DM particles need only very little momentum to hit the resonance and annihilate efficiently, implying that resonant annihilation depletes the low-momentum tail of the DM distribution and shifts the average DM momentum (and therefore the temperature) to larger values.During chemical freeze-out, since the average DM momentum is already large, only a small fraction of DM particles can annihilate resonantly.The increase in the DM temperature further decreases the available phase space for resonant annihilation and therefore decreases ⟨σv⟩ neq .This results in the small dip in ⟨σv⟩ neq for ϵ R = 0.001 around chemical freeze-out x f ∼ 20 seen in Fig. 3 and corresponds to reducing the efficiency of DM annihilations and increasing its abundance.After DM has chemically decoupled, its temperature now redshifts as matter and therefore DM cools much faster than the SM, increasing the relative number of DM particles with low momentum.As a result, resonant annihilation which is active at T DM ∼ ϵ R m χ happens at slightly earlier times (since T DM ≪ T ) compared to the kinetically coupled case when they occur at T ∼ ϵ R m χ .Hence, resonant annihilation is more efficient in depleting the DM (due to the x dependence of Eq. ( 15)).Since these two effects change the relic density in opposing ways, the final relic density is only slightly different from the kinetically coupled case.This is especially true given that the 1/x 2 weighting in the integral of Eq. ( 15) ensures that the most substantial contributions to the relic abundance come from early times before the difference between ⟨σv⟩ eq and ⟨σv⟩ neq becomes too large.3. The thermally averaged annihilation cross-section for kinetically decoupled (solid) and kinetically coupled (dashed) DM as a function of x (defined with respect to the SM temperature, T ) for ϵR = 0.1 (dark green) and ϵR = 0.001 (light green) and mχ = 100 MeV.We fix gχ = 0.01 and choose κ such that we reproduce the observed DM abundance after solving Eqs. ( 12)-( 13).Early kinetic decoupling can suppress or enhance resonant annihilation at given temperature.
For ϵ R ∼ 0.1, on the other hand, DM particles need larger momentum to annihilate resonantly and therefore resonant annihilations shift the average DM momentum (and temperature) to smaller values.Additionally, resonant annihilation is active exactly during chemical freezeout, x DM ∼ ϵ −1 R ∼ x f .This implies that if DM is kinetically decoupled, the depletion in the large-momentum DM phase space effectively turns off resonant annihilations quite quickly as shown by the dark green line in Fig. 3. Furthermore, the x 2 scaling of Eq. ( 15) in this regime enhances the difference in the total relic abundance since ⟨σv⟩ eq and ⟨σv⟩ neq differ substantially around the time of chemical freeze-out.As a result, the relic density is significantly altered: DM is over-produced and larger couplings are required to obtain the observed DM abundance.

E. Late-time Abundance of the Excited State
The relative fraction of excited-state particles, f * = n χ2 /n χ1 is a key quantity in determining the impact of late-time DM behaviour in cosmological environments as well as terrestrial experiments, as discussed further in Sections III and IV.Even if DM is symmetrically produced in the ground and excited states, the ground and excited states can inter-convert through processes within the dark sector or through scattering processes with the SM as long as their rates exceed the Hubble expansion rate.In particular, as long as chemical equilibrium is maintained within the dark sector, the excited state number density is given by n χ2 ∼ n χ1 e −δ/TDM .Chemical equilibrium in the dark sector can be maintained through DM-SM scattering, which also maintains kinetic equilibrium with the SM, χ 1 e ↔ χ 2 e, or through DM up/downscattering, χ 1 χ 1 ↔ χ 2 χ 2 .The fractional abundance of a cosmologically stable excited state at late times is determined the DM temperature when it chemically decouples, f * = n χ2 /n χ1 ≈ e −δ/T chem , where T chem is determined by whichever of the processes listed above decouples last, T chem = min[T χe , T χχ ], where Here, n e is the electron number density, n χ2 = n χ1 e −δ/TDM is obtained by scaling back the present-day DM abundance n χ1 ∼ T eq T 3 /m χ1 , and the relevant cross-sections are from Ref. [41].If χ 2 χ 2 ↔ χ 1 χ 1 decouples after χ 2 e ↔ χ 1 e, two temperature scales enter in Eq. ( 17), the SM temperature T that largely determines the Hubble rate, and the DM temperature which is a function of the temperature at which DM decouples from the SM, T DM ∼ T 2 /T χe .
In the standard thermal history for inelastic DM, one finds that T chem ≲ δ owing to the large DM-DM and/or DM-SM couplings which ensures chemical equilibrium in the dark sector is maintained until late times.This results in a strong suppression of the excited state at late times, f * ∼ O(10 −4 ) [42,57].However, the small cou-plings present in our parameter space result in T chem ≫ δ, and therefore a similar abundance of the ground and excited state, f * ∼ 1 in much of the parameter space.This is represented in Fig. 4, in which we show the relative abundance of the ground and excited states as a function of the DM mass for different values of ϵ R and g χ .Note that f * ∼ 1 for g χ ≪ 0.1 for all DM masses of interest.For g χ ∼ 0.1, we find a suppression of the excited state to f * ∼ 0.01 when m χ ≲ 100 MeV.

III. COSMOLOGICAL AND ASTROPHYSICAL CONSTRAINTS
Due to the high late-time abundance of the excited state, this thermal history has unique signatures compared to other pseudo-Dirac DM thermal histories.In this Section, we determine the qualitatively new behavior in astrophysical systems caused by the presence of the excited state and estimate the resulting constraints on the model as inferred from existing measurements.

A. BBN
Sub-GeV DM may significantly impact the abundance of light elements in the Universe produced during BBN.If DM has a thermal abundance and is relativistic at a SM temperature of a few MeV, it contributes to the number of relativistic degrees of freedom in the early universe modifying N eff and changing the abundances of light elements like helium.Measurements of the helium fraction can thus be used to place a lower bound on the mass of thermal DM, m χ ≳ 10 MeV [58,59].In this work, we conservatively consider DM above this scale to ensure that the parameter space is not ruled out.However, we note that this constraint could be slightly weaker in parts of our parameter space due to the early kinetic decoupling of DM from the SM.In particular, if the DM decouples early enough, then some of the SM degrees of freedom in the bath at that time, given by g * ,0 in total, heat the SM bath at later times and raise the SM temperature relative to T DM by a factor of a few in order to conserve entropy, T SM /T DM ∼ (g * ,0 /g * ,BBN ) 1/3 .This means that the DM contribution to the energy density, and hence N eff , could be diluted by a factor of (g * ,0 /g * ,BBN ) 4/3 .Furthermore, once the DM becomes non-relativistic at a temperature T DM ∼ m χ , the DM temperature will drop even further relative to the SM temperature, T DM ∼ T 2 SM (g * ,0 /g * ,BBN ) −2/3 /m χ for T SM < m χ .Therefore, by the time of BBN, MeV-scale DM may have had its energy density diluted and may be non-relativistic, thus not contributing substantially to N eff .In concert, these effects could increase the range of allowed masses for this model; we leave a more detailed exploration to future work.
In addition to the modification of N eff , DM annihilating into SM states at a temperature of a few keV can inject energy into the SM plasma causing the photo-disassociation of light nuclei like deuterium.The corresponding bound on the annihilation cross-section not only depends on the DM mass and the relevant final states but also depends on the temperature of the kinetic decoupling, T kd , in the case when the thermally averaged cross-section has a temperature dependence [58].In our model, the cross-sections show this dependence around the keV-scale temperatures relevant for photo-disassociation for ϵ R ≪ 1. Accurately evaluating this bound for such small ϵ R 's is therefore non-trivial, and we leave a detailed study for future work.For the purposes of this work, we note that for a velocity-independent annihilation cross-section, the bound from photo-disassociation corresponds to σv ≲ few × 10 −25 cm 3 /s, and the constraint for velocity dependent cross-section gets weaker with increasing T kd [58].For resonant annihilations, even the small couplings considered in this work can result in significantly larger cross sections at the late times relevant to BBN.Since the resonant annihilation cross section peaks near T ∼ ϵ R m χ , in this work we consider ϵ R ≥ 0.001, which along with the conservative lower bound on m χ discussed above, ensures that the DM annihilation cross section is always below the upper bound during temperatures relevant for photo-disassociation.

B. Cosmic Microwave Background
DM annihilation into SM final states during recombination can inject energy into the SM plasma.This energy injection can alter the ionization history and affect CMB anisotropies due to the scattering of CMB photons on additional free electrons that would have otherwise recombined into neutral atoms in the standard cosmology.This injection is usually described in terms of an effective parameter [28], where f * is the fraction of DM in the excited state as described in the previous Section, f eff is the efficiency fraction of injected energy that gets deposited in the plasma (where we adopt the values from Ref. [45] using the spectra from Refs.[60,61] and which depend on the DM mass and SM final state), and ⟨σv⟩ CMB is the total annihilation cross section into SM states at recombination.The bound on p ann includes annihilation channels into all visible SM final states (i.e.excluding neutrino final states).In our case, the final states are dominantly electrons (and muons for DM masses above the muon mass), and we include the relevant branching fractions as appropriate when computing the total annihilation cross section.In the case of multiple possible final states, we weight the cross sections by the relevant deposited energy efficiency fraction, f eff .Since the DM particles are highly non-relativistic during and after recombination, to obtain ⟨σv⟩ CMB we evaluate Eq. ( 6) in the limit ϵ → 0, To compute the CMB limit on our parameter space, we require that the total f eff -weighted cross section in Eq. ( 19) not exceed the one implied by the bound on p ann .The excluded parameters are depicted in Fig. 2 as dotted lines which, despite being excluded by the CMB, would yield the observed amount of DM, as described in the previous Section.We find that -in contrast to the case of non-resonant pseudo-Dirac DM -very small mass splittings and large excited state fractions remain unconstrained by the CMB, particularly in the part of the plane corresponding to small couplings.Additionally, for the sub-keV values of δ that we consider here, the bound from the CMB is independent of δ because we do not get a substantial enough suppression in the excited state abundance as a result of the early kinetic decoupling (i.e.f * ∼ 1 near the boundary of the CMB exclusion).

C. Self-interacting DM
Models of inelastic DM can also have unique SIDM behaviour in DM halos, affecting density profiles and subhalo mass functions in a way that is distinct from purely elastic SIDM [62,63].In this model, elastic scattering between the ground and excited states can occur at tree level, which is especially relevant in the thermal histories we consider where it is possible to have a high abundance of the excited state at late times.As we are in the Born regime, α χ m χ /m A ′ ≈ α χ /2 ≪ 1, we can use results previously derived in the literature for the relevant SIDM cross sections.The s-channel resonance is not typically relevant in astrophysical environments (in contrast to e.g.Ref. [64]), since the lowest value of ϵ R we consider is ∼ 10 −3 due to the strong BBN constraints described in Sec.III A, which corresponds to a minimum resonant velocity of ∼ 10 4 km/s, in contrast to the ∼ 10 2 km/s velocities that are typical in galaxies like the Milky Way.
For tree-level elastic scattering between the ground and excited states, we employ the Born cross section of Ref. [65], assuming that the t-channel dominates.We find that the elastic scattering cross section is generally much less than the σ/m χ ∼ 1 cm 2 /g characteristic of SIDM constraints from merging galaxy clusters [66].The exception to this lies at the edge of perturbativity g χ ∼ 1 for the lightest DM masses we consider, m χ ∼ 10 MeV.However, this part of the parameter space is excluded by the CMB, as evident from Fig. 2.Even if the CMB constraint could be circumvented, in this part of the parameter space the abundance of excited state particles (which would be required for tree-level elastic scattering) is generally more suppressed as shown in Fig. 4, which would also weaken the bounds from merging clusters.Elastic scattering between particles in the same mass state (i.e. two ground state particles) only occurs at the 1-loop level in this model; we confirm that the cross sections for these processes fall well below the ∼ 1 cm 2 /g level (due to the α 4 χ scaling) using the expressions in the Supplemental Materials of Ref. [57].We therefore conclude that merging cluster constraints on elastic SIDM do not have any impact in the viable parameter space for this model.
Upscattering from the ground state to the excited state is kinematically forbidden below the velocity threshold of v ∼ δ/m χ , which can take on a wide range of values in the parameter space we consider, some of which are relevant to astrophysical systems.Above this velocity threshold, the cross section for upscattering in the Born regime saturates to the value of the elastic cross section between the ground and excited states [67]; while this cross section is generally below ∼ 1 cm 2 /g in our parameter space, the endothermic nature of the scattering can lead to substantial qualitative differences in the DM distribution compared to the elastic scattering case [63], and therefore strong conclusions cannot be drawn about observational prospects without dedicated simulation work.Similarly, downscattering from the excited state to the ground state takes on the same value as upscattering once the velocities are above threshold; below threshold, there is an enhancement to the downscattering cross section, rendering it potentially quite large at low velocities (with σ/m χ ≫ 1 cm 2 /g).The effects of inelastic scattering has only begun to be explored in simulation, with no direct analog of this situation having been analyzed.In this thermal history, up to 50% of the DM begins in the excited state similar to Ref. [62], which showed that even a few percent of the DM downscattering can have a significant effect on the structure of DM halos.On the other hand, the velocity thresholds in our parameter space can be easily accessible in astrophysical systems, potentially leading to some upscattering, which has a highly nontrivial interplay with the effects of downscattering as studied in Ref. [63].Clearly futher exploration of the parameter space in simulation will be fruitful for connecting lateuniverse halo observables to DM parameters motivated by self-consistent DM thermal histories.

D. Indirect Detection
DM annihilation in astrophysical environments produces cosmic rays and high-energy photons, which can be employed to search indirectly for DM.A null detection of annihilation byproducts can thus be used to constrain the annihilation cross section, ⟨σv⟩.The strength of the expected signal additionally depends on the integrated DM density (the J factor).As discussed in the previous Subsection, SIDM effects, particularly those from inelastic interactions, may alter the expected density profiles of DM halos.As this has yet to be quantified via simulation, we assume in this discussion that resonant inelastic DM has the same J factor as what is typically considered in the literature.Additionally, as shown in Sec.III, the CMB constraints on DM annihilation push us toward g χ ≲ 0.1, which corresponds to the ground and excited states being equally populated, or f * ∼ 1, for the δ ranges of interest (see also Fig. 4).Therefore, as opposed to standard thermal histories of inelastic DM, we do not necessarily have a late-time suppression of annihilation from the thermal depletion of the excited state.
Since the DM velocity in present-day halos is too small (by about two orders of magnitude) for the dark photon mediator to be produced on shell, the DM annihilation rate at late times can be calculated in the heavy mediator (non-resonant) limit.In Fig. 5, we show the constraints on the total present-day DM annihilation crosssection to SM states, ⟨σv⟩ 0 , from measurements of the gamma-ray flux from the galactic center using Fermi IN-TEGRAL, COMPTEL, and EGRET [61].We also show constraints from Voyager [68] and XMM-Newton [69], which bound the DM annihilation cross-section to electrons and therefore have been cut off at the muon threshold, m µ ∼ 100 MeV beyond which annihilation to other final states becomes relevant.For reference, we show the resonant pseudo-Dirac DM parameter region allowed by the CMB in grey.Even though current constraints do not yet probe the target parameter space that is allowed by the CMB, future telescopes such as GECCO [70,71], MAST [72], GRAMS [73,74] and AMEGO [75,76] will be able to explore the parameter space [91].e + e − LHC π 0 → Xγ e-brem FIG. 6.Some of the strongest constraints in the κ − m A ′ plane from accelerator experiments, including BaBar [77][78][79], NA64 [80,81], LHCb [82,83], CMS [84], NuCal [85,86], E141 [87], NA48 [88], and E137 [89], as computed with DarkCast [90].The different colours correspond to different dark photon production channels.We show two representative cases where the dark photon decays primarily invisibly (left) and visibly (right).The solid black lines represent the target for ϵR = 0.001 and ϵR = 0.1, with the dashed portions corresponding to exclusions from the CMB.Also shown for comparison is the thermal prediction assuming thermal freeze-out with m A ′ = 3mχ (dotted line) [48].

IV. TERRESTRIAL SEARCHES A. Accelerator searches
Accelerator experiments provide a complementary probe to search for DM and any associated mediators.Dark photons with masses in the MeV-GeV range can be produced at collider and beam dump experiments, resulting in either missing transverse energy or displaced vertex signatures.The production cross-section for the dark photons depends on their coupling to the SM, κ, whereas the decay signature depends on their lifetime, τ A ′ = 1/Γ A ′ which is in general a function of κ, g χ , ϵ R and δ.Since δ ≪ m A ′ , m χ , one can assume to a very good approximation that Γ A ′ is independent of δ (i.e.we can take δ → 0).In this limit, the bounds presented in Ref. [45] which were obtained using a modified version of Darkcast [90] can be directly applied to our parameter space.In particular, as was pointed out in Ref. [45], for a given m A ′ , the constraints on κ depend only on a combination of g χ and ϵ R through the dark photon's reduced invisible decay width, In Fig. 6, we display the bounds on dark photons in the κ−m A ′ plane for two fixed values of the dark photon's reduced invisible width, γ inv = 10 −5 (left) and γ inv = 10 −13 (right).The shaded regions correspond to different dark photon production channels.For γ inv = 10 −5 , the invisible decay width is larger than the visible one for κ ≲ 0.1 and therefore the beam dump experiments that search for A ′ decays to leptons lose sensitivity.The solid black lines in both panels correspond to the κ values that reproduce the observed DM abundance following the production mechanism outlined in Sec.II, for ϵ R ∈ [0.001, 0.1].The dashed parts of the black lines correspond to the points excluded by the CMB.We note that for Γ DM > Γ SM , (or equivalently g χ ≫ κ), the CMB constrains the part of the parameter space which cannot be probed by accelerator experiments making the accelerator and cosmological bounds highly complementary.Additionally, we find that for a given γ inv and ϵ R , or equivalently, a fixed g χ , there exists a maximum m A ′ beyond which DM is always over-produced.From Eq. ( 14), we see that the DM relic density is proportional to m A ′ / min(g 2 χ , κ 2 ) meaning that once the couplings are fixed, a larger m A ′ results in a larger DM abundance.This causes the relic lines to be vertical in the right panel of Fig. 6 (in the left panel, the maximum m A ′ lies outside the plotted region).Finally, we note that for a given γ inv , κ-values bounded by the two curves corresponding to ϵ R = 0.001 and ϵ R = 0.1 also reproduce the observed DM abundance for different values of ϵ R resulting in a much broader thermal target.For reference, we show the usual thermal target assuming non-resonant, s-wave thermal freeze-out (dotted black line) in the left panel of Fig. 6.Future accelerator experiments are poised to more fully explore the parameter space of this model, as shown in Fig.   6 except with some of the strongest projections for future experiments including Belle II [92], FASER [93], HPS [94], LDMX [95], LHCb [96,97], SeaQuest [98], SHiP [99,100], and Yemilab [101].

B. Direct Detection
Most thermal histories for pseudo-Dirac DM result in a relic excited state fraction suppressed by several orders of magnitude such that most of the DM in the halo of the Milky Way (MW) is in the ground state.In this case, because of the off-diagonal coupling, the only tree-level scattering process on a SM target would be upscattering to the excited state.In contrast, for the thermal history considered in this work, around half of the DM is in the excited state at late times for most parts of the parameter space.In particular, the largest values of g χ that suppress the abundance of χ 2 in this thermal history are already excluded by the CMB (see Figs. 2 and 4).Therefore, the primary signature of pseudo-Dirac DM in direct detection experiments is downscattering of the excited state χ 2 e → χ 1 e, which deposits an energy ∼ δ.The deposited energy can ionize electrons in the target which can be detected either directly through charge-coupled devices [102], or by detecting secondary scintillation photons using photomultiplier tubes [103].The absence of a kinematic barrier for downscattering implies an enhancement in the event rate for sub-GeV DM.
In the following discussion, we consider DM-electron scattering in semiconductor and Xenon targets, and we place constraints on the fiducial DM-electron scattering cross-section [20], where µ χ,e is the reduced mass of the DM-electron system.The recoil energy of the electron is E Re = ∆E e − ∆E B where ∆E B is the electron binding energy and the energy deposited by sub-GeV DM downscattering is for momentum transfer q and relative velocity v.For downscattering, the minimum velocity required to transfer a momentum with magnitude q and an energy ∆E e to the electron is therefore given by v min (q, ∆E e ) ≡ ∆E e − δ q which corresponds to the differential event rate for atomic targets [20,104], χe n,l where ∆E n,l is the binding energy of the electron in the nl shell, ρ χ2 = f * ρ DM ≈ 0.5 × 0.4 GeV cm −3 is the density of the excited state [51], and f nl→∆Ee−E nl (q) is the electron ionization form factor which we evaluated numerically following the prescription of Ref. [105] using DarkART [106].The function η(v min ) in Eq. ( 24) can be related to the DM-velocity distribution f χ (v) by In this work, we assume the Standard Halo Model for f χ (v) [107].The magnitude of the momentum transferred to the electron for a given DM mass and energy deposition ∆E e is bounded by where v max is the largest DM velocity relative to the detector frame as determined by the escape velocity of the halo at the Earth's position and the velocity of Earth in the Galactic rest frame [107].A similar expression as Eq. ( 24) can be obtained for semiconductor targets which depends instead on the crystal form factor as described in Ref. [21].The relevant energy threshold for such a target is the band gap between the valence band to the conduction band.Note that in Eq. ( 24), we have set the DM form factor to unity, F DM = 1, which corresponds to the heavy mediator limit, m A ′ ≫ αm e as determined by the typical Fermi momentum (for semiconductors targets) or inverse Bohr radius (for atomic targets).
In contrast to the case of elastic scattering, the minimum momentum transfer, q min for a given energy transfer ∆E e , can be zero if the mass splitting is above the energy threshold for ionization.This results in much larger event rates since the inelastic scattering kinematics have substantial overlap with peaks in the electron ionisation form factor [108].A similar effect also occurs for semiconductor targets where the peaks in the crystal form factor become more kinematically accessible for δ > 0 [21].Inelastic scattering also results in a characteristic spectrum of events peaked around ∆E e ∼ δ.As a result, bounds on σ e derived under the assumption of elastic scattering cannot be directly applied to the parameter space of this model.
In order to re-cast the bounds, we use the prescription outlined in Ref. [109] for calculating event rates in the XENON10 [109,110], XENON1T [103] and SEN-SEI [111] experiments.For XENON10 and XENON1T, we calculate the electron ionization form factors using DarkART [106].The crystal form factors for SENSEI are obtained from QEDark [21].We use the publicly available data for the three experiments [103,[109][110][111] to derive 90% confidence level exclusions on σ e , shown in Fig. 8.For XENON10 and XENON1T, we show the exclusion for δ = 10 eV and δ = 100 eV, while for SENSEI, we show the exclusion for δ = 10 eV.We verify that our analysis reproduces the elastic scattering (δ = 0) bounds from these experiments, shown as shaded regions in Fig. 8.We note that these bounds are conservative since we do not model any backgrounds or place any cuts in the observed events.In other words, we treat any event as a potential DM signal, resulting in the weakest possible limits on the cross section.
As shown in Fig. 8, the allowed parameter space of this model (represented by the grey band) is an attractive target for upcoming experiments probing light DM.In particular, we plot the sensitivity curves for Oscura [112] for δ = 0 eV (light purple) and δ = 10 eV (purple) assuming an exposure of 30 kg-year.Furthermore, for the O(eV) values of δ we consider here, up-scattering inside the Earth may result in an enhanced density of the excited state at the detector, which would also result in stronger bounds and forecasts than the ones presented here [108].Finally, even stronger bounds and forecasts may be obtained by considering electron ionization caused by DM-nucleon scattering through the Migdal effect [113][114][115].These bounds were evaluated for δ ≳ O(keV) in Ref. [116].We leave an analysis of these bounds for sub-keV values of δ for future work.

V. DISCUSSION
Pseudo-Dirac DM is a minimal modification of standard vector portal DM that can result in qualitatively new cosmological, astrophysical, and experimental phenomenology.In this work, we examine the parameter space of this well-studied model in the regime of small (sub-keV) mass splittings and in the presence of resonant annihilations, χ 2 χ 1 → A ′ → SM SM, where m A ′ ≈ 2m χ .The resonantly enhanced annihilations imply that tiny couplings are able to reproduce the observed DM relic density.Additionally, because of these small couplings, DM can kinetically decouple from the SM before its final relic abundance is reached.Therefore, in order to accurately predict the relic density, one must solve a coupled system of Boltzmann equations for the densities and temperatures of the relevant species.We used the numerical Boltzmann solver DRAKE to properly account for this ef-fect and found that the predicted DM abundance can have corrections as large as an order of magnitude, depending on the underlying parameters of the theory.
The early kinetic decoupling ensures that the excited state is not thermally depleted.Despite the presence of the excited state, this model is consistent with strong bounds coming from BBN and the CMB owing to the strong velocity suppression in the annihilation crosssection at sub-keV temperatures (when the dark photon can no longer be produced on-shell).As a result, as shown in Fig. 2, most of the parameter space of this model is unconstrained, in contrast to sub-GeV Dirac DM that freezes-out through an s-wave process.
The presence of the long-lived excited state can have unique astrophysical signatures that are usually not relevant for pseudo-Dirac DM.For instance, tree-level elastic scattering could cause SIDM behavior, with the caveat that cross sections exceeding σ/m χ ∼ 1 cm 2 /g either have couplings that are excluded by the CMB or lie in a region of parameter space with f * ∼ 0.01.More notably, exothermic downscattering can be relevant, especially in low-velocity environments where there is an enhancement to the cross section.The extent to which exothermic scattering matters in situ is difficult to quantify without further simulation of inelastic SIDM halos, however previous work has found that small mass splittings with δ/m χ ∼ 10 −6 can have a dramatic impact on the properties of a DM halo and its subhalos [62,63].The relic excited state can also result in signals in indirect detection experiments.The ground and excited state present in the MW can annihilate into various SM states and therefore gamma-ray and X-ray telescopes can be used to look for signatures of this model.Since DM annihilation at late times mimics the off-resonance Dirac case, we use previous analyses of gamma-ray and X-ray data [61,68,69] to place bounds on the total DM annihilation cross-section.We find that despite the small couplings, the resonant inelastic parameter space is an attractive target for the next generation of telescopes such as GECCO, MAST, GRAMS and AMEGO (see Fig. 5).
The excited state can downscatter in a range of direct detection experiments.Because of the absence of kinematic barrier for this process, the event rate is significantly enhanced compared to elastic scattering and also compared to pseudo-Dirac DM thermal histories with an exponentially suppressed abundance of the excited state.Using state-of-the-art numerical codes DarkART and QEDark to obtain the electron ionisation and crystal form factors respectively, we calculate the event rates for inelastic DM-electron scattering at Xenon-and Siliconbased experiments.We use the analysis procedure described in Ref. [109] to place bounds on the DM-electron scattering cross-section for different values of δ.We find that future semiconductor-based experiments such as Oscura will begin to probe the resonant inelastic DM pa-rameter space as shown in Fig. 8.
Simultaneously, accelerator searches for dark photons can also explore relevant parameter space for this model.We use the publicly available code Darkcast, to depict bounds on the kinetic mixing, κ for fixed values of g χ .The presence of a resonance implies a broadening of the thermal target as shown in Fig. 6.We find that accelerator bounds are complementary to those set by the CMB, as the accelerator experiments probe the parts of the parameter space that are harder to constrain using earlyuniverse probes.Future experiments will constrain large parts of the parameter space of this model, as shown in Fig. 7.
In summary, sub-GeV resonant pseudo-Dirac DM is an attractive thermal target for a variety of terrestrial DM experiments and astrophysical searches.The complete exploration of the phenomenology of this model leaves a lot of promising directions for future work.In particular, a more accurate treatment of the photo-disassociation bounds coming from BBN may further constrain this parameter space.Additionally, the early kinetic decoupling present in this model may loosen the m χ ≥ 10 MeV lower bound on the mass of thermal DM coming from N eff , despite the thermal equilibrium between the dark and visible sectors at early times.Furthermore, it will be necessary to perform additional cosmological simulations in order to understand the effect of up-and downscattering on halo properties which have immediate consequences for direct and indirect detection experiments.Finally, a more rigorous analysis of the direct detection bounds needs to be undertaken, including (1) the upscattering of the ground state as it passes through the Earth before downscattering in the detector and (2) scattering on electrons through the Migdal effect.Such an analysis may result in even stronger bounds and forecasts on resonant pseudo-Dirac DM than the ones presented here.

1 R
FIG.2.The couplings that yield a DM abundance that matches the observed relic density for various DM masses and values of ϵR.Dashed lines indicate CMB annihilation constraints on those couplings, while solid lines are consistent with the CMB.In this parameter space, the smaller of the two couplings (corresponding to whichever decay channel is the bottleneck) determines the relic abundance, with the larger coupling being irrelevant.Small deviations from this behaviour occur at large values of mχ and ϵR where early kinetic decoupling has a significant effect on the relic abundance.These trends contrast with the parameter space for standard thermal freeze-out where the product of couplings is the most relevant for setting the relic abundance.For comparison, we show the couplings for thermal freeze-out for m A ′ = 3mχ (off resonance) as dot-dashed lines.
FIG.3.The thermally averaged annihilation cross-section for kinetically decoupled (solid) and kinetically coupled (dashed) DM as a function of x (defined with respect to the SM temperature, T ) for ϵR = 0.1 (dark green) and ϵR = 0.001 (light green) and mχ = 100 MeV.We fix gχ = 0.01 and choose κ such that we reproduce the observed DM abundance after solving Eqs.(12)-(13).Early kinetic decoupling can suppress or enhance resonant annihilation at given temperature.

1 FIG. 4 .
FIG.4.The relative number density of the excited and ground states, f * = nχ 2 /nχ 1 , as a function of the DM mass.The different lines correspond to gχ = 0.1 (solid) and gχ = 0.01 (dotted).The two colours correspond to two different values of ϵR and the shaded areas correspond to varying δ between 1−100 eV from top to bottom.Due to the early decoupling of the processes that would deplete the excited state, the excited state remains abundant at late times in most of the parameter space, in contrast to most other pseudo-Dirac DM thermal histories.

7
photon production channels.For γ inv = 10 −5 , the invisible decay width is larger than the visible one for κ ≲ 0.1 and therefore the beam dump experiments that search for A ′ decays to leptons lose sensitivity.The solid black lines in both panels correspond to the κ values that reproduce the observed DM abundance following the production mechanism outlined in Sec.II, for ϵ R ∈ [0.001, 0.1].The dashed parts of the black lines correspond to the points excluded by the CMB.We note that for Γ DM > Γ SM , (or equivalently g χ ≫ κ), the CMB constrains the part of the parameter space which cannot be probed by accelerator experiments making the accelerator and cosmological bounds highly complementary.Additionally, we find that for a given γ inv and ϵ R , or equivalently, a fixed g χ , there exists a maximum m A ′ beyond which DM is always over-produced.From Eq. (14), we see that the DM relic density is proportional to m A ′ / min(g 2 χ , κ 2 ) meaning that once the couplings are fixed, a larger m A ′ results in a larger DM abundance.This causes the relic lines to be vertical in the right panel of Fig.6(in the left panel, the maximum m A ′ lies outside the plotted region).Finally, we note that for a given γ inv , κ-values bounded by the two curves corresponding to ϵ R = 0.001 and ϵ R = 0.1 also reproduce the observed DM abundance for different values of ϵ R resulting in a much broader thermal target.For reference, we show the usual thermal target assuming non-resonant, s-wave thermal freeze-out (dotted black line) in the left panel of Fig.6.Future accelerator experiments are poised to more fully explore the parameter space of this model, as shown in Fig.7

FIG. 8 .
FIG. 8.Constraints on the DM-electron scattering crosssection, σe as a function of the DM mass from XENON10 (blue), XENON1T (green), and SENSEI (orange) for different values of δ, with f * = 1.The shaded regions correspond to exclusions on elastic scattering, i.e., δ = 0.The dashed purple lines are sensitivity curves for Oscura assuming δ = 0 (light) and δ = 10 eV (dark) respectively.In grey, we show the resonant iDM parameter space allowed by CMB.