{'fraction_non_alphanumeric': 0.03891526479058028, 'fraction_numerical': 0.018946174629471086, 'mean_word_length': 4.278478783966195, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Inelastic Dirac Dark Matter

Feebly interacting thermal relics are promising dark matter candidates. Among them, scenarios of inelastic Dark Matter evade direct detection by suppressed elastic scattering off atomic nuclei. We introduce inelastic Dirac Dark Matter, a new model with two Dirac fermions in the MeV-GeV mass range. At feeble couplings, dark matter can depart from chemical as well as kinetic equilibrium with the Standard Model before freeze-out. In this case, the freeze-out is driven by conversion processes like coscattering, rather than coannihilation. We show that inelastic Dirac relics are consistent with cosmological observations, in particular with nucleosynthesis and the cosmic microwave background. Searches for dark sectors at colliders and fixed-target experiments, in turn, are very sensitive probes. Compared to the strongly constrained pseudo-Dirac scenario, inelastic Dirac Dark Matter offers a new search target for existing and upcoming experiments like Belle II, ICARUS, LDMX and SeaQuest.


Introduction
Direct detection has put pressure on thermal WIMPs. Current searches for dark matternucleon scattering are extremely sensitive, and null results have ruled out many scenarios of weakly interacting massive particles (WIMPs) as candidates for cold thermal dark matter (DM) [1,2]. 1 An elegant option to evade direct detection constraints is to consider inelastic Dark Matter, where elastic nucleon scattering is absent or parametrically suppressed and inelastic up-scattering into a heavier dark partner is kinematically suppressed [6]. Inelastic dark matter has become a benchmark target for searches in particle physics and astrophysics [7]. In a minimal realization, commonly dubbed iDM, the dark sector consists of a pseudo-Dirac state with two Majorana fields, ξ 1 and ξ 2 , interacting with quarks and leptons via an abelian dark force by exchanging a dark photon A [8]. This interaction drives the dark matter freeze-out in the early universe and the relic abundance is set through coannihilation ξ 1 ξ 2 → A * → ff into Standard Model (SM) fermions f . The abundance measured today [9] favors iDM candidates in the MeV-GeV mass range. The phenomenology of this predictive scenario has been investigated in detail [8,[10][11][12][13][14][15][16][17][18][19][20][21][22], specifically for direct detection [11,19,21,22], at colliders [16][17][18], and at fixed-target experiments [12][13][14][15]20]. Taken together, searches for light dark particles in all three areas have excluded most of the parameter space of iDM. Inelastic dark matter from thermal freeze-out via coannihilation appears to be strongly constrained.
In this work, we introduce inelastic Dirac Dark Matter (i2DM) as a new model for feebly coupling dark matter. We promote the two dark fermions in iDM to Dirac fields, one being charged and the other one uncharged under a dark U (1) gauge symmetry. The symmetry is spontaneously broken by a Higgs-like mechanism, which causes the two dark fermions to mix. As a result, the dark matter candidate χ 1 interacts only feebly through the small mixing, while the coupling of the dark partner χ 2 through the dark force is unsuppressed.
This moderate modification of iDM leads to a very different cosmology: In i2DM, the relic abundance can be set by partner annihilation χ 2 χ 2 → ff or coscattering χ 1 f → χ 2 f , which is not an option in iDM where thermal freeze-out is necessarily driven by coannihilation. In contrast to the well-studied freeze-out through coannihilation and partner annihilation [23,24], the role of conversion processes like coscattering has been investigated in specific dark matter scenarios only recently [25][26][27][28][29][30][31][32]. Our goal is to demonstrate that i2DM is a new cosmologically viable candidate for feebly coupling inelastic dark matter in the MeV-GeV range that can be probed at current and future experiments.
In Sec. 2, we introduce inelastic Dirac Dark Matter and explain the main characteristics of the model. The formalism is described in more detail in App. A. In Sec. 3, we analyze the density evolution of the dark sector in the early universe before freeze-out. We pay special attention to deviations from chemical and kinetic equilibrium, which occur for small couplings. In this regime, computing the relic abundance requires to solve a coupled system of Boltzmann equations and to keep track of the dark matter momentum distribution. Details on our calculations can be found in App. B. In Sec. 4, we investigate possible effects of i2DM on astrophysical and cosmological observables. In particular, we discuss the impact of the QCD phase transition, effects on the formation of light elements, imprints on the cosmic microwave background, as well as constraints from supernova cooling. The resulting bounds from astrophysics and cosmology set a clear search target for i2DM. In Sec. 5, we test this new dark matter target at laboratory experiments. Signatures of i2DM strongly depend on the lifetime of the dark partner. We find that fixed-target and flavor experiments are most sensitive to dark fermions with long lifetimes through searches for displaced decays, scattering, and missing energy. In Sec. 6, we conclude with an outlook to future experiments that can conclusively test inelastic Dirac Dark Matter.

Inelastic Dirac Dark Matter
We introduce a dark sector consisting of two Dirac fermions χ 1 and χ 2 with masses m 1 and m 2 , interacting with the SM particles via a dark photon A with mass m A . The lighter of the dark fermions, χ 1 , serves as a dark matter candidate. The dark partner, χ 2 , will play a crucial role for the interactions between the dark and visible sectors.
Such a scenario can be constructed from a renormalizable theory with two fermion fields, both of them SM gauge singlets and one of them charged under a new abelian gauge symmetry U (1) D . This symmetry is spontaneously broken by the vacuum expectation value of a dark scalar, which induces mixing between the dark fermions and gives the dark gauge boson a mass. The dark gauge boson kinetically mixes with the hypercharge field, thus acting as a mediator between the dark fermions and the SM fermions. For details on the model we refer the reader to App. A. A FeynRules model for i2DM is available at [33].
A priori, such a dark sector could be realized at any mass scale. Throughout this work we focus on dark particles in the MeV-GeV range, which could be resonantly produced at colliders and fixed-target experiments. The relevant interactions of the dark sector with visible matter are described by the Lagrangian s W c W Z µ sin 2 θJ µ 1 − sin θ cos θJ µ 12 + cos 2 θJ µ 2 , (2.1) the mixing is much heavier than the other dark particles and does not affect the observables we consider. In general, the presence of a light dark scalar could lead to interesting effects [17,34] and deserves a dedicated analysis. The phenomenology of i2DM is described by six independent parameters m 1 , ∆, m A , α D , , θ (2.3) with α D = g 2 D /(4π). The relative mass splitting between the dark fermions is defined as Throughout this work we focus on the mass hierarchy so that decays of the dark photon into dark fermions are kinematically allowed. 2 The total decay width of the dark photon is given by where α e = e 2 /(4π) is the fine structure constant, while Γ(A → χχ) and Γ(A → SM) denote the (normalized) decay rates into dark fermions and into leptons and hadrons, respectively. Dark photon decays into pairs of dark fermions, χ i χ j , and leptons, i¯ j , are described by the kinematic function Dark photon decays into hadrons can be computed by rescaling the leptonic decay rate with e + e − data [36]. For 2 α e α D , decays into SM particles are suppressed and the dark photon mostly decays into dark fermions. The corresponding decay rate is Γ(A → χχ) = sin 4 θ Γ(m 1 , m 1 ) + sin 2 (2θ) Γ(m 1 , m 2 ) + cos 4 θ Γ(m 2 , m 2 ) . (2.8) For m 1 ≈ m 2 , the branching ratios are determined to a good approximation by the dark fermion mixing, so that B(A → χ 1 χ 1 ) ≈ sin 4 θ , B(A → χ 1 χ 2 , χ 2 χ 1 ) ≈ sin 2 (2θ) , B(A → χ 2 χ 2 ) ≈ cos 4 θ . (2.9) The freeze-out dynamics relies on dark fermion annihilation into SM particles. Annihilations into leptons via χ i χ j → A * → + − can be calculated in perturbation theory. Annihilations into hadrons can be predicted by re-scaling the cross section for annihilation into muons with the measured ratio [36,37] R(s) = σ(e + e − → hadrons) σ(e + e − → µ + µ − ) . (2.10) The total cross section for dark fermion annihilation at a center-of-mass energy √ s is then given by (2.11) In dark fermion interactions with SM fermions the dark photon acts as a virtual mediator. As a consequence, for m A m 1,2 the scattering rates and decays of dark fermions scale as [10] (2.12) As long as the dark photon is heavy compared to the momentum scale probed in observables, the dark sector interactions are described in terms of the four parameters The phenomenology of i2DM crucially relies on the properties of the dark partner. For m 2 1 GeV, the dark partner decays to almost 100% into leptons [17]. The decay rate via a heavy virtual dark photon is given by where we have neglected the lepton mass in the final state. We neglect hadronic decays in our analysis. We will refer to this model as inelastic Dirac Dark Matter (i2DM) to distinguish it from the widely studied scenario with Majorana fermions, often called pseudo-Dirac Dark Matter or simply inelastic Dark Matter (iDM) [8, 15-18, 20, 21]. iDM builds on a single pseudo-Dirac fermion charged under a new U (1) D force. The dark gauge symmetry is spontaneously broken by Majorana mass terms for the chiral components of the Dirac spinor. The resulting dark sector contains two Majorana fermions, ξ 1 and ξ 2 , that couple to the dark photon mostly through inelastic interactions In both iDM and i2DM models, elastic scattering off atomic nuclei is suppressed, which strongly reduces the sensitivity of direct detection experiments. In i2DM, elastic dark matter interactions are additionally suppressed for small dark fermion mixing θ. We will come back to direct detection in Sec. 5.1.
Despite similar predictions for nucleon scattering, i2DM and iDM feature very distinct dark matter dynamics in the early universe. The difference lies mostly in efficient A χ 2χ2 interactions in i2DM, which are suppressed or completely absent in iDM. At first sight, this appears as a small modification. However, as we will show, the presence of partner interactions in i2DM has drastic effects on the dark matter freeze-out. We obtain new dark matter candidates with couplings that would be too small to explain the relic abundance from coannihilation, as in iDM.

Freeze-out at feeble couplings
Dark matter relics in the MeV-GeV range must be feebly coupled to the thermal bath in order to account for the observed DM abundance, Ω χ h 2 = 0.12 [9]. Moreover, viable scenarios of inelastic dark matter typically require a compressed spectrum of dark-sector particles. For i2DM, this leads to the parameter region of interest Within this regime, the relic abundance can be set by various mechanisms. We consider scenarios where dark matter is in kinetic equilibrium with the thermal bath over a period of time before the relic abundance is set. The relic abundance should therefore be set by a freeze-out process, and we expect the dark fermions to be non-relativistic around the freeze-out temperature 3 In general, the freeze-out dynamics is determined by the evolution of the density distribution functions f i (x, q i ) for all relevant species i, expressed in terms of the dimensionless time and momentum variables Here T is the temperature of the thermal bath, and we have chosen the mass of the lightest dark fermion, m 1 , to normalize x. The norm of the three-momentum of species i is denoted as p i = | p i | and scales with the scale factor a as 1/a. The comoving number density, Y i , is obtained from the phase-space integration of the distribution function as where g i is the number of degrees of freedom of species i. Both the number density, n i (x), and the entropy density of the universe, s(x), scale as 1/a 3 . In what follows, we denote the 3 Other possible mechanisms for feebly coupled dark matter do exist. For instance, in the case of freezein [38,39], dark matter was never in kinetic equilibrium and is usually relativistic at the time where its abundance freezes.
number density of χ i in kinetic equilibrium and with zero chemical potential as n eq i . For later convenience, we define the ratios of number densities, where n eq = n eq 1 + n eq 2 is the total dark sector equilibrium number density. The freeze-out temperature T fo is determined by the time x fo = m 1 /T fo at which the comoving dark matter density approaches the dark matter abundance Y 0 observed today, In our numerical analysis, we determine x fo by requiring that the dark matter density at freeze-out satisfies As we will see, i2DM dark matter candidates are not always in chemical and/or kinetic equilibrium with the thermal bath around the freeze-out temperature. As a result, freezeout as defined in Eq. (3.6) does not necessarily coincide with chemical decoupling, as in the case of vanilla WIMP dark matter [40]. We therefore define the times x 1 and x 2 , where χ 1 and χ 2 chemically decouple from the bath, corresponding to the decoupling temperatures T 1 and T 2 . The time evolution of the density distribution function of a particle species i is described by the Boltzmann equation (3.8) Here E i = E i (x, q i ) is the energy associated with a species of mass m i and momentum q i , and H = H(x) is the Hubble rate at time x. In a radiation-dominated era, the Hubble rate scales as H(x) ∼ x −2 . The collision term C describes interactions of species i with all other involved species j. We emphasize that Eq. (3.8) holds even if species i is not in kinetic and/or chemical equilibrium with the thermal bath, but is only valid as long as the number of relativistic degrees of freedom in the universe is constant. 4 In what follows, we will refer to Eq. (3.8) as the unintegrated Boltzmann equation. By making several simplifications, the set of N partial integro-differential equations for species i, j = {1, . . . N } can be reduced to a set of ordinary differential equations. See App. B for details. If all species are in kinetic equilibrium the relic abundance can be calculated from the time evolution of the number densities n i (x) from Eq. (3.4) [23,41]. The corresponding evolution equations are referred to as integrated Boltzmann equations.
In the remainder of this section, we will discuss the freeze-out dynamics for i2DM. In Sec. 3.1, we analyze all relevant annihilation, scattering and decay processes that can play a role in setting the DM relic abundance. Depending on the relative importance of these processes, we encounter different phases of freeze-out. In Sec. 3.2, we discuss these phases in detail and explain how to account for deviations from chemical and kinetic equilibrium when computing the relic abundance.

Dark matter interactions
For i2DM, the relevant interactions of the dark fermions with the thermal bath entering the collision term C in Eq. (3.8) are (co)annihilation: (co)scattering: At temperatures below the GeV scale, f denotes all (hadronized) quarks and leptons that are in equilibrium with the thermal bath. To determine the relevance of the various processes for DM freeze-out, we investigate the thermally averaged interaction rates Γ of the dark fermions with the bath. We define the interaction rates for annihilation, scattering and decays of particle i as (3.10) Here and below, i, j = {1, 2} denote the dark fermions χ 1 and χ 2 , and k, l label the SM fermions f . The brackets indicate the thermal average. The reaction densities γ can be expressed as 5 γ ab→cd = n eq a n eq b σ ab→cd v ab , γ a→bcd = n eq a Γ a→bcd where a, b, c, d can denote particles from both the visible and the dark sector, σ ab→cd v ab is the thermally averaged cross section with particles a and b in the initial state, Γ a→bcd is the decay rate of particle a in its rest frame, and K 1 , K 2 denote the modified Bessel functions of the first and second kind. The Moeller velocity is defined by where p a , p b denote the 4-momenta of particles a and b, and E a , E b are their energies in the center-of-mass frame. 6 For simplicity, we label the reaction densities only by the involved dark fermions as 5 Reaction densities have been widely used in the context of leptogenesis, as in Ref. [42]. They allow us to define the relevant interaction rates in a thermal bath for the production of feebly coupled dark matter in a compact and unambiguous way [30,43]. For more details see App. B. 6 In the text, by default pa refers to the norm of the 3-momentum. Only in the Moeller velocity or in a delta function that enforces 3-momentum and energy conservation, pa denotes the 4-momentum.
The relevance of the various interaction rates at a certain temperature can be inferred from their scaling with the time variable x = m 1 /T and with the model parameters. For nonrelativistic dark sector particles χ 1 and χ 2 and relativistic involved SM fermions around freeze-out, the reaction densities for (co)annihilations in thermal equilibrium scale as γ 11 = n eq 1 n eq 1 σ(χ 1 χ 1 → ff ) v 11 ∝ x −3 e −2x y sin 4 θ (3.14) The (co)scattering reaction densities scale as Three-body decays and inverse decays yield the reaction density given that For processes with leptons in the final state we have calculated all rates analytically. To compute annihilations into hadrons, we rescale the cross section as described in Sec. 2, following Ref. [17]. For scattering processes we only include the dominant scatterings χ ± → χ ± off leptons = {e, µ}, which are still relativistic at sub-GeV temperatures. We neglect scatterings off hadrons, whose number densities are Boltzmann-suppressed for masses above the muon mass.
In the expressions above, the variable y determines the effective strength of the dark force and thereby the overall efficiency of dark sector interactions. As we will discuss in Secs. 4 and 5, y is constrained by cosmology and laboratory searches. In addition, the dark fermion mixing θ must be small to circumvent bounds from direct detection experiments. Therefore dark matter annihilation γ 11 and scattering γ 1→1 must be strongly suppressed. In the absence of further interactions, such a suppression would lead to an overabundance of dark matter today. As a consequence, feebly interacting i2DM candidates cannot be thermal WIMPs in the classical sense, where the relic abundance is determined by the WIMP pair-annihilation rate at freeze-out.
Indeed, for feebly coupling i2DM, dark matter pair annihilation and scattering with the thermal bath play no role for the temperature evolution of the dark matter density. Instead, the evolution of χ 1 is driven by interactions with the dark partner χ 2 . The dark partner is kept in chemical and kinetic equilibrium with the bath via efficient annihilation and scattering, driven by the reaction densities γ 22 and γ 2→2 . However, at freeze-out all interaction rates of χ 2 are exponentially Boltzmann-suppressed by powers of e −x∆ , compared to the interactions of χ 1 . For χ 2 interactions to impact the freeze-out of χ 1 , the mass difference ∆ must be small, as in Eq. (3.1). Viable scenarios of i2DM with small dark matter couplings require a compressed spectrum of dark fermions.
Conversions χ 1 ↔ χ 2 play an essential role in i2DM freeze-out; they keep dark matter in equilibrium with the thermal bath and impact the evolution of the number density. Both coscatterings and decays contribute to the conversion rate of χ i , Due to the respective scaling of the reaction densities with x, see Eqs. (3.15) and (3.16), the suppression of the thermal rates at low temperatures is stronger for coscattering (γ scat 2→1 ∝ x −9/2 ) than for decays (γ dec 2→1 ∝ x −3/2 ). Depending on their relative amplitude at a given temperature, either process can dominate the conversion rate Γ i→j and thus the evolution of the number densities.
To illustrate the impact of the various processes on i2DM, in Fig. 1 we show the time evolution of the dark sector interaction rates Γ and the yield Y for three i2DM benchmarks. The three benchmarks belong to different phases of freeze-out, which we will discuss in detail in Sec. 3.2.

Phases of freeze-out
From the discussion in Sec. 3.1, it becomes clear that the freeze-out dynamics should be very sensitive to the parameters y, θ and ∆. When successively decreasing the dark interaction strength y, we identify three different phases of freeze-out, distinguished by the processes that set the dark matter relic abundance: 1. coannihilation phase: Ω χ h 2 set by χ 1 χ 2 ↔ ff and χ 2 χ 2 ↔ ff , 2. partner annihilation phase: Ω χ h 2 set by χ 2 χ 2 ↔ ff , 3. conversion phase: Ω χ h 2 set by χ 1 f ↔ χ 2 f and/or χ 2 ↔ χ 1 ff .
In addition, the thermal history of the dark matter candidate depends on whether departures from chemical or kinetic equilibrium with the bath have occurred prior to freeze-out. This critically depends on the efficiency of χ 1 ↔ χ 2 conversions at the time x 2 at which χ 2 chemically decouples. We distinguish between three regions of parameter space, where the following conditions are satisfied: where Γ 1→2 denotes the conversion rates from Eq. (3.17) and H(x 2 ) is the Hubble expansion at the time x 2 . These three regions allow us to systematically study the effects of dark matter chemical and kinetic decoupling before freeze-out on the relic abundance. 10 -3 10 -3 Chemical decoupling of χ 2 occurs when its (co)annihilation rate drops below the Hubble rate. If conversions χ 1 ↔ χ 2 are efficient around x 2 , the decoupling time can roughly be estimated using and χ 2 decouples at the same time as χ 1 . However, if conversions are absent, χ 2 chemically decouples around Numerically we determine the time x i where χ i chemically decouples by requiring that the density yield deviates from equilibrium by 20%, The classification made above allows us to understand the density evolution of the dark fermions shown in Fig. 1. The three i2DM benchmarks correspond to the freeze-out phases of partner annihilation (upper left plot) and conversion (upper right and lower plots). Partner annihilation is relevant in region (A), while conversion can prevail either in region (B) in kinetic equilibrium (upper right plot), or in region (C) beyond kinetic equilibrium (lower plot). The interplay between the different freeze-out phases and decoupling regions is shown in Fig. 2 as a function of the model parameters y and tan θ for two scenarios with fixed dark matter masses. In the upper plot of Fig. 2, the three benchmark i2DM scenarios from Fig. 1 are marked as green bullets. Below we first discuss Fig. 1 in detail and then turn to Fig. 2.
In Fig. 1, all benchmarks correspond to fixed parameters m 1 = 60 MeV, ∆ = 0.05 and α D = 1/4π. In each of the plots, the top panel shows the evolution of the various interaction rates Γ , normalized to the Hubble rate. The conversion rate Γ 1→2 that distinguishes regions (A), (B), and (C) is driven by y tan 2 θ, see Eq. (3.15). It decreases when going from partner annihilation in region (A) (upper left plot) to conversion beyond kinetic equilibrium in region (C) (lower plot). The relative scaling of the interaction rates is determined by the reaction densities from Eqs. (3.14), (3.15) and (3.16) discussed in Sec. 2. In particular, the reaction densities for (co)annihilations, γ ij ∼ e −2x , drop faster at low temperatures than scattering and decays, γ i→j ∼ e −x . The relative strength of (co)annihilations depends exponentially on the mass splitting ∆ and also on the dark fermion mixing θ. Due to the small splitting and mixing in the three benchmarks, partner annihilation dominates (light green), followed by coannihilation (dark green) and suppressed dark matter annihilations (blue). As mentioned in Sec. 2, conversions through scattering (purple) decrease faster with time than inverse decays (orange). Around freeze-out, however, conversions dominate over decays in all three benchmarks.
In the bottom panels of Fig. 1, we show the time evolution of the dark fermion comoving number densities, Y i (x), (solid) and the equilibrium yield, Y eq i (x), for comparison (dashed). In the lower plot, we also indicate the evolution that is obtained when neglecting the kinetic decoupling of dark matter (dotted). To guide the eye, we highlight the times for freeze-out and chemical decoupling with black vertical lines, determined by Eqs. (3.7) and (3.23).
We now turn our attention to Fig. 2, which illustrates the different phases of freezeout for i2DM as a function of the dark interaction strength y and the dark fermion mixing tan θ for two fixed dark matter masses m 1 = 60 MeV (top) and 150 MeV (bottom). The regions (A), (B) and (C), corresponding to decreasingly efficient conversions as suggested by Eqs. (3.18)- (3.20), are delineated with dashed gray lines. The exact relations between the conversion rate and the Hubble rate along these lines are given in Eqs. (3.32) and (3.33). The observed relic abundance Ω χ h 2 = 0.12 is obtained along the solid colored contours in the (y, tan θ) plane for fixed values of the mass splitting ∆. When increasing the dark matter mass m 1 , the contours shift to the right, meaning that the observed abundance is obtained for larger values of y. This is easily understood, as all (co)annihilation and conversion rates scale as Γ ij , Γ i→j ∝ y/m 2 1 . In the upper part of the plots, all contours converge and the mass splitting ∆ plays no role in setting the relic abundance. Here the abundance is set by pair annihilations χ 1 χ 1 → ff , see Eq. (3.14). In Sec. 5, we will see that laboratory searches exclude this region of parameter space. As a result, we focus on dark matter candidates with small couplings, corresponding to the phases of coannihilation and partner annihilation in region (A), and on conversions in regions (B) and (C).
In region (A), efficient conversion rates satisfying Eq. (3.18) keep χ 1 in chemical and kinetic equilibrium with the thermal bath until freeze-out. In particular, efficient conversions ensure that χ 1 and χ 2 have equal chemical potentials, so that their number densities are related by In this case, the freeze-out conditions are similar to the ones of a thermal WIMP and the observed relic dark matter abundance is obtained for a freeze-out time [23,24,41] x WIMP 25 . , conversion processes are less efficient and the rates Γ 1→2 (x 2 ) are about ten to one hundred times larger than the Hubble rate. As a result, the dark matter density departs from chemical equilibrium prior to freeze-out. The effect is visible in the second benchmark displayed in the upper right plot of Fig. 1, corresponding to the second green bullet in the upper plot of Fig. 2.
Once the conversion rate Γ 1→2 (x 2 ) is further suppressed, dark matter cannot be expected to be kept in kinetic equilibrium until freeze-out. We enter region (C), defined by Eq. (3.20) and illustrated by a benchmark in the lower plot of Fig. 1, corresponding to 10 -10 10 -9 10 -8 10 -7 10 -6 10 -5 10 -5 the lower green bullet of the upper plot of Fig. 2. To quantify the effect of kinetic decoupling, in Fig. 1 we display the physical dark matter yield obtained by respecting kinetic decoupling (solid purple curve) compared to the same yield obtained when neglecting kinetic decoupling (dotted purple curve).
In what follows, we will discuss the physics of the three freeze-out phases in detail, paying particular attention to non-equilibrium effects in the conversion phase. Further details about the computation of the dark matter relic abundance in the presence of chemical and kinetic decoupling can also be found in App. B.

Coannihilation
Coannihilation sets the relic abundance when the rate of χ 1 χ 2 ↔ ff is larger than the Hubble rate and dominates over dark matter pair annihilations χ 1 χ 1 ↔ ff around the freeze-out time, i.e., when According to Eqs. (3.10) and (3.14), in i2DM the ratio of these rates scales as One might deduce that coannihilation sets the relic abundance if e x fo ∆ tan 2 θ 1. However, this is only a necessary condition, because in this regime the relic abundance could also be driven by partner annihilation χ 2 χ 2 → ff . To determine the relative impact of partner annihilation, one needs to consider the weighted ratio of coannihilation and partner annihilation rates, r 1 Γ 12 (x)/r 2 Γ 22 (x). 8 In i2DM, this ratio scales as e x∆ tan 2 θ, exactly as in Eq. (3.27). As a result, in the coannihilation phase the comoving number density Y 1 (x) freezes once the weighted sum of coannihilation and partner annihilation drops below the Hubble rate, 9 The relic dark matter abundance is set around the freeze-out temperature of a thermal WIMP. Freeze-out through coannihilation is realized in the upper part of region (A) in Fig. 2. Due to the scaling Γ 12 ∼ y sin 2 (2θ), the relic abundance contours tend to larger interaction strength y at smaller mixing θ in this regime.
In the coannihilation phase, i2DM resembles iDM, where the ξ 1 − ξ 2 couplings of Eq. (2.15) prevail and coannihilations ξ 1 ξ 2 → ff are the dominant number-changing interactions of the dark matter candidate ξ 1 with the bath. However, lowering the dark interaction strength y leads to inefficient coannihilation around freeze-out. In iDM, this results in an overabundance of dark matter. In i2DM, suppressed coannihilations can be compensated by efficient partner annihilations and conversions, thus explaining the observed dark matter abundance even if the dark sector is feebly coupled. 8 The effective coannihilation and partner annihilation rates entering in the computation of the dark matter relic abundance have to be weighted by ri, the ratio of the dark species equilibrium densities to the total one, when the dark sector species are in chemical equilibrium, see e.g. [23]. 9 This relation assumes that χ1 and χ2 are both in kinetic and chemical equilibrium prior to freeze-out.

Partner annihilation
If the coannihilation rate is suppressed compared to the Hubble rate around x = x WIMP , the relic DM abundance can be set by partner annihilations χ 2 χ 2 → ff . In this phase, the freeze-out condition reads 10 As in the coannihilation phase, the freeze-out time is fixed to x fo x WIMP , up to a moderate logarithmic dependence on the model parameters [23]. In Fig. 2, the phase of partner annihilation lies in region (A) and is characterized by vertical lines. Due to the scaling r 2 Γ 22 ∼ y cos 4 θ e −2x∆ , the relic dark matter abundance is essentially independent of the mixing as cos 4 θ ≈ 1. To obtain the observed abundance, variations of the dark interaction strength y can be compensated by the mass splitting ∆, as illustrated by the various contours. In the partner annihilation phase, χ 1 ↔ χ 2 conversions have to be efficient enough to satisfy the decoupling condition from Eq. (3.18). As a consequence, χ 1 and χ 2 chemically decouple from the bath around the same time. This is visible in the upper left plot of Fig. 1. Notice that χ 2 → χ 1 decays happen on a time scale shorter than the period of chemical decoupling and do not affect the freeze-out time x fo ≈ x WIMP .
Whenever coannihilation or partner annihilation set the relic abundance, dark matter is in chemical and kinetic equilibrium until decoupling. In this case, the evolution equations from Eq. (3.8) can be reduced to one single integrated Boltzmann equation, written in terms of the total number density of all dark particles, n = i n i , as commonly used for thermal WIMPs [23,24]. The dark matter abundance can be computed with any of the available Boltzmann solvers [44][45][46], which explicitly make use of Eq. (3.24). For the computations in this work we have used our own Boltzmann solver. We have verified that our results for the relic abundance in region (A) in Fig. 2 and in the left panel of Fig. 1 agree with the results obtained from micrOMEGAs [44].

Conversion
As we discussed above, partner annihilation can set the relic dark matter abundance even if χ 1 χ 1 → ff and χ 1 χ 2 → ff annihilations are suppressed. If the dark fermion mixing θ is very small, χ 1 ↔ χ 2 conversion rates can be comparable to or even fall below the Hubble rate around χ 2 chemical decoupling, as in Eqs. (3.19) and (3.20). The dark matter abundance is now set by conversion processes, i.e., by coscattering and/or (inverse) decays.
In i2DM, conversion-driven freeze-out can occur in regions (B) and (C) in Fig. 2, Here coscattering dominates over decays in the thermal history around freeze-out. Due to the scaling Γ 1→2 ∝ sin 2 (2θ), the contours of constant Ω χ h 2 are sensitive to the mixing angle θ. For a fixed interaction strength y and mass splitting ∆, the relic abundance in this regime is generally larger than what one would expected from partner annihilation. The reason is that χ 1 ↔ χ 2 conversions are less efficient and the dark matter yield Y 1 (x) can start to deviate from the equilibrium yield Y eq 1 (x) well before the freeze-out time x = x fo [25,28]. The latter effect is illustrated in the upper right and lower plots of Fig. 1. The increased yield has to be compensated by a larger interaction strength y, causing the contours in Fig. 2 to bend towards the lower right corner. Notice that efficient partner annihilation is essential for conversions to explain the observed relic abundance. Suppressed partner annihilations would result in an overabundance of dark matter.
To describe deviations of the dark matter density from chemical and kinetic equilibrium within the conversion phase, it is convenient to distinguish three key moments: 1. the time x 1 at which χ 1 chemically decouples from the bath; 2. the time x 2 at which χ 2 chemically decouples from the bath; 3. the dark matter freeze-out time x fo , which can differ from x WIMP 25 in the conversion phase.
Numerically, we determine x 1 , x 2 and x fo using Eqs. (3.23) and (3.7). Deviations of the dark matter density from chemical equilibrium around freeze-out typically occur for Deviations from chemical and kinetic equilibrium can occur for Thanks to efficient elastic scattering, the dark partner χ 2 is kept in kinetic equilibrium with the bath throughout the dark matter freeze-out process and in particular for x > x 2 after χ 2 chemical decoupling.
Deviations from chemical equilibrium If x 1 x 2 , conversions are barely efficient during χ 2 chemical decoupling, as in Eq. (3.19). This scenario corresponds to region (B) in Fig. 2. Numerically we find that the boundary between regions (A) and (B) corresponds to using Eq. (3.23) to evaluate x 2 . Below this boundary, the χ 1 density departs from chemical equilibrium before χ 2 . This effect is illustrated for a benchmark in the upper right plot of Fig. 1, where the freeze-out process for Y 1 (x) terminates around x fo > x 1 . At that time, we expect that conversions are still sufficiently active to keep χ 1 in kinetic equilibrium with the thermal bath via χ 2 . In particular, we assume that the dark matter distribution function f 1 (x, q) is well approximated by the Boltzmann distribution f eq 1 (x, q) throughout the entire evolution process. 11 To account for deviations from chemical equilibrium in region (B), the Boltzmann equations commonly used for (co)annihilating dark matter [23] have to be supplemented 11 Conversion-driven freeze-out with deviations from chemical equilibrium was studied before in Ref. [25], where f1(x, q) was observed to depart from f eq 1 (x, q) prior to dark matter freeze-out. In i2DM, however, we expect that the distribution function of χ1 resembles f eq 1 (x, q) before freeze-out, because χ1 is kept in kinetic equilibrium at early times via χ2 → χ1 decays and coscatterings (with χ2 being in kinetic equilibrium with the bath), see Fig. 1. This was not the case for the dark matter model studied in Ref. [25]. by explicitly including the conversion rate in the coupled system of evolution equations for Y 1 (x) and Y 2 (x). For details on our implementation of the Boltzmann equations we refer the reader to App. B.1.
Deviations from chemical and kinetic equilibrium For x 1 x 2 , conversions become inefficient before χ 2 chemically decouples, and the condition of Eq. (3.20) is satisfied. This scenario corresponds to region (C) in Fig. 2. Numerically, we find that the boundary between regions (B) and (C) is given by with x 2 evaluated with Eq. (3.23). Below this line, the mixing θ is so small that conversions fail to keep dark matter in kinetic equilibrium. This effect is visualized in the lower panel of Fig. 1: The dark matter yield Y 1 (x fo ) including departures from kinetic equilibrium (solid curve) is larger than under the assumption of kinetic equilibrium (dotted curve). The time between dark matter chemical decoupling and freeze-out is now stretched over a larger range between x 1 and x fo .
In region (C), the phase-space density f 1 (x, q) is expected to deviate significantly from a Maxwell-Boltzmann distribution for x > x 1 . To account for deviations from kinetic equilibrium, the unintegrated Boltzmann equations from Eq. (3.8) have to be solved for the momentum-dependent phase-space density of χ 1 [25,26]. Further details are given in App. B.2. Using our own implementation, we obtain the solid purple curve in the lower plot of Fig. 1. For comparison, we show the evolution of Y 1 (x) obtained using the integrated Boltzmann equations relevant in region (B), but neglecting deviations from kinetic equilibrium (dotted line). The deviations are modest, ranging around 35%. 12 Throughout this work and in particular in Figs. 2, 3 and 4, we use the set of integrated Boltzmann equations from App. B.1, unless specified otherwise. This method is computationally much faster than solving the unintegrated Boltzmann equations and reproduces the exact results to a good approximation.

Bounds from cosmology and astrophysics
The parameter space of i2DM is constrained by several observables in cosmology and astrophysics. In general, the presence of light new particles with masses in the MeV-GeV range affects the thermal history of the universe. Two main effects can be distinguished. First, the presence of additional particles in the thermal bath changes the evolution of the Hubble expansion and of the entropy density of the universe. This may significantly affect the QCD phase transition, Big Bang Nucleosynthesis (BBN), the Cosmic Microwave Background (CMB) or supernova cooling. Second, the annihilation or decay of new particles injects energy into the thermal plasma by inducing excitations, ionisation and heating. Such effects can modify BBN and the CMB compared to the predictions of the cosmological standard model.
In this section, we study the impact of i2DM dark fermions on cosmological and astrophysical observables. In Sec. 4.1, we consider bounds on particles that freeze-out around the QCD phase transition, at temperatures T QCD ∼ 200 MeV. In Sec. 4.2, we derive constraints from BBN around T BBN ∼ 0.1 MeV. In Sec. 4.3, we discuss effects on the CMB around T ∼ 1 eV, while in Sec. 4.4 we report on the constraints from the measurement of the effective number of neutrinos, N eff , at the CMB and BBN times. Finally, in Sec. 4.5, we consider constraints on dark sector particles escaping supernovae.

QCD phase transition
If new particles freeze-out around the GeV scale, the QCD phase transition around T QCD ∼ 200 MeV affects the relic abundance [47]. The confinement of quarks and gluons into hadrons reduces the effective number of relativistic degrees of freedom contributing to the entropy density, h eff [48]. Calculations of h eff around the phase transition are subject to significant uncertainties, leading to variations of about ten percent in the relic abundance [49][50][51]. However, the effect of the QCD phase transition on the thermal evolution of light new particles can be much larger than the mentioned uncertainties. In particular, it can affect the relic abundance of dark matter candidates that freeze out during or shortly after the phase transition.
In Sec. 3.2, we have investigated the freeze-out dynamics of i2DM separately from effects of the QCD phase transition. If the freeze-out occurs in thermal equilibrium, dark matter candidates with masses decouple from the bath around T WIMP ≈ m 1 /x WIMP T QCD , late enough to neglect effects of the QCD phase transition on dark matter decoupling [51]. In i2DM this holds in the phases of coannihilation and partner annihilation. However, in the conversion phase, Eq. (4.1) is not necessarily satisfied because chemical and kinetic equilibrium are not guaranteed until freeze-out, see Sec. 3.2.3. In this case, the dark fermions decouple from chemical equilibrium at earlier times x 1,2 x WIMP , which might be affected by the QCD phase transition if the corresponding decoupling temperature is similar to T QCD . On the other hand, dark matter freeze-out is expected to be unaffected by the phase transition, because the dark states decouple at temperatures T T WIMP . In any case, bounds from laboratory searches exclude i2DM dark matter candidates with masses near 1 GeV, see Sec. 5. Even if non-equilibrium effects can change the decoupling times of the dark states, we do not expect effects from the QCD phase transition to affect the thermal history of dark matter candidates with masses well below the GeV scale.

Big Bang Nucleosynthesis
The formation of primordial light nuclei starts around T BBN ∼ 0.1 MeV. New physics can affect the nuclei abundances in multiple ways. First of all, new dark particles can affect BBN by modifying the Hubble rate or the entropy density of the universe. Second, efficient annihilation of MeV-GeV-scale dark particles to electrons or photons may alter the rate at which light elements form [52]. Third, dark particles decaying to electrons or photons at later times can destroy the already formed primordial nuclei [53,54]. In i2DM, dark matter annihilation and dark partner decays around T ∼ T BBN are not strong enough to cause observable effects. Late decays of dark partners with lifetimes τ 2 > t BBN , however, can destroy the newly formed elements through photodisintegration and constrain parts of the i2DM parameter space. Below we will discuss photodisintegration in i2DM in detail and briefly argue why annihilation and decays during BBN are not efficient.
For annihilating dark matter during BBN, a lower mass bound of m DM 10 MeV has been derived in [52]. This bound applies for vanilla WIMP annihilation with σv ≈ 10 −26 − 10 −28 cm 3 /s. In the coannihilation and conversion phases of i2DM, however, the cross section for annihilations χ 1 χ 1 → {e + e − , γγ} is much smaller than for vanilla WIMPs, resulting in weaker bounds on the dark matter mass. We expect these bounds to be superseded by constraints on ∆N eff from CMB measurements, see Sec. 4.4.
On the other hand, decays χ 2 → χ 1 + − and χ 2 → χ 1 γγ of dark partners with MeV-GeV masses and lifetimes τ 2 > 10 3 s can produce an electromagnetic cascade of photons with energies above the binding energy of light nuclei, consequently disintegrating them. 13 Disintegration is efficient only if the emitted photons do not lose their energy too rapidly before reaching the target nucleus. This condition is satisfied if the energy of the photons lies below the di-electron threshold, E th e + e − m 2 e /(22 T ) [55]. 14 On the other hand, for photodisintegration to take place, E th e + e − should lie well above the binding energy of light elements, corresponding to temperatures below a few keV. Photodisintegration is thus efficient for photon energies We calculate the effects of photodisintegration on the nuclei abundances starting from a continuous electron spectrum originating from χ 2 → χ 1 e + e − decays. We neglect loop-suppressed direct photon production, χ 2 → χ 1 γγ, as well as final-state radiation, resulting in a conservative bound on the electromagnetic flux obtained from partner decays [56]. To investigate photodisintegration for i2DM, we use the public code ACROPOLIS [54]. The code calculates the modified primordial abundances of light elements induced by photodisintegration, accounting for the redistribution of the energy injected by the decaying dark particle in the plasma. In particular, ACROPOLIS includes the exponential suppression in the photon spectrum for energies E γ > E th e + e − , which was 13 For lifetimes τ2 ∼ tBBN ∼ 10 2 s, we have estimated from the analysis of [53] that no further constraints arise. This expectation is justified by comparing our i2DM predictions with BBN bounds on decaying dark scalars with a certain lifetime, shown in Fig. 4 in [53]. Compared to [53], in i2DM we expect a smaller branching ratio to electrons and photons and softer spectra of the 3-body decays products; freeze-out temperatures significantly lower than 10 −2 GeV for MeV-GeV particles; and an exponentially suppressed abundance of the decaying dark partner χ2. 14 Qualitatively this condition can be understood from the requirement that the center-of-mass energy of the injected photon and the thermal bath photon scaling as EγEγ th is of the order of m 2 e , where Eγ is the energy of the injected photon and Eγ th ∼ T is the average energy of a photon from thermal bath. neglected in previous studies, 15 up to energies E γ < E 0 , where E 0 is a model-dependent upper limit on the photon spectrum. We set E 0 to the maximum kinematically allowed energy. For the initial abundances of light nuclei, we use the Standard BBN prediction extracted from the code AlterBBN [57]. The input number density of χ 2 is extracted from our system of Boltzmann equations.
The resulting constraints from photodisintegration on partner decays in i2DM are shown in Fig. 3. Bounds on the displayed parameter space are only visible in the lower right corner of the left panel with tan θ = 10 −4 , ∆ = 0.05, far from the correct relic abundance (the green dotted line). At larger tan θ, photodisintegration is not efficient because the dark partners decay too early. The lifetime of χ 2 also determines the upper edge of the excluded region. At dark matter masses, m 1 120 MeV, photodisintegration becomes inefficient due to the small absolute mass splitting between the dark fermions, m 2 − m 1 , which causes too soft decay products. For larger mass splittings ∆ = 0.1, in the right panel, photodisintegration is sensitive to smaller dark matter masses. However, the lifetime of χ 2 is generally smaller so that smaller couplings would be necessary for photodesintegration to take place. As a consequence, the excluded region lies below the plotted area in Fig. 3 and BBN bounds are irrelevant for viable i2DM relics. In summary, in i2DM photodisintegration excludes dark partners with lifetimes much larger than in cosmologically viable scenarios.

Cosmic Microwave Background
Dark sector particles that annihilate or decay around the time of recombination can affect the overall shape of the CMB black body spectrum, as well as its temperature and the polarization anisotropy spectra. Measurements of the CMB temperature and the polarization anisotropy spectra set strong constraints on the annihilation cross section of WIMP-like dark matter and also on decaying new particles with lifetimes τ > 10 13 s. For particles with shorter lifetimes, extra constraints can be obtained by studying deformations of the blackbody spectrum before recombination, usually referred to as spectral distortions. In addition, dark matter with couplings to neutrinos, photons or electrons can shift the effective number of neutrinos around the time of recombination ∆N eff (T CMB ), and affect the CMB anisotropies. We will discuss the constraints arising from ∆N eff in Sec. 4.4.
Charged SM particles that can arise from decays or annihilations of dark particles can inject energy into the plasma in the form of heat. Heat injections at redshifts z 2 · 10 6 induce spectral distortions in the CMB. Measurements of spectral distortions are sensitive to dark particles with lifetimes τ 10 4 s. 16 However, existing bounds from the COBE-FIRAS experiment [59] are largely superseded by the BBN bounds discussed in Sec. 4.2. 17 Future CMB missions similar to PiXie could strengthen the bounds on particles with 15 Previous studies used the universal photon spectrum, which is an analytic approximation working very well for energies Eγ < E th e + e − , but neglects photon with energies above this threshold. 16 In [58] it has been shown that spectral distortions cannot set competitive bounds on dark matter annihilation compared to bounds from CMB anisotropies. 17 See also [60] for FIRAS/Planck constraints on decays into low-energy photons.
lifetimes τ 10 4 s by up to two orders magnitude [61,62]. In contrast, BBN bounds are not expected to improve as much in the future.
Charged particles from decays or annihilations of dark sector particles can also ionize the plasma. Ionization at redshifts z 10 3 modifies the CMB anisotropy spectra. 18 This ionized fraction of the energy deposit induces a broadening of the surface of last scattering; for instance, it attenuates the CMB power spectrum on scales smaller than the width of the surface [64]. Measurements by the Planck collaboration constrain this effect and set strong upper limits on the cross section for s-wave dark matter annihilation [9]. In i2DM, dark matter annihilation χ 1 χ 1 → ff is suppressed well below these limits.
Compared to ionization effects on the CMB anisotropies, searches for dark matter annihilation in indirect detection experiments impose much weaker bounds in the MeV-GeV mass range [3,65]. Future missions such as e-astrogram [65] can provide competitive bounds on the annihilation cross section, which however are still far above the suppressed annihilation rates in i2DM.

Effective number of neutrinos N eff
As mentioned above, stable dark matter coupled to neutrinos, photons or electrons can change the effective number of neutrinos by ∆N eff . This modification can affect the abundance of light nuclei (set at T BBN ∼ 0.1 MeV) and the CMB anisotropies (set around T CMB ∼ 0.4 eV). Both measurements can thus set bounds on ∆N eff .
Efficient dark matter scattering with the electromagnetic or neutrino bath can induce an entropy transfer to these species after neutrino decoupling (at T ν ∼ MeV), thus modifying the neutrino-to-photon temperature ratio compared to the SM prediction. In practice, this information is encapsulated in ∆N eff . For dark matter that only couples to neutrinos, the entropy transfer reheats the neutrino bath and induces a positive shift ∆N eff > 0. For dark matter coupling to either electrons or photons, the electromagnetic bath gets reheated, which induces a negative shift ∆N eff < 0.
To estimate the effect of light dark particles on N eff , we follow the detailed analysis of [66], which uses the constraint on ∆N eff by the Planck collaboration [9], This result agrees with the estimates in [52,67], for instance. We emphasize that a modification of N eff can only be observed in CMB data if [67] (i) the dark matter is in kinetic equilibrium with either electrons, photons or neutrinos at temperatures above and below T ν ∼ MeV; 18 Any energy release into the plasma at redshifts earlier than z ∼ 1400 hardly affects the ionization history and has little impact on the CMB anisotropies [60,63]. As a result, MeV-GeV dark partners can only affect the CMB through spectral distortions.
(ii) the dark matter becomes non-relativistic at temperatures T T ν , typically for masses below a few tens of MeV.

As a result, in i2DM the bound from Eq. (4.4) only holds in regions (A) and (B), defined in
Eqs. (3.18) and (3.19), where the dark fermions are in kinetic equilibrium prior to freezeout. For simplicity, we restrict ourselves to dark matter candidates with m 1 > 10 MeV in all three regions.

Supernova cooling
Feebly interacting dark-sector particles can be produced in proto-neutron stars and freely escape them. The energy carried by the dark particles speeds up the supernova cooling and therefore can constrain models with light dark sectors. In i2DM, dark fermions are mostly produced through decays of dark photons produced in bremsstrahlung during neutronproton collisions or directly in collisions of light SM fermions [68]. Existing bounds on the cooling of the supernova SN1987A [69,70] constrain the i2DM parameter space only at very small couplings that are irrelevant for our purposes. Based on the results for iDM in [68], we estimate that in i2DM supernova cooling constrains dark couplings up to at most y ≈ 10 −13 for masses m 1 200 MeV; a parameter space where the dark matter relics would be overabundant. The bounds could potentially be even weaker if different core-collapse simulations were used [71]. Valid i2DM relics are thus not subject to bounds from supernova cooling.

Laboratory searches
In this section we discuss the phenomenology of inelastic dark matter in direct detection experiments, at particle colliders and at fixed-target experiments. Despite the original motivation of inelastic dark matter to evade direct detection, recent investigations reveal sensitivity to certain scenarios of iDM. Colliders and fixed-target experiments set strong bounds on the parameter space of iDM. We show that i2DM can evade some of these bounds, but can be conclusively tested at future experiments. Indirect detection searches are not sensitive to i2DM, because the pair-annihilation rate of dark matter is suppressed by small kinetic mixing and small dark fermion mixing, well below the reach of current and projected future experiments.

Direct detection
In scenarios of feebly interacting inelastic dark matter, elastic scattering χ 1 X → χ 1 X off nucleons or electrons, X = n or e, is suppressed below the sensitivity of current direct detection experiments. In the MeV-GeV mass range, the recoil energy in nucleon scattering is typically too small to be observed, but electron recoils are a promising road to detection [72]. The strongest current bound on dark matter-electron scattering by Xenon 1T lies around σ(χ 1 e → χ 1 e) ∼ 10 −40 cm 2 [73] for m 1 ≈ 100 MeV. For comparison, in i2DM with m 1 = 100 MeV, m A = 3m 1 and large mixings sin θ = 0.1 and = 0.02, the predicted cross section is σ(χ 1 e → χ 1 e) = 6 × 10 −41 cm 2 . Future direct detection experiments might reach a higher sensitivity to dark matter-electron scattering. Whether they can reach i2DM sensitivity will depend on the progress in collider searches (see Sec. 5.3), which can probe kinetic mixing down to ≈ 10 −3 and potentially suppress the target cross section as σ ∼ 2 .
If the mass splitting ∆ between the dark states is smaller than a few 100 keV, inelastic up-scattering χ 1 n → χ 2 n can produce an observable nuclear recoil signal, provided that the threshold of the experiment is low enough to detect the recoil energy [6,11].
For larger mass splitting, up-scattering can be observable if dark matter is accelerated through interactions with cosmic rays in the atmosphere around the earth [74]: The dark matter candidate χ 1 scatters inelastically off cosmic rays, mostly protons, via χ 1 p → χ 2 p. Subsequent decays χ 2 → χ 1 ff produce a relativistic component of dark matter, which passes the energy threshold for χ 1 n → χ 2 n scattering in experiments [22]. 19 Dark partners produced from cosmic-ray reactions can induce nuclear recoils via down-scattering χ 2 n → χ 1 n, provided that their lifetime is long enough to reach the experiment [21,22,76]. Down-scattering off electrons χ 2 e → χ 1 e is an interesting alternative, which can also address the current excess of electron recoils at Xenon1T [77]. For inelastic dark matter, efficient down-scattering only occurs if the lifetime is longer than several years [22], corresponding to mass splittings much smaller than those considered in this work.
For i2DM, we expect that the rates for up-scattering and down-scattering are generally suppressed compared to iDM, due to the smaller A χ 1 χ 2 coupling proportional to tan θ. On the other hand, in i2DM scenarios with small mass splitting, dark partners produced from cosmic-ray up-scattering are long-lived enough so that a substantial fraction of them could reach the detector before decaying. In this case, elastic scattering via χ 2 n → χ 2 n should dominate and leave an interesting characteristic signature of i2DM. In such scenarios, also down-scattering is expected. If the dark states are heavy and compressed enough to induce an observable nuclear recoil, up-scattering of i2DM off cosmic rays followed by elastic scattering off nuclei in the detector could be directly probed at experiments with a low energy threshold, cf. Ref. [21]. A dedicated analysis of i2DM at direct detection experiments goes beyond the scope of this work, but is a promising direction for future research.

Electroweak precision observables
A general bound on kinetic mixing of a dark photon is obtained from electroweak precision tests. In electroweak observables measured at LEP, Tevatron and the LHC, kinetic mixing modifies the Z boson's mass and couplings to SM fermions at O( 2 ). For dark photons with masses well below the Z resonance, a global fit to electroweak precision data yields a 95% CL upper bound of [78] 0.02 for m A 10 GeV; (5.1) stronger bounds apply for dark photons with masses closer to the Z pole. For m A 10 GeV, slightly stronger bounds have also been obtained from e ± scattering off protons at HERA [79].

Collider searches
At e + e − colliders, dark fermions coupling via a dark photon can be produced via three main processes: In models of inelastic dark matter, the first process is suppressed by construction. The second process dominates in iDM scenarios, which typically rely on the coupling of the dark photon to χ 1 and χ 2 . The third process is characteristic of i2DM, since the coupling A D χ 1 χ 2 is suppressed by tan θ and the dark photon mostly decays via A D → χ 2 χ 2 , see Eq. (2.9). Dark photon decays into SM fermions are suppressed as 2 α/α D , so that resonance searches at BaBar [80] and LHCb [81,82] are not sensitive to the scenarios investigated in this work. Collider signals of inelastic dark matter depend on whether the dark partners decay within or outside the detector. Below the hadronic threshold, the lifetime of χ 2 is mostly determined by χ 2 → χ 1 + − decays, which strongly depend on tan θ and the mass splitting ∆, see Eq. (2.14). If one or two dark partners decay within the detector, the signature consists of one or two prompt or displaced vertices of charged leptons, in association with a photon and missing energy. In iDM, the phenomenology of this signature has been investigated in detail for the Belle II experiment [16,18,83]. For sufficiently large mass splitting ∆, Belle II will be able to probe scenarios of iDM in the GeV range. For a smaller mass splitting, the decay products of χ 2 are too soft to be detected, leading to a signal with a photon and missing energy. The same signature is expected if χ 2 decays outside the detector. A search for mono-photon signals at BaBar has set an upper bound on the kinetic mixing of invisible dark photons [84], In iDM, this bound excludes most of the parameter space for dark matter candidates below the GeV scale. A similar search at Belle II can probe even smaller dark sector couplings and thereby improve the sensitivity to iDM [16]. In i2DM, the dark photon decays close to its production point due to efficient A → χ 2 χ 2 decays. However, the dark partners have larger decay lengths than in iDM, because χ 2 → χ 1 ff decays are suppressed by tan 2 θ, see Eq. (2.14). Therefore the dark photon does not leave any trace in the detector and the bound of Eq. (5.3) from BaBar's monophoton search applies. In Fig. 3, we display the bounds on kinetic mixing from collider searches in the parameter space of i2DM for ∆ = 0.05 (left) and 0.1 (right). The observed relic abundance is obtained along the contours for fixed values of the dark fermion mixing

Bounds from fixed-target experiments
Fixed-target experiments with a large separation of the particle source and the detector are particularly sensitive to particles with a long decay length. Searches for long-lived particles with sub-GeV masses have been performed at various fixed-target experiments and have been reinterpreted for iDM scenarios, for instance in Refs. [12,15,85]. In general, a beam of particles is dumped on a target material, producing light mesons such as pions or kaons. In iDM and i2DM, dark photons can be produced either in meson decays, through the Primakoff process or from bremsstrahlung. Subsequently the dark photons decay into pairs of dark fermions χ 1 and/or χ 2 .
Depending on the model parameters, inelastic dark matter can be detected via three signatures: (displaced) decays of dark partners, χ 2 → χ 1 + − ; (up-)scattering of dark fermions off the detector material, χ i N → χ j N ; or missing energy from dark fermions that are stable at the scales of the experiment. The relative sensitivity to each signature depends on the lifetime of the dark partner, as well as on the experimental setup: shortlived dark partners are mostly observed through decays, while long-lived partners scatter inside the detector material or decay after passing through the detector. In what follows, we discuss the relevant signatures for i2DM and derive bounds from searches at fixed-target experiments.
Partner decays For partner decays, the expected event rate in the far detector is given by 20 4) where N A ∝ 2 is the total number of dark photons produced in a given experiment, and N 2 is the total number of χ 2 states resulting from N A dark photon decays. The branching ratios are B(A → χ 1 χ 2 ) ≈ 1 for iDM and B(A → χ 2 χ 2 ) ≈ 1 for i2DM with tan θ 1. 21 Finally, P dec (d k ) is the probability to detect the decay products of particle k with decay length d k . The decay length d k = (βγ) k cτ 2 depends on the boost, (βγ) k , and on the lifetime, τ 2 , of the dark partner. For decay lengths longer than the baseline, the decay probability scales as P dec (d) ∝ 1/d. In this regime, the expected event rates for iDM and i2DM depend on the model parameters as In i2DM, the lifetime of the dark partner scales as τ 2 ∼ 1/ tan 2 θ. For small tan θ, the dark partner tends to decay after passing the detector, which reduces the event rate. The factor of 2 accounts for the two dark partners produced in A → χ 2 χ 2 decays. For fixed parameters α D , m 1 , m A , ∆, a bound on obtained from searches for dark partner decays in iDM translates into a bound on · (2 tan 2 θ cos 4 θ) 1/4 in i2DM.
Decays at CHARM Strong bounds on long-lived dark particles decaying into electrons have been set at the neutrino experiment CHARM. At CHARM, dark fermions can be efficiently produced from π 0 or η meson decays, π 0 (η) → γA → γχ i χ j . Null results of a beam-dump search for heavy neutrinos decaying into electron pairs [86] have been reinterpreted for χ 2 → χ 1 e + e − decays in iDM [15]. In Fig. 4, we show the resulting bounds in the parameter space of i2DM, using Eq. (5.5) to rescale the predicted event rates. The three scenarios are distinguished by the mass splitting ∆, while we have fixed = 10 −3 to evade the collider bounds from Sec. 5.3. The observed relic abundance is obtained along the colored contours. The lifetime of χ 2 varies strongly with the mass splitting, τ 2 ∝ 1/∆ 5 , see Eq. (2.14). For ∆ = 0.1, most of the dark partners decay within the considered decay volume and the search is sensitive to even small mixing tan θ. For ∆ = 0.05, the sensitivity decreases due to the longer lifetime and softer e + e − momenta, which are less likely to pass the analysis cuts. For ∆ = 0.01, the dark partner is essentially stable compared to the length of the decay volume and the search becomes insensitive to i2DM. 20 Here we neglect dark partners produced via upscattering χ1 → χ2. 21 In our numerical analysis of i2DM, we neglect A → χ1χ2 decays. Decays at LSND For light dark sectors, even stronger bounds have been obtained from the neutrino experiment LSND. At LSND, dark fermions can be abundantly produced from pion decays via π 0 → γA → γχ i χ j . In Ref. [12], LSND data has been interpreted in terms of χ 2 → χ 1 e + e − decays in iDM, under the conservative assumption that the e + e − pair is not resolved in the calorimeter and can mimic elastic neutrino-electron scattering. In Fig. 4, we show the resulting bounds rescaled for i2DM (labelled 'LSND decay'). LSND is very sensitive to dark partners with m 2 < m π 0 /2. The lower cutoff is determined by the kinematic threshold for χ 2 → χ 1 e + e − decays. The sensitivity could be improved with a dedicated analysis of three-body decays, rather than a re-interpretation of neutrino-electron scattering.
Dark fermion scattering In addition to decays, long-lived dark fermions can be detected through up-scattering χ 1 → χ 2 , down-scattering χ 2 → χ 1 , or elastic scattering χ 2 → χ 2 inside the detector. We neglect up-scattering, which is typically sub-dominant for suppressed χ 1 − χ 2 couplings. The expected scattering rate is then given by where P scat (d k ) is the probability for particle k to scatter off the material inside the detector, with d k → ∞ for χ = χ 1 , and N 2 is the number of produced dark partners. Neutrino experiments are particularly sensitive to dark fermion scattering, which mimics neutrinoelectron scattering. Among various experiments, LSND sets the currently strongest bounds on inelastic dark matter. For iDM and sub-GeV masses, the relevant process is downscattering χ 2 N → χ 1 N [12]. For i2DM, elastic scattering χ 2 → χ 2 dominates for tan θ 1. The respective event rates scale as for decay lengths d k larger than the distance between the target and the detector.
Scattering at LSND The LSND bounds on scattering are derived from the same analysis as for partner decays. Again, we translate the results for iDM from Ref. [12] to i2DM, shown in Fig. 4 as 'LSND scat'. For tan θ 0.3, the bounds are insensitive to dark fermion mixing and exclude small dark matter masses. The bounds disappear for small kinetic mixing 10 −4 , where LSND loses its sensitivity due to the low dark photon production rate.
It is interesting to compare these results for i2DM with iDM. For tan θ = 1, the A χ 1 χ 2 interaction strength is similar in iDM and i2DM. For the benchmark scenarios shown in Fig. 4, iDM is excluded by CHARM and LSND, unless partner decays into electrons are kinematically forbidden. Viable scenarios of sub-GeV iDM require a stronger interaction strength y for efficient coannihilation, while keeping the kinetic mixing small to evade bounds from mono-photon searches [20], see Sec. 5.3. In turn, i2DM scenarios with small tan θ can evade the CHARM bounds due to the longer lifetime of the dark partners. In this regime, the relic abundance is set by partner annihilation and/or coscattering. In summary, null searches at current fixed-target experiments are a severe challenge for iDM, while i2DM is a viable option due to the impact of partners on dark matter freeze-out.

Prospects of future fixed-target experiments
As we discussed in Sec. 5.4, existing fixed-target experiments are very sensitive to sub-GeV inelastic dark matter. However, in the multi-GeV range i2DM scenarios with a compressed dark sector currently escape detection. In order to fully explore the parameter space of i2DM, we study the discovery potential of proposed searches at fixed-target experiments, which could be realized in the near future. While a number of experiments can be sensitive to i2DM, here we focus a few promising proposals. At SBN, similarly to other fixed-target experiments discussed in Sec. 5.4, dark partners can be produced from decaying dark photons, which are created in the interaction of the proton beam with the target. The SBN detectors can be used to search for χ 2 → χ 1 + − decays or scattering χ i → χ j off the detector material.
Dark sector searches at neutrino experiments inevitably feature a large neutrino background. At SBN, there are two proposals to reduce this background. One option is to deflect the proton beam around the target into an iron absorber placed 50 m downstream, which has previously been done at MiniBooNE to study light dark sectors [88,89]. This mode is referred to as "off-target". The second option is to use the NuMI 120 GeV proton beam, impacting on a graphite target. ICARUS and MicroBooNE are placed at angles of 6 and 8 degrees against the NuMI beam direction. This option is known as "off-axis". Both off-target and off-axis options have been studied for iDM [20], where SBND has the best sensitivity in the off-target mode, while ICARUS can set the strongest bound in the off-axis mode. In Fig. 4, we show our reinterpretation of these predictions for i2DM, following the procedure described in Sec. 5.4. We present the results for ICARUS, assuming that the off-axis mode can be realized with much less technological effort than the off-target option.
Compared to existing fixed-target experiments, ICARUS is substantially more sensitive to i2DM, provided that the dark partners are sufficiently long-lived to induce enough signal in the detector.
SeaQuest Originally developed to study the sea quark content of the proton with a 120 GeV proton beam and various targets, the Fermilab experiment SeaQuest has a good potential to probe dark sectors [13,90]. It is already equipped with a displaced muon trigger to study exotic long-lived particles decaying to muons, and could be supplemented by an electromagnetic calorimeter to also probe electron signals. In Ref. [13], the sensitivity of SeaQuest to dark partner decays in iDM has been studied for three different decay volumes. In Fig. 4, we show the corresponding predictions for i2DM for the largest possible decay volume. Compared with CHARM and ICARUS, SeaQuest can improve the sensitivity to i2DM for dark matter with masses near the η resonance. As discussed in Ref. [13], the reach of SeaQuest could be further enhanced by running the experiment without the magnet, which however would require a dedicated analysis of the experimental setup.
LDMX The proposed Light Dark Matter eXperiment (LDMX) [91] is an electron beam-dump experiment designed primarily for probing light dark matter models. Its search strategy relies on dark sector particles being produced in the beam dump that escape the detector, which extends up to about 1 m downstream from the target. This gives rise to a signature of missing energy. All charged particles in an event are vetoed, except for the soft remnant of the incoming electron. A signal of inelastic dark matter is detected if a substantial fraction of dark partners decay after passing through the detector. The LDMX collaboration has investigated the projected sensitivity to many MeV-GeV dark sector models, including iDM [14]. In Fig. 4, we show our interpretation for i2DM for the most conservative design option. LDMX is well suited to probe i2DM scenarios with very small mass splitting ∆, which are difficult to detect in experiments that detect the decay products of the dark partner.
We summarize our projections for future fixed-target experiments for the three i2DM benchmarks in Fig. 4. For ∆ = 0.1, the dark matter target is already probed by existing experiments. We therefore do not show projections for future experiments, but note that they could be sensitive to other regions of the parameter space. For ∆ = 0.05, SeaQuest alone can improve the sensitivity to dark partner decays, but cannot fully probe the dark matter target. ICARUS can complement SeaQuest by also detecting dark fermion scattering, which is predominant for small dark matter masses. LDMX, searching for missing energy, can extend the reach of ICARUS at larger masses and long lifetimes. Either ICARUS or LDMX could conclusively probe this scenario. For even smaller mass splitting ∆ = 0.01, dark partner decays cannot be observed due to the long lifetime and the softness of the SM decay products. The sensitivity of ICARUS and LDMX through scattering and missing energy, however, is kept and allows to conclusively probe the dark matter scenario. Larger dark photon masses m A > 3m 1 or smaller kinetic mixing reduces the rate of produced dark photons and thus the sensitivity of any experiment. For smaller dark couplings α D , the lifetime of the dark partners is enhanced and dark partner decays close behind the target are less abundant. In this case, scattering and missing energy are more promising signals, especially for small dark matter masses. All in all, fixed-target experiments have a high potential to conclusively test i2DM in the near future, provided that they are built and successfully run.

Conclusions and outlook
In this work we have introduced a new model for feebly coupling dark matter, called inelastic Dirac Dark Matter. Compared to the widely studied model of inelastic Dark Matter with Majorana fermions, i2DM has a different thermal history and is less constrained by current laboratory searches. The main difference is due to the variable interaction strength of the dark matter candidate χ 1 , parametrized by a mass mixing θ with its dark partner χ 2 .
At small mass mixing, the dark matter candidate decouples from the SM bath before freeze-out and the relic dark matter abundance cannot be set by dark matter annihilation, coannihilation or partner annihilation anymore. Instead, coscattering and decay processes are crucial to explain the observed abundance even for tiny dark matter interactions y with the thermal bath. At such feeble couplings, dark matter can decouple from chemical and even kinetic equilibrium before freeze-out. We have computed the resulting effects on the relic abundance by numerically solving a coupled set of Boltzmann equations for the time evolution of the dark fermions. For particles in the multi-MeV range, we find that the relic abundance can be obtained with interactions as feeble as y 10 −11 and is only mildly affect by deviations from kinetic equilibrium.
Requesting that the QCD phase transition should not affect the freeze-out dynamics and that ∆N eff contributions for i2DM is in agreement with CMB data be suppressed, we identify the cosmologically viable parameter region for i2DM candidates with mass and interaction strength in the range 10 MeV < m 1 < 1 GeV , y > 10 −11 . (6.1) In this region, we have investigated possible effects of i2DM on Big Bang Nucleosynthesis, the Cosmic Microwave Background and supernova cooling, but find them much too small to modify current observations. Laboratory searches for i2DM mostly rely on the production of dark photons, subsequently decaying into dark fermions. Due to the feeble interaction, the dark partner typically appears stable at the scales of current colliders. Searches for mono-photons and missing momentum at flavor experiments set a strong bound on the overall coupling of the dark photon to the Standard Model, 10 −3 . These bounds imply a minimum dark fermion mixing and require efficient partner annihilation to satisfy the relic abundance. The parameter range for viable i2DM is thereby confined to light and compressed dark sectors with At direct detection experiments, χ 1 N → χ 2 N up-scattering off atomic nuclei is not efficient enough to further test i2DM, even if dark matter can be accelerated through interactions with cosmic rays. On the other hand, long-lived dark partners can be produced in cosmicray up-scattering might leave an observable signature of χ 2 N → χ 2 N elastic scattering in the detector material. We leave this interesting avenue for future work. Very promising probes of i2DM are fixed-target experiments with a long baseline, which can search for decays or scattering of long-lived dark partners in far detectors. Scattering is a prominent signal for light dark partners, while decays can be observed for heavier partners up to the GeV scale. Current searches for such signatures at CHARM and LSND already probe a significant portion of the i2DM parameter space, but lose steam for small mass splitting, where the momentum deposit in the detector is soft.
In the near future, neutrino experiments like ICARUS and SeaQuest or the beam-dump experiment LDMX can significantly improve the sensitivity to inelastic dark matter. Together with mono-photon searches at Belle II, they can conclusively test if inelastic Dirac dark matter is at the origin of the observed relic abundance.

A Formalism of i2DM
We provide further details about inelastic Dirac dark matter introduced in Sec. 2. The interactions of mass eigenstates described by the Lagrangian from Eq. (2.1) can be derived from gauge-invariant interactions of the dark sector with the Standard Model in the unbroken phase of the underlying theory, Here D µ = ∂ µ + ig D Q DÂDµ is the covariant derivative for the gauge fieldÂ D of the dark U (1) D symmetry. The field strength tensors for the dark and the hypercharge gauge fields are denoted asF µν D andB µν respectively. The dark fermion field χ D is charged under U (1) D , while χ 0 is a singlet under all gauge interactions.
The scalar φ D with potential V (φ D ) and charge Q D = 1 facilitates spontaneous U (1) D breaking, once it develops a vacuum expectation value v D , so that φ D = (v D + ϕ D )/ √ 2. Upon symmetry breaking, the dark photonÂ D acquires a mass mÂ D = g D v D , and the dark fermions χ D and χ 0 mix through the dark Yukawa coupling λ, resulting in two mass eigenstates The corresponding masses are where ∆ is the relative mass difference defined in Eq. (2.4). In order to obtain canonical kinetic terms in Eq. (A.1), we first redefine the U (1) Y and U (1) D gauge bosons to absorb the kinetic mixing via whereˆ = / cos θ W . Subsequently, we transform the fields to the physical eigenstates of weak interactions, A µ , Z µ , A µ , via two rotations: one rotation, R ξ , mixing the SU (2) gauge field W 3 µ andÃ µ D with an angle ξ; and a second rotation, R W , mixing W 3 µ andB µ with the Weinberg angle θ W . The overall transformation takes the form and mẐ = gv/(2c W ), with the SU (2) gauge coupling, g, and the SM Higgs vacuum expectation value, v. The masses of the Z boson and the dark photon are while the photon remains massless. In the limit of small kinetic mixing, { , η, ξ} 1, one recovers m Z mẐ and m A mÂ D [78,[92][93][94]. To derive the Feynman rules and transition amplitudes, we use FeynRules [95]. In particular, we have adapted the Hidden Abelian Higgs Model [92], which is available from the FeynRules model data base to i2DM. Our FeynRules model for i2DM is available at [33]. In terms of mass eigenstates, the Lagrangian finally reads where t W refers to the tangent of the Weinberg angle θ W . The dark currents J µ i with (i = 1, 2) and J µ 12 are defined in Eq. (2.2). The SM fermion currents are given by where f are the SM fermions and As a result, the photon has no fundamental couplings to the dark fermions χ 1 , χ 2 . In particular, the dark matter candidate χ 1 carries no millicharge and is not subject to otherwise strong constraints [96].

B Boltzmann equations beyond thermal equilibrium
On general grounds, in an isotropic and homogeneous universe the evolution of a particle species i is described in terms of a distribution function f i (t, | p i |), expressed in terms of the physical time t and of the norm of the physical 3-momentum p i , denoted as | p i | ≡ p i . 22 If a species i interacts with other species j, the time evolution of species i is described by the Boltzmann equation In the text, pi usually refers to the norm of the 3-momentum, except when it appears in a 4-dimensional delta function which enforces both 3-momentum and energy conservation.
where C[f i (t, p i ), f j (t, p j )] is the collision term involving all decay and scattering processes with the other species j, and E i = p 2 i + m 2 i is the energy of particles of species i with mass m i . In this appendix, and in particular in Sec. B.2, we give more details about the collision terms for coscattering, C coscat , and 3-body inverse decays, C decay .
The total time derivative in Eq. (B.1) can be re-expressed in terms of partial derivatives with respect to time and momentum as In Eq. (B.2) we have used the rescaled time x and momentum variable q i introduced in Eq. (3.3). In this appendix, we will make use of the variables (t, p i ) or (x, q i ) whenever convenient. Also, H = H(x) is the Hubble expansion rate and h eff = h eff (x) is the effective number of relativistic degrees of freedom contributing to the entropy density s = h eff 2π 2 /45 T 3 . 23 In this paper we consider dark matter production in a radiation dominated era, in which case the Hubble rate reduces to where M P is the Planck mass and g * = g * (x) denotes the number of relativistic degrees of freedom in the thermal bath at time x contributing to the radiation energy density ρ R = g * π 2 /30 T 4 . For our numerical analysis, we use tables available in the public code micrOMEGAs [44] to evaluate h eff (x) and g * (x). When the number of relativistic degrees of freedom, h eff and g * , can be considered constant around freeze-out, the unintegrated Boltzmann equations of Eq. (B.2) simplify to as reported in Eq. (3.8). To compute the dark matter freeze-out beyond kinetic and chemical equilibrium in i2DM, we will use the latter evolution equation, see Sec. B.2 for details.
In what follows, we will encounter the equilibrium number densities n i eq (t), i.e., the number densities obtained by integrating over the kinetic equilibrium distribution functions f eq i (t, p i ), where f eq i (t, p i ) are the Fermi-Dirac, Bose-Einstein or Boltzmann distributions for zero chemical potential. 23 The insertions of dh eff /dx in Eqs. (B.2) and (B.3) are due to the usual choice of rescaled momentum q = p/T , instead of the time-independent comoving momentum q = pa, where a is the scale factor. Another convenient choice of time-independent rescaled momentum would be q = p/s 1/3 see e.g. [97,98].

B.1 Deviation from chemical equilibrium before freeze-out
In i2DM, as long as deviations of χ 1 from kinetic equilibrium with the heat bath can be neglected up until freeze-out, the evolution of χ 1,2 can be described by the number densities n 1,2 (x), or equivalently by the comoving number densities Y 1,2 (x). Within this framework, it is useful to introduce the reaction densities for decays and scattering, which can be rewritten in terms of the thermally averaged cross section and the decay rate in the rest frame of species i, as in Eq. (3.11). In Eq. (B.6), |M| 2 is the squared scattering amplitude, averaged over initial-and final-state degrees of freedom, and dφ i = g i d 3 p i /(2E i (2π) 3 ) denotes the phase-space element.
As long as χ 1,2 can be assumed to be in kinetic and chemical equilibrium before freeze-out, the set of partial integro-differential equations for the dark fermions from Eq. (B.1), integrated over the 3-momenta p i , reduces to one single well-known ordinary differential equation for the total dark sector number density n(t) = i=1,2 n i (t) [23,41]. This simplification applies to freeze-out in the phases of coannihilation and partner annihilation in region (A). In this regime, one can compute the dark matter relic abundance with any publicly available Boltzmann solver, such as micrOMEGAs [44], DarkSUSY [45] or MadDM [46].
In the conversion phase, the very small interaction strength y implies strongly suppressed coscatterings and decays and chemical equilibrium between χ 1 and χ 2 is lost. As a result, the relation n 1 /n 2 = n eq 1 /n eq 2 cannot be assumed. In order to take deviations from chemical equilibrium into account, we have to solve a coupled system of Boltzmann equations including the conversion rates γ 1→2 [25,27,28,30,31], 24 where γ ij are the reaction densities defined in Eq. (3.13) and γ i→j = γ dec i→j + γ scat i→j . For i2DM, we have checked that we can neglect quantum statistical effects in the relevant density distributions f i and use the Maxwell-Boltzmann equilibrium distributions f eq i . In i2DM, Eq. (B.7) can safely be used for freeze-out scenarios in regions (A) and (B), where kinetic equilibrium is maintained up to freeze-out. 25 24 Similar Boltzmann equations exist for the antiparticlesχ1 andχ2. In i2DM these equations are completely equivalent and we obtain the total abundance for χ1 +χ1 by simply doubling the yield for χ1. 25 For the scenarios considered in this paper, the results we obtain by carefully accounting for kinetic decoupling, see Sec. B.2, are very close to those derived using Eq. (B.7), see the discussion around Fig. 1.

B.2 Deviation from chemical and kinetic equilibrium around freeze-out
When the conversion rate is around 10 to 100 times the Hubble rate at the time of χ 2 chemical decoupling, dark matter can no longer be expected to be kept in chemical or kinetic equilibrium with χ 2 around the DM freeze-out time x fo . The assumptions for Eq. (B.7) do not hold anymore, and a priori one needs to solve Eq. (B.5) including the full momentum and time dependentce in the χ 1 distribution function f 1 (x, q 1 ). For i2DM, we apply two levels of simplification to Eq. (B.5). First, we neglect the time dependence in h eff for the dark matter density evolution, setting h eff = h eff (x fo ) when integrating over f 1 (x, q 1 ). For dark fermions in the MeV-GeV mass range, we obtain a good approximation to the exact solution of Eq. (B.5). In this mass range, freeze-out is expected to happen between the QCD phase transition and neutrino decoupling, a period in which h eff 10.75 is fairly constant. Second, we follow Refs. [25,26] and only take into account the dominant interaction processes driving the dark matter distribution f 1 (x, q 1 ) towards kinetic equilibrium. For i2DM, this means that we include conversion processes but neglect coannihilation.
With these simplifications, the unintegrated Boltzmann equation for f 1 (x, q 1 ) from Eq. (B.5) can be rewritten as where q = q 1 = p 1 /T to simplify the notation. The contributions to the collision operator C 1→2 (x, q) =C coscat (x, q) +C decay (x, q) from coscattering and decays are spelled out in Secs. B.2.1 and B.2.2. The above description assumes that χ 2 and the light SM fermions involved in the conversion processes are in kinetic equilibrium. We also neglect all spin statistics effect. The differential equation in Eq. (B.8) can be solved iteratively with Y 2 (x) as an input. The latter is obtained from the integrated Boltzmann equation for Y 2 in Eq. (B.7), which in turn involves Y 1 (x), or equivalently the zeroth moment of f 1 (x, q) in q, obtained by integrating over Eq. (B.8). More details on the integration of Eq. (B.8) will be discussed in Sec. B.2.3.
In the second equation we have used the relation of detailed balance, f eq 2 f eq f = f eq 1 f eq f . With Eq. (B.10), we can write the collision term as with the coscattering collision operator Here the cross section for coscattering process is defined as [27,28] σ coscat (s) = 1 with the modulus of the dark matter 3-momentum in the centre-of-mass frame,p 1 , the squared center-of-mass energy, s, and the energy variables (B.14)

B.2.2 Three-body decays
As mentioned in Sec. 3.1, the inverse decay process χ 1 f f → χ 2 can play a role in keeping χ 1 in kinetic equilibrium. The corresponding collision reads Assuming again that all involved SM fermions f, f are in chemical and kinetic equilibrium with the thermal bath and that χ 2 is in kinetic equilibrium with the bath, the collision term reduces to 1 E 1 C decay =C decay f eq 1 n 2 n eq 2 − f 1 , (B. 16) with the collision operator for decays, after integrating over dφ f . Applying the methodology from [99] for 3-body decays and neglecting the SM fermion masses, the collision operator can be rewritten as where λ(x, y, z) = (x 2 − (y + z) 2 )(x 2 − (y − z) 2 ), (B.19) The integration boundaries are E ± 2 = m 2 2 + (p ± 2 ) 2 , with (B.20) h(x, q) = f eq 1 (x, q) .
Multiplying both sides of Eq. (B.21) by u(x, q) = exp dx g(x, q) and using ∂u(x,q) ∂x = g(x, q)u(x, q), we can simplify this equation to ∂(u(x, q)f 1 (x, q)) ∂x = ∂u(x, q) ∂x h(x, q). where we set the initial time to x 0 = 1. In i2DM, at early times we can use the boundary conditions that χ 1 is in kinetic equilibrium, f 1 (x 0 , q) = f eq 1 (x 0 , q), and that χ 2 is in chemical equilibrium, Y 2 (x 0 ) = Y eq 2 (x 0 ), such that h(x 0 , q) = f 1 (x 0 , q). The dark matter phase-space distribution then reduces to f 1 (x, q) = f eq 1 (x, q) In order to solve this equation, we have to specify the comoving number density of χ 2 , which remains in kinetic equilibrium throughout the whole evolution of χ 1 . Hence, we can use the integrated Boltzmann equation for χ 2 from Eq. (B.7). Solving Eqs. (B.25) and (B.7) together is numerically difficult, so we choose to solve the two equations iteratively. In a first step, we solve the system of integrated Boltzmann equations from Eq. (B.7). This gives us an initial value for Y 2 (x). We then solve the unintegrated Boltzmann Eq. (B.21) for χ 1 and feed again the integrated Boltzmann equation for χ 2 from Eq. (B.7) to obtain the next iteration for Y 2 (x). We stop this iteration once the difference in the relic dark matter abundance between the last two iterations is less than one percent.