Dark Matter abundance via thermal decays and leptoquark mediators

We explore a new mechanism for reproducing the Dark Matter (DM) abundance: scatterings of one DM particle on light Standard Model particles. Strong bounds on its decays can be satisfied if DM undergoes freeze-in and has a mass around or below the pion mass. This happens, for example, in theories with a right-handed neutrino interacting with charged fermions through a leptoquark exchange. These leptoquarks can be linked to the ones motivated by the B-physics anomalies if assumptions about the flavour structure are made. DM signals are unusual, with interesting possibilities for direct and indirect detection. Achieving thermal freeze-out instead requires models with more than one DM flavour, and couplings parametrically smaller than what needed by the usual pair annihilations.


Introduction
The measured value of the Dark Matter (DM) energy density can be understood in terms of various hypothetical scattering processes occurring in the early universe. The most studied scenario considers 2 ↔ 0 annihilation processes, where we count the number of DM-sector particles plus anti-particles in the initial and final state, without indicating the Standard Model (SM) particles. Examples of such processes that change the number of DM A Z 2 symmetry acting on the dark sector only usually justifies the presence of 2 DM (or, more generically, dark sector) particles and implies DM stability. These annihilations (or, more generically, co-annihilations) can lead to DM around the TeV scale as a freeze-out thermal relic. A freeze-in non-thermal relic and other possibilities can also arise through the same interactions.
Processes more complicated than 2 ↔ 0 received less attention: 2 ↔ 1 scatterings, known as semi-annihilations, can arise in models where DM is stable, and behave similarly to annihilations (see [1][2][3] for recent references). 3 ↔ 2 or 4 ↔ 2 scatterings give qualitatively different behaviours sometimes dubbed 'cannibalistic DM' [4] or 'pandemic DM' [3], depending on the regime and on whether DM is in thermal equilibrium with the SM sector. These scatterings can reproduce the DM cosmological density for DM masses in the eV or MeV range, and can be justified by Z 3 or higher symmetries.
Little attention has been received by the simplest possibility: 1 ↔ 0 processes that could be dubbed 'decadent DM'. Examples of such transitions that change the DM number by one unit are DM SM ↔ SM SM, DM ↔ SM SM, SM ↔ DM SM.
These processes arise in absence of a symmetry that forbids couplings between one DM particle and the SM sector. This is what can happen, for instance, in presence of interactions under which both DM and SM fields are charged, where imposing a symmetry that separates the dark sector from the SM might be not be a viable option. Prototypical examples of such a scenario include "sterile neutrino" DM, that can couple to leptons via a Higgs boson or a right-handed vector boson, or models where DM interacts with quarks via a leptoquark mediator. In both cases, the mediator carries a conserved SM charge and is allowed to couple to a DM-SM current. The latter case of DM interactions mediated by a leptoquark is particularly interesting in the context of models that aim at explaining the anomalies observed in B-meson decays (see e.g. [5][6][7][8][9][10][11][12] and references therein). The possibility of DM coupled to leptoquarks has been considered before [13][14][15][16][17], but usually in combination with an additional Z 2 symmetry that forbids DM decay. Here we shall explore the minimal possibility that allows the mediator to simultaneously couple to DM-SM and pure SM currents.
Given that 1 ↔ 0 processes are unavoidable in some theories, we show that, on the other hand, they can play a useful role in cosmology. Indeed, since DM needs to scatter on SM particles which are abundant in the primordial plasma, in order to reproduce the DM relic abundance these processes need a smaller cross section than that of the usual DM annihilations in eq. (1), where DM needs instead to find a rare DM particle. We refer to the interactions in eq. (2) as "thermal decays" because at finite temperature such transitions are conveniently described by Boltzmann equations with a thermal width Γ T of DM. This terminology emphasises the obvious and potentially insurmountable problem: DM risks decaying too fast, since there is no separate dark sector with a conserved DM number. Interactions that lead to the thermal decays in eq. (2) also lead to DM decays. Let us be quantitative: • On one hand, the DM life-time τ DM = 1/Γ DM must surely be longer than the age of the Universe τ DM > T U (∼ 10 18 s). Many DM decay modes are more strongly constrained, requiring τ DM > ∼ 10 26 s [18][19][20].
• On the other hand, thermal decays can give non-relativistic DM freeze-out 1 reproducing the cosmological DM abundance if their rate Γ DMT is comparable to the Hubble rate H ∼ √ d SM (T / GeV) 2 /(5 10 −6 s) at temperatures around the DM mass M , where d SM is the number of SM degrees of freedom.
As a result, the two requirements conflict easily and by many orders of magnitude: one needs Γ DM /Γ DMT < ∼ 10 −30 .
This gap by 30 orders of magnitude is such a strong constraint on possible models that this scenario is widely ignored. We nevertheless explore this possibility: while the simplest freeze-out picture in models with one DM particle is clearly ruled out, more complex possibilities with multiple DM particles or with DM freeze-in emerge and can lead to a peculiar phenomenology. A common feature of all viable models is the presence of some small coupling between DM and SM, which could be related to an underlying flavour structure.
The paper is organized as follows. In section 2 we write the generic form of Boltzmann equations that describe DM with thermal decays and study their features, deriving approximate solutions. We find that the DM relic abundance depends on a different power of the cross-section than usual DM annihilations. As a result DM thermal decays can be more efficient than DM annihilations, and can match the cosmological DM abundance with rates not much stronger than the minimal ones needed to reach thermal equilibrium. DM freeze-in needs smaller sub-equilibrium rates.
In section 3 we discuss possible models that only involve one DM particle. DM has to be lighter than hundreds of MeV in order to close the phase space for tree-level decays into SM particles. We then focus on models where DM carries lepton number (e.g. a right-handed neutrino) and interacts with the SM via TeV-scale leptoquarks motivated by the flavour anomalies. We will find that such models can reproduce the DM abundance via freeze-in, but not via freeze-out.
In section 4 we study possible signatures. We first explore possible values of the DM couplings motivated by simple flavour patterns, and the possible relations with B-physics anomalies. We then discuss direct and indirect detection. Interactions involving one DM particle give the unusual direct detection signals discussed in section 4.2, suggesting a tentative explanation of the Xenon1T anomaly. Indirect detection signals are discussed in section 4.3, where we show that DM decays DM → νγ could explain the line claimed at 3.5 keV, compatibly with cosmological production via freeze-in of thermal decays.
Finally, in section 5 we show that freeze-out is compatible with DM stability in models with two or more DM particles (for example three generations of right-handed neutrinos). 1 Freeze-in needs a smaller Γ DMT . We do not consider relativistic DM freeze-out (analogous to neutrino decoupling), because the resulting DM abundance would be of order unity Y DM ≡ n DM /s = 45ζ(3)g DM /g SM (T dec )π 4 and thereby the DM cosmological abundance would be reproduced for warm DM with mass M ≈ 0.7 eV g SM (T dec )/g DM . To avoid problematic warm DM, the DM decoupling temperature T dec must be mildly below the DM mass M . So we need a sizeable Γ DMT ∼ H at T ∼ M .
Conclusions are given in section 6.

Boltzmann equation for thermal decays
As usual, it is convenient to write the Boltzmann equation for Y ≡ n/s in terms of z = M/T , where M is the DM mass, n is the DM number density, and s the total entropy. Including scatterings that involve two DM (the usual DM annihilations) and one DM (our thermal decays) the Boltzmann equation is where the factor Z = 1/(1 + 1 3 d ln d SM d ln T ) can often be approximated as 1. In order to match the DM cosmological density, DM must freeze-out while non-relativistic, so that its thermal equilibrium abundance is Boltzmann suppressed, n eq i d i (m i T /2π) 3/2 e −m i /T , where d i is the number of degrees of freedom for the particle i. The space-time density of scatterings γ eq i can thereby be simplified going to the non-relativistic limit. The initial state of the usual DM DM ↔ SM SM annihilations becomes non-relativistic at low temperature. When considering more generic processes an unusual step is needed: identifying which side of the reaction becomes non-relativistic, namely the side that involves heavier particles. Defining as A, B the particles that can scatter at rest (they can be in the initial state or in the final state, depending on the kinematics of each process), the rate is approximated as The initial-state symmetry factor S is only present when the non-relativistic particles are in the initial state. In the presence of multiple processes, the total γ eq i are computed by summing their rates. In the non-relativistic limit the Boltzmann equation becomes where we defined such that DM is in thermal equilibrium at T ∼ M if λ i 1. In the annihilation case, λ ann is the usual dimension-less combination that helps obtaining approximate solutions in the limit λ ann 1. The thermal scattering term contains an extra unusual but important factor B(z), since it contains a Boltzmann suppression. Its non-relativistic approximation is  Figure 1: Evolution of the DM abundance Y = n/s for thermal decays, assuming a heavier initial-state involving a SM particle with mass m B = δ M and DM mass M = 1 GeV. The case δ = 1 corresponds to DM annihilations. For δ = 0 we considered scatterings with heavier initial state, interaction strength λ dec = 10 6 , degrees of freedom d DM = 2 and d B = 4. For δ = 0 (red curve) we instead assumed d B = 2 and a smaller λ dec = 10 4 : despite this it produces a lower DM relic abundance.
The ∆ factor that controls the exponential suppression is a combination of the masses M of DM, and m A and m B of the particles that become non-relativistic. The important observation is that, in general, ∆ = M . The DM mass M is the only relevant mass that controls how Boltzmann-suppressed the usual DM annihilations into light SM particles are. One has m A = M and thereby ∆ = m B for processes where the initial state is heavier and hence non-relativistic. As we are considering DM A SM B ↔ SM SM processes, ∆ is the mass of the initial-state SM particle. If instead the final-state is heavier and hence non-relativistic, DM SM ↔ SM A SM B , ∆ is given by a more complicated combination of masses.
In the presence of multiple process, the one with smallest ∆, i.e. least Boltzmann suppression, tends to dominate. The explicit expressions for ∆ imply that the Boltzmann suppression is only avoided at the kinematical border where DM decay starts being kinematically allowed.
We see that the 'thermal decay' terms corresponds to an effective decay width Γ DMT = Bσ dec0 that gets Boltzmann suppressed at low T . Indeed the rate of a true decay with width Γ DMT would be γ dec = Γ DMT n eq 1 K 1 (z)/K 2 (z) z 1 n eq 1 Γ DMT .

Approximate solution to the Boltzmann equations: freezeout
The freeze-out abundance produced by DM annihilations or by thermal decays can be approximatively computed as follows. These processes reach thermal equilibrium if λ ann,dec 1. Long before freeze-out, i.e., at early z z f , we can expand eq. (6) to first order in small Y − Y eq , finding for λ i 1 where δ ≡ ∆/M so that ∆/T = zδ. We see that δ < 1 i.e. ∆ < M enhances the effect of thermal decays compared to thermal annihilations. In general, the values of λ dec /λ ann and of δ decide whether freeze-out is determined by annihilations or decays. The final abundance is roughly approximated as Y ∞ ∼ Y eq (z f ) where the freeze-out temperature z f is estimated imposing Y − Y eq = cY eq with c ≈ 1. In the usual case where only DM annihilations are present, λ dec = 0, one gets the usual estimate for the freeze-out temperature z f ≈ ln 2λ ann /(d SM √ z f ) , and the usual freeze-out abundance • If instead only DM thermal decays are present, λ ann = 0, one gets where δ = ∆/M and W (x) is the 'product log' function, that solves x = W e W and is roughly approximated by ln x as x → ∞. The final DM abundance is approximated showing that thermal decays need a smaller value of λ 1 if δ < 1, i.e. ∆ < M . Fig. 1 shows examples of numerical solutions to the Boltzmann equations for fixed interaction strength λ dec and varying δ: we see that smaller δ leads to a lower DM final abundance in agreement with eq. (13).
• Thermal decays are maximally efficient for δ → 0. Let us discuss if/how δ = 0 can be realised. A look at kinematics shows that this limit corresponds to processes where DM scatters with a much lighter SM particle B, such as the photon. Thereby this case only arises near to the dangerous kinematical border where the phase space for DM decays starts opening up. Cancellation of IR divergences between real and virtual effects starts being relevant. In this special case B = Y eq B is constant (so that 2 This expression is accurate up to order unity factors, as we neglected evolution at z > z f . A more precise expression is obtained matching at the appropriate z f the approximation of eq. (11) valid at z z f with the approximation at z . However the matching does not result in a simple analytic expression. the non-relativistic limit used in previous approximation no longer holds); decoupling happens at z f ≈ λ dec Y eq B (no longer log-suppressed); subsequent evolution is non-negligible leading to This case is also illustrated in fig. 1.
Let us next discuss what these results mean in practice. We can approximate a nonrelativistic cross section at T ∼ M as where the interaction of DM with the SM sector is described by some dimension-4 interaction with dimension-less coupling g DM , or by some dimension-6 operator with coefficient G DM with mass dimension −2. Depending on how many particles are involved in this DM interaction, extra couplings can be needed to obtain a cross-section, so that we add extra powers p of a generic SM coupling (such as e or g 2 or f π /m π ). Since its power is model-dependent, for the purpose of the present parametrisation we assume g SM ∼ 1.
The condition λ > ∼ λ min translates into (thermal equilibrium) (16) where G F is the Fermi constant. Both annihilations and thermal decays need at least to reach thermal equilibrium, corresponding to λ min ∼ 1. DM annihilations give the freeze-out abundance Y ∞ ∼ 1/λ ann , so the cosmological DM abundance is reproduced for λ ∼ M/ eV 1, corresponding to For δ < 1, DM thermal decays need DM couplings smaller than in eq. (17), possibly down to the minimal value in eq. (16) for λ min ∼ 1, reached in the limit δ 1 i.e. ∆ M where DM thermal decays become maximally efficient. Fig. 2 shows the values of g DM and G DM needed to reproduce the DM relic abundance via freeze-out of thermal decays, assuming g SM = 1, considering different values of 0 ≤ δ ≤ 1 (solid curves). The needed values range between those characteristic freeze-out of DM annihilations (δ = 1, eq. (17) and black curves in fig. 2) down to those needed to achieve thermal equilibrium (δ = 0, eq. (16) and red curves in fig. 2). For δ < 1 the maximal DM mass M allowed by unitarity is above 100 TeV (the upper limit in case of DM annihilations): but model-building considerations will point to much lighter sub-GeV DM with very small couplings.

Approximate solution to the Boltzmann equations: freezein
If DM couplings are so small that it never thermalized, one can instead assume a vanishing initial DM abundance (needless to say, this is a more questionable assumption than thermal equilibrium) and match the small amount of generated DM with the observed DM abundance. The final DM abundance produced by freeze-in is obtained by setting Y = 0 in the right-handed side of eq. (6). One gets The DM abundance suggested by cosmology, Y ∞ ∼ eV/M , is thereby reproduced for a comparably small value of Γ DMT /H. Unlike in the freeze-out case, replacing DM annihilations with DM thermal decays does not lead to significant changes in the dynamics. What changes is that the needed scattering can be partially mediated by SM couplings, that can have significant values g SM ∼ 1. Freeze-in can be dominated at low T ∼ M or at high T M depending on the dimension of the interaction: • If DM has a dimensionless coupling g DM , the rate γ eq dec ∼ g 2 This matches the cosmological DM abundance for (freeze-in of thermal decays).  Precise order one factors depend on the process (such as a scattering or a decay, see [21] for examples). A general expression, valid up to these process-dependent order one factors, is obtained by inserting in eq. (18) the universal non-relativistic limit This is plotted in the left panel of fig. 2 assuming the cross section in eq. (15) with g SM ∼ 1.
• If DM has a non-renormalizable interaction G DM one has γ eq dec ∼ g 2 SM G 2 DM T 8 so that DM production is dominated at T M : the result depends on the full model. For example, the non-renormalizable operator could be produced by a mediator with mass M med and renormalizable couplings. If the reheating temperature satisfies T RH M med the DM abundance is dominated at T ∼ M med and can be approximated in a way similar to eq. (18) Finally, electro-weak symmetry breaking can dominantly contribute to the masses of some particles. While this phenomenon can enhance freeze-in above the weak scale (for example, making the relevant particles all nearly massless), it cannot provide a qualitatively new way of getting long-lived DM with desired cosmological abundance.

DM with lepton number and freeze-in
We seek theories where DM number is significantly violated, but where DM is very long lived. We consider DM coupled to two (or more) SM particles with masses m 1 and m 2 : To avoid the most immediate danger, the tree-level decay DM → SM 1 SM 2 , we assume M < m 1 + m 2 so that this process is kinematically blocked. We also need to assume M ∼ m 1,2 , so that the SM particles take active part in DM decoupling. Finally, if the coupling (21) is large, it might be necessary to assume M > |m 1 −m 2 |, to avoid potentially problematic decays of SM particles into DM.
The next danger is virtual decays such as DM → SM 1 SM * 2 . The SM 2 particle with mass m 2 and d 2 degrees of freedom is off-shell: the virtual SM * 2 has quadri-momentum k 2 with x ≡ k 2 2 /m 2 2 < 1 and decays into lighter SM particles SM 3,4,... not directly coupled to DM. The DM decay width averaged over polarizations is The integral is dominated by the highest available x max < 1, so that some decay modes of SM 2 might be closed for SM * 2 . The off-shell decay width is thereby approximated as where Γ i is the total decay width of SM i , and BR i (possibly equal to 1) its branching ratios into light enough particles (such as photons, neutrinos, electrons). Each SM i particle has a value of i , and its smallness is a figure-of-merit for coupling long-lived DM to SM particles. Table 1 lists the most narrow SM particles. All of them are lighter than about a GeV (heavier weak-scale SM particles decay too fast for our purposes), so DM too has to be similarly light. At sub-GeV energy, physics is described by baryons and mesons, rather than by quarks. In order to possibly reproduce the DM abundance via freeze-out, DM interactions cannot be suppressed by a scale much above the weak scale. The electron and proton are stable in view of conservation of charge and baryon numbers. The other SM particles undergo weak decays suppressed by (m/v) 4 and possibly by phase space (the neutron) or by flavour (the K, Ω, Σ in table 1). DM decays can be suppressed by ∼ 10 −15 , far from the needed < 10 −30 . One dose of off-shell vaccine does not protect DM stability well enough. Two doses of the off-shell vaccine can give enough protection. Let us assume that DM is coupled to two SM particles as in eq. (21) and its mass is M < m 1 , m 2 . Then doubly-virtual tree-level decays DM → SM * 1 SM * 2 are suppressed by 1 2 ∼ 10 −30 , corresponding to two weak interactions. An example in this sense is DM coupled to µπ, see fig. 3a.
However, even if DM only directly couples to massive SM particles, so that its tree-level decay channels are suppressed as discussed in the previous section, one-loop effects can give DM decays into light SM particles. DM that couples to a particle/anti-particle pair (such as π + π − ) decays at one loop level into photons, DM → γγ, with rate suppressed by a QED loop factor, ∼ (e/4π) 4 ∼ 10 −5 . This is far away from the needed 10 −30 suppression. DM that couples to π ± ∓ decays as DM → νγ at one loop level (see fig. 3b), with rate suppressed by an electroweak loop factor, ∼ (f π /4πv) 4 ∼ 10 −18 , but still below the needed 10 −30 . Loop decays are present because this DM has lepton number, which is carried by nearly-massless neutrinos. In section 5 we discuss non-minimal models that achieve the level of stability needed for freeze-out via thermal decays. We here discuss a minimal model that achieves the level of stability needed for freeze-in.

Mediator models
We consider DM to be a fermion χ with no SM gauge interactions and lepton number one. This corresponds to a sterile neutrino in presence of a Yukawa coupling χLH. This is a well-known allowed DM candidate if it has keV-scale mass: it can be cosmologically produced, compatibly with bounds on χ → νγ decays [22], as freeze-in, via resonant oscillations or via the decay of an extra scalar.
We instead explore the possibility that the SM neutral lepton singlet χ has a heavier mass M and has gauge interactions in an extension of the SM. The natural possibilities allowed by group theory are the following: • A (χγ µ R ) W µ R interaction in models where the right-handed W R boson is part of an extra SU(2) R . The full couplings of the W R mediator to fermions are where i, j are flavour indices. Integrating out the W R generates [9,5], that can arise in models with a SU(4) Pati-Salam symmetry [23,24,11,12]. Its couplings to fermions are We will discuss in section 4.1 how a flavour structure can arise. The latter term has been considered as a possible interpretation of the flavour anomalies in decays of B mesons, adding interest to the possibility of obtaining DM in this context. Integrating out the leptoquark at tree level generates [25]. The second operator in eq. (27) can also be generated by the exchange of a vector leptoquark A vector leptoquarkŪ 1 ∼ (3, 1) −1/3 can couple to χ and d R , but not to other SM currents, and thus cannot generate any DM-SM 3 interaction.
• Alternatively, interactions can be mediated by a variety of scalar fields. An example is the already mentioned Yukawa coupling with a Higgs doublet, H cL L χ, in which case the DM can be identified with a right-handed neutrino [22]. Scalar leptoquarks with suitable quantum numbers, instead, can couple χ to SM quarks [5,10,26,27]. The only two possibilities that are able to induce an interaction with SM leptons In eq. (29) and (28) above we omitted di-quark couplings for S 1 andṼ 2 that would violate baryon number. A scalar leptoquarkS 1 ∼ (3, 1) 2/3 can also couple to χ and have di-quark couplings, thus inducing a mixing between χ and the neutron.
The full list of possible DM couplings to one lepton and two quarks, all of which can be generated by the tree-level exchange of one of the previous mediators, is [28] L eff We expect that the various operator coefficients satisfy C i < ∼ (TeV) −2 because the charged or colored mediators that generate the effective interactions are generically constrained to be heavier than the TeV scale by collider bounds. We will sometimes collectively denote these couplings as 4G DM / √ 2, including a normalization factor for ease of comparison with the Fermi constant G F .
Below the QCD scale L eff DM becomes effective couplings to pions, χ π, and to heavier mesons (such as ρ and excited pions), as well as couplings to baryons, χ pn. Heavier mesons give negligible effects in cosmology and are not subject to significant bounds, despite that kinematical space for their tree-level decays is open. The coupling to baryons can have more significant effects, since baryons remain as relics at T m p,n thanks to the baryon asymmetry (DM couplings to quarks never lead to DM coupled to baryons but not to pions). Using chiral perturbation theory (see appendix A) the coupling to pions is obtained substituting the quark terms in the effective operator eq. (31) with pion terms with the same transformation properties: where f π = 93 MeV, B 0 ≈ 22f π is the quark condensate, D is the gauge-covariant derivative and we will not need higher orders. We assumed DM couplings to right-handed quarks, but this information is lost in DM coupling to pions (the suppression is of order f π /Λ QCD ∼ 1). This will remove a suppression of weak loops.

DM decays
We start discussing loop decays into massless particles, which are always kinematically allowed, and then move to tree-level decays.
The electroweak interactions couple pions to the charged lepton current (the SM Lagrangian is written in eq. (59)). As a result leptonic DM with the interactions of eq. (31) acquires at one loop a negligible mass mixing between DM and neutrinos 3 , and a mixed magnetic moment with neutrinos: 4 A precise estimate is not possible, as the exact loop integral is dominated by momenta around the QCD scale where all resonances contribute. The resulting DM decay rates are 3 In the pion effective theory the mass mixing µ χν + h.c. is UV-divergent and estimated as The UV cut-off has been estimated as Λ 2 UV ∼ M 2 W /(4π) 2 . Indeed, in the high-energy theory above the QCD and EW scales, the analogous effect is a log RG mixing between two dimension 6 operators: eq. (31) at two loop induces χLHH † H, so µ ∼ C χR v 3 y y u y d g 2 2 /(4π) 4 . The extra HH † and thereby y u y d terms appear because, in the full theory, DM is coupled to right-handed quarks. 4 Above the QCD and EW scales, the analogous effect is negligible: the dimension-6 operator of eq. (31) induces at 2 loops the operator (χσ µν LH)B µν HH † . DM  π * µ DM  π * µ + πµ * DM  π * µ * + π 0* ν DM  πµ π  µ DM DM mass M Figure 4: DM with lepton number coupled to µπ: kinematical thresholds for most dangerous decays. A * denotes a virtual particle, and thereby a Γ/M suppression. We here neglected the mass difference between charged and neutral pions.
Tree decay DM  ν ℓ π 0 Loop decay DM  ν ℓ γ 10 -2 10 -1 1 1 10 1 10 2 10 3 10 -10 The dashed curve is the contribution from DM loop decays into ν γ; the continuous curve is the contribution from tree-level decays into ν π 0 → ν γγ. Right: DM coupled to right-handed leptons with coupling C χR = 4G F / √ 2. Curves can be rescaled for different couplings taking into account that life-times scale as 1/C 2 χL,χR , while the bounds from cosmological DM stability and from CMB do not depend on DM couplings, and are discussed in the text.
Coming to tree-level decays, they can a) dominate over loop decays if DM → π is kinematically open, or b) be comparable to loop decays if one of the two SM particles needs to be off-shell to kinematically open the decay, or c) be negligible compared to loop decays if both * π * need to be off-shell. Let us discuss in more detail the case of DM coupled to = µ as a concrete example. Fig. 4 shows the main kinematical regimes. In the co-stability range 34 MeV = m π − m µ < M < m π + m µ = 245 MeV (36) all tree-level decays among χ, π, µ are closed. DM heavier than m π + m µ decays at tree level as DM → µ ± π ∓ . One of these particles must be off shell if the DM mass M is in the range m µ < M < m π + m µ , suppressing the DM tree-level decay rate by a µ , π ∼ 10 −17 factor. Furthermore, π → µ * χ → eν e ν µ χ decays are suppressed by µ ∼ 10 −16 , and µ → π * χ = eν e χ decays by π ∼ 10 −19 . So bounds BR(π → µχ) < ∼ 10 −7 are easily satisfied. If m π − m µ < M < m µ (in green in fig. 4), the tree-level DM → µ * π * decay rate is suppressed by π µ ∼ 10 −33 as both µ and π must be off-shell (see appendix A.1 for details). However, in this case, loop decays dominate. Fig. 5 shows the life-time of DM coupled to left-handed (left panel) and right-handed (right panel) leptons and quarks. We fix a reference value for C χR and C χL equal to The solid curve in the left panel is the contribution from tree-level decays DM → ν π 0 . For M < m π 0 this channel is kinematically closed and therefore the relevant contribution is into ν γγ via the off-shell π 0 . The dashed curve is the contribution from DM loop decays DM → ν γ which is independent on the lepton mass, as shown in eq. (35), and is the dominant decay channel for M < 2 MeV.
In the right-panel of the same figure, the solid curves indicate the contributions from tree-level decays, DM → π. These predictions depend on the leptons (blue, red and green are for e, µ and τ respectively) and the different features are due to the kinematical thresholds as discussed in details above. The dashed curves show the contribution from DM loop decays into ν γ, which depend on the lepton mass since we need a mass insertion to flip the chirality of the lepton (see eq. 35). For a given lepton the loop decay is dominant at M < m + m π and therefore the co-stability regime is not important from the phenomenological point of view.
For life-times longer then the age of the Universe, limits from Planck and Voyager II apply. These bounds are shown as grey areas in the figure. For DM coupled to lefthanded leptons, since there is always a neutrino that does not release energy into the plasma, we rescale the Planck bound for the channel DM → γγ given in fig. 7 of [29]. In particular for M > 2 MeV we rescale the limit by a factor 2/3 (the final state is ν γγ for both on-and off-shell π 0 ), while for lighter DM by a factor 1/2 since the loop decay into ν γ dominates. The bound is more complicated for DM coupled to right-handed leptons. For M < m + m π the dominant channel is into ν γ, so we follow the same procedure described above. For M > m + m π the final state after the decay of the pion and lepton always involves a e − e + pair with a given number of neutrinos (3, 5 and 7 for e, µ and τ respectively). We use the Planck and Voyager II limits given in [19,20] for the decay channel DM → e + e − rescaling them by a factor 2/5, 2/7 and 2/9 for e, µ and τ respectively. This is a good approximation for Planck, since the deposited energy matters, while for Voyager II this is a rough approximation given that the energy spectrum of the e ± changes dramatically as we have a chain of multiple decays.

DM thermal decays
We can neglect scatterings that involve two DM particles, as their cross sections σ 2 arise at second order in the small new physics couplings C of eq. (32). On the other hand, 2 ↔ 2 scatterings that change DM number by one unit arise at first order, by combining the C χ π interaction with either process Boltzmann factor ∆ Figure 6: Assuming DM coupled to π, we list the main processes that change by one unit the DM number. Their contribution to cosmological DM freeze-out is controlled by the exponent δ = ∆/M , plotted in the figure. Standard DM annihilations correspond to exponent δ = 1, and more efficient scatterings correspond to δ < 1.
• an electro-magnetic interaction with coupling e ∼ 0.3, adding an extra γ to the process (rates in Appendix A.2); or • a strong interaction with coupling ∼ m π /f π ∼ 1 (see eq. (32); the SM pion Lagrangian contains no π 3 interactions), adding an extra massive π to the process. The rates are given in Appendix A.3.
So these processes have a larger cross section σ 1 σ 2 , that needs a much smaller critical value σ 1 > σ cr 1 to affect DM cosmological freeze-out, if their kinematical factors ∆ that control the Boltzmann suppression are smaller than the DM mass M . The most important processes are shown in fig. 6. It shows that some off-shell decays can be less Boltzmann suppressed than DM annihilations, but only when M > m , so that tree-level DM decays are not suppressed by two off-shell factors. 5 Fig. 7 shows the final result. Its left panel is reported the thermal decay rates as function of z = M/T : as expected the process with smallest δ (χγ ↔ πµ in the example) dominates at low temperature. The green curves in the right panel show the coupling needed to reproduce the cosmological DM abundance as freeze-out via thermal decays (one curve is the numerical result, one nearby curve is the analytic approximation). It shows that DM interactions to µπ needs to be a few order of magnitude stronger than weak interactions. As expected, freeze-out is excluded by too fast DM decays [19,20] (red curve).   Figure 7: Dark Matter as a fermion coupled to µ ± π ∓ . Left: rates of main thermal decays compared to the Hubble rate. As expected, χγ ↔ πµ dominates at T M . Right: parameter space (mass, coupling). The red region is excluded by DM decays at tree (dashed) or loop (continuous) level, estimated as in eq. (35). The blue region is excluded by collider data constraining the quark-level effective operators of eq. (31) as explained in section 3.4. The DM abundance is reproduced thermally along the upper green curve, or as freeze-in below the lower green curve.

Collider bounds
Let us discuss collider bounds on the effective (χγ µ q)(qγ µ ) interaction. We find that the strongest constraints are set by the ATLAS analysis in [30], which perform a search for a heavy charged boson in events with a charged lepton and missing transverse momentum. In order to derive the bound, we have performed a recast of the ATLAS analysis. In particular, we have considered the reported bounds from the model-independent analysis in the muon channel, based on cuts on the variable m T , defined as m T = [2p T E miss T (1 − cos φ( , E miss T ))] 1/2 . In order to identify the fiducial region, we selected: We considered the model where a leptoquark, with mass M med √ s, mediates the effective interaction in eq. (31), obtaining the constraint: These bounds surpass the limits from HERA [31].

Calculations have been made with
MadGraph [32]. Note that the limits are obtained in the effective field theory approximation and we thus expect O(1) corrections in the lower mass region where M med ≈ √t ≈ TeV, with √t denoting the energy exchanged in the t-channel propagator of the leptoquark. This result also excludes thermal freeze-out, and is plotted as a blue curve in fig. 7. This conclusion still holds taking into account DM scatterings with baryons (such as χp → n + ), that do not get Boltzmann suppressed at low T thanks to the baryon asymmetry Y B ≈ 0.8 · 10 −10 . Similarly to weak interactions of neutrinos, these processes maintain DM in thermal equilibrium down to T ∼ (M Pl Y B G 2 DM ) −1/3 . However, if C < ∼ G F , this decoupling temperature cannot be below the desired DM mass M .

DM relic abundance via freeze-in
While freeze-out is excluded, freeze-in is allowed. Three different regimes are possible, depending on the highest temperature T RH attained by the SM plasma: 1. At T < ∼ Λ QCD in the confined phase, DM interacts with pions according to eq. (32).
The dimension-5 interaction with right-handed leptons contributes to the DM freezein abundance as Y ∼ e 2 C 2 χR f 2 π M Pl T : the DM abundance is dominated by the maximal temperature T RH . The minimal amount of freeze-in (upper boundary of the green region in fig. 7b) arises if T RH ≈ M . This confined contribution alone allows to reproduce the DM abundance for C χR ∼ T 0 /M M Pl /ef 3/2 π ∼ 10 −6 G F , nearly compatible with DM stability for M < m π + m µ . For the renormalizable interaction with left-handed leptons, freeze-in happens around T ∼ M , and the abundance The pion description breaks down at T ∼ Λ QCD when the QCD phase transition happens.
2. At temperatures above the QCD phase transition and below the mass M med of the particle that mediates the dimension-6 effective interaction with quarks in eq. (31), it leads to a higher Y ∼ C 2 χR M Pl T 3 . This is again dominated by the highest temperature T < ∼ min(T RH , M med ), where M med < ∼ C −1/2 χR . Since very small DM couplings are needed in order to get the cosmological DM abundance, we can neglect annihilation processes that involve two DM particles, described by extra effective operators of the type |ū R γ µ χ| 2 .
3. At temperatures above the mediator mass, DM has renormalizable interactions and freeze-in is dominated at T ∼ M med . The dimension-less DM coupling to the mediator must be small, in order not to over-produce DM. Other interactions (such as gauge interactions with SM particles) can be large and keep the mediator in thermal equilibrium, so that the freeze-in DM abundance is produced by its decays. In Boltzmann approximation, the space-time density of mediator decays into DM is where Γ med is the mediator partial decay width into DM, and d med is its number of degrees of freedom. 6 Inserting γ eq in eq. (18) gives If the mediator is the U leptoquark one has Γ med = g 2 χ M med /24π and d med = 18, and the DM abundance is matched for We remind that the coefficients C of the effective Lagrangian depend also on other couplings of the mediator. We now explore the possible flavour structure of these couplings, to see if some natural explanation can be found for the small size of g χ .

Signals and anomalies 4.1 Flavour expectations for DM couplings and the B anomalies
The fact that a leptoquark such as U µ could mediate at the same time the semi-leptonic interactions responsible for the B-physics anomalies, and the interactions responsible for DM freeze-in, although with couplings of different sizes, makes it interesting to speculate about a common origin of the two processes. The B-physics anomalies can be explained, compatibly with other bounds from flavour and collider physics, in models where new physics couples dominantly to the third generation of fermions [24,5,[35][36][37][38]. Couplings to lighter generations of quarks and leptons then arise from the flavour rotations that diagonalize the Yukawa interactions, and are suppressed similarly to the masses and mixings of light fermions. This can applied to vector leptoquarks assuming that only 3rd-generation fermions are charged under an extended Pati-Salam group [39][40][41]. In the left-handed quark sector these rotations are of the order of the CKM matrix.
Given the tiny values of the coupling g χ needed to reproduce the cosmological DM abundance via freeze-in, for DM masses in the MeV range, it is then natural to consider the DM χ as a right-handed neutral lepton of the first family. Other generations of χ, with larger couplings to SM, can be so heavy that they decay back fast to SM particles, without contributing to DM. 6 We can neglect scattering rates because suppressed by extra powers of couplings, such as g 2 . At T M med the decay rate is suppressed by a Lorentz factor T /M med compared to the scattering rates γ eq ∼ g 2 2 g 2 χ T 4 . This enhancement of scattering rates is not important for the final DM abundance, because it is dominated at T ∼ M med . Furthermore this enhancement can be partially included as a NLO contribution to the decay rate, describing it as the thermal contribution to the mediator mass, δM 2 med ∼ g 2 2 T 2 . The thermal mass avoid an apparent IR-enhancement at T M med of t-channel scattering rates, that would anyhow cancel with loop corrections in a full NLO computation (KNL theorem, see [33,34] for a similar computation).
We parametrize the SM Yukawa couplings as where f i are small parameters carrying the flavour suppression of the couplings to the i-th family of f = {Q, u, d, L, e}, and the relations hold up to O(1) factors. Here and in subsequent estimates we also ignore an overall coupling, assumed to be of order one. Such a pattern of Yukawa couplings is easily obtained in Froggatt-Nielsen scenarios [42], or in models with partial composite fermions [43][44][45]. More stringent relations between the i f can be obtained in models with a larger flavour symmetry, such as Minimal Flavour Violation [46][47][48] or U(2) models [49,50].
The values of Q,u,d can be estimated from the CKM mixings and the quark masses. Indeed, from which one gets where λ ∼ 0.2 is the Cabibbo angle, and η q is an overall O(1) parameter. The parameters in the lepton sector can not be estimated as easily at this level, since the neutrino mixings are not necessarily directly related to the Yukawa couplings.
We can now make the assumption that the flavour structure of the new physics couplings is determined by the same flavour parameters that control the SM Yukawas. Let us focus on the case of the vector leptoquark U ∼ (3, 1) 2/3 for definiteness. The couplings of U to fermions are defined in eq. (26). In terms of the flavour spurions those couplings read where we have introduced an additional (small) parameter χ that controls the DM couplings. The exchange of a vector leptoquark contributes to semi-leptonic left-handed interactions via the couplings g L . These couplings can thus be partially determined by fitting the b → s anomalies [5]. The contribution to the Wilson coefficient of the semi-leptonic operator (b L γ µ s L )(μ L γ µ µ L ) reads which fits the observed deviations if ( L 2 ) 2 η 2 Q /M 2 U ≈ (6 TeV) −2 . For M U ≈ TeV the muon couplings can be L 2 ∼ O(0.1). The first-generation coupling is required to be small, L 1 L 2 , in order to suppress new physics effects in b → se + e − . The right-handed couplings e are then related to the charged lepton Yukawa couplings, L i e j ∼ (y e , y µ , y τ ), and one has where we assumed L 1 O(10 −2 ) (an order of magnitude smaller than the muon coupling) and L 3 ∼ 1.
Notice that all the right-handed mixings u,d,e estimated above are smaller than the left-handed ones Q,L (except for the case of the electron where no hint on the size of L 1 exists). In particular, this is consistent with the observation of sizable LQ-induced effects in left-handed currents, but the absence of large right-handed currents.
Finally, reproducing the cosmological DM abundance via freeze-in as in eq. (41), fixes the LQ coupling to DM, and in turn χ ∼ g χ / u . Two limiting scenarios can be envisaged, depending on which up-quark flavour couples dominantly to DM. If the dominant decay is with the top quark, given the hierarchy of (44), for M 10 keV and M U ∼ TeV one needs χ 10 −8 , much smaller than the corresponding electron spurion. If instead χ couples dominantly to the first family of quarks, due to the additional u 1 suppression one gets χ ∼ 10 −5 , which is comparable to the analogous factor for electrons. The latter scenario could be motivated if χ is identified with the first-family member of a flavour triplet of right-handed neutrinos χ i : in this case its dominant coupling can be to other first-family fermions if flavour violation is suppressed in the right-handed sector, as is the case for instance in MFV or U (2) models.
We finally discuss why the Higgs doublet H, unlike the leptoquark U, cannot generate DM via freeze-in. The Higgs doublet too can couple to DM as where we estimate the Yukawa coupling as y i χ ≈ χ L,i . Thus Higgs decays h → ν L χ can in principle contribute to DM freeze-in production. However DM dominantly produced in this way is excluded, because eq. (48) implies a too large mixing of DM with neutrinos, θ i ≈ y i χ v/M χ . This is problematic because it gives a contribution to neutrino masses and, more importantly, because it contributes to DM decays as Higgs freeze-in instead needs a Yukawa coupling y χ ≈ 7 × 10 −8 keV/M , which is slightly larger for masses M few eV. And ad-hoc scalar coupled to two DM particles was introduced in [51] to achieve freeze-in: we have shown that this role can be played by the vector leptoquark motivated by flavour anomalies, as its Uuχ couplings do not directly induce DM decay.
In any case, we need to assume that y χ is small enough so that eq. (49) is satisfied and the Higgs plays no role in freeze-in. In this case, the corrections to neutrino masses do not exceed the experimental values and do not impose further constraints. We point out that the required Yukawa coupling is smaller than the naive expectation y χ ≈ χ L even in the case where flavour violation is suppressed in the lepton sector and the dominant coupling is to electrons. The need for small mixings between left-and right-handed neutrinos is a common issue in models where the leptoquarks arise from partial unification at the TeV scale, and most likely requires a more complex neutrino sector [40,41].

Direct detection signals
Direct detection experiments are sensitive to DM scattering with light degrees of freedom such as light quarks, gluons and electrons. We consider DM that interacts with SM leptons and quarks as described by eq. (31), and in particular on the first two 4-fermion interactions. After matching the interaction Lagrangian to the nucleon level, DM chargedand neutral-current interactions with nucleons arise. We now briefly discuss the novel signatures in direct and indirect detection searches of DM induced by these interactions. We anticipate that, although these signatures are potentially interesting, the effective scale of the dimension 6 operators in eq. (31) must be large to avoid bounds from DM decay as shown in fig. 5. Therefore, one generally expects the signals to be much below the experimental sensitivities -with a few notable exceptions that we discuss.
Concerning charged currents, the leading order amplitudes are obtained from the (axial) vector interaction with right-handed leptons C χR χγ µ e R nγ µ 1 + g A γ 5 2 p , where g A ≈ 1.27 [52], and the scalar interaction with left-handed leptons (χe L ) (np), with a coefficient proportional to C χL . Both these interactions lead to a DM-induced β decay. At the nuclear level, the relevant process is the β − transition DM + A Z N → e − + A Z+1 N , which is energetically possible if the DM mass is larger than the capture threshold Q+m e , where Q = m A,Z+1 − m A,Z is the mass difference of the two nuclei. Following [53] the thermally averaged cross section arising from the charged current operator is where f (E R ) is a function of the detected energy that encapsulates the nuclear details of the β − transition. For light DM particles with kinetic energies much less than δ = M − (m e + Q), the spectrum of recoiling electrons exhibits a peak at E R ≈ δ. The left panel of fig. 8 compares the estimated sensitivities of current direct detection experiments, taken from [53], with bounds from DM decays and from colliders.
One could wonder whether there is a possibility to explain the anomalous counting rate of electromagnetic recoils in the Xenon1T detector peaked around 3 keV [54]. Considering Xenon as target, the process DM + A 54 Xe → e − A 55 Cs can have capture thresholds as low as a few hundreds keV: in particular, the 131 54 Xe isotope has an abundance of about 21% and a capture threshold of only 355 keV [55]. One therefore can have a signal peaked at E R ≈ 3 keV if the DM mass is M ≈ (m e + Q) + E R ≈ 358 keV. This needs some tuning among the electron and DM masses, and a slightly larger M would give a peak in the few hundreds of keV range. We show the best fit to the Xenon1T data in the right panel of fig. 8, obtained using the rate reported in [53]. Assuming local DM density ρ = 0.4 GeV/ cm 3 , the anomaly can be fitted for values of the DM coupling C χR ≈ 1/(3.4 TeV) 2 . For this value of the coupling to right-handed electrons, the loopinduced decay rate into ν e γ gives a DM life-time longer than the experimental bound from Planck, as shown in the left panel of fig. 8. However, the LHC bounds on the lepton plus missing energy final state, obtained in the EFT approximation in eq. (38), constrain C χR to be a factor 2 smaller. 7 A similar coupling to left-handed electrons is instead largely excluded by DM decays. The Xenon signal can thus most likely not be explained in our framework, but future experiments will have the sensitivity to probe a parameter space relevant for DM freeze-in via coupling to right-handed leptons, provided that the reheating temperature is low, T RH ∼ Λ QCD . Furthermore a novel specific signature from the 136 55 Cs decay arises. Indeed, following the analysis done in [55] for solar neutrino capture in LXe, one could detect DM with lepton number by using a delayed coincidence signature from long-lived states of 136 55 Cs. Concerning neutral currents, DM couplings to left-handed leptons induce an interaction with nucleons N described at leading order by the scalar operators (χν L )(nn) and (χν L )(pp).
This is spin-independent and leads to a sort of exothermic DM-nucleus collision DM N → ν N with mass splitting δ = M − m ν M . For light DM particles with kinetic energies much less than the splitting, exothermic spectra exhibit a peak at nuclear recoil energy where m N is the mass of the target nucleus. Hence, for Xenon target we expect detectable nuclear recoil in the keV range when the DM mass is around 20 MeV. Nuclear recoil affects atomic electrons, inducing extra bremsstrahlung and Migdal emission. It is worth noticing that the window of DM masses probed in such scenario is much smaller than the one achieved with standard DM-nucleus collisions. However, the spectrum is analogous to the one of exothermic collisions and therefore we cannot discriminate further among them.
DM absorption signals have been studied previously in [53]. However, when discussing neutral currents these authors considered a (χν R )(qq) effective interaction involving ν R (rather than ν L or ) in order to avoid too fast DM decay χ → ν R γ. Since ν R has no SM gauge interactions their signals are equivalent to those of inelastic DM with two states χ, Xenon1T data Xe + χ  131 Cs + e -Background Figure 8: Left: Estimated experimental sensitivities for the DM absorption cross section (from [53]), compared to the indirect detection bound on DM life-time (these are stronger than in [53], where the dominant QCD contribution to DM → νγ was not considered) and to LHC bounds estimated in the EFT limit (red line). The magenta star shows the parameters that can fit the Xenon1T excess. Right: Best fit to the Xenon1T electron-recoil data from DM absorption on 131 Xe.
χ and a (χχ )(qq) interaction. When discussing charged currents [53] considers our same (χ )(du) interaction, but does not include the DM decay channel χ → νγ mediated by QCD states at one loop. We report in the left panel of fig. 8 the projected constraints for the DM-induced β − signals taken from [53]. One can see that, while DM coupled to right-handed electrons can give observable signals in various experiments, for DM coupled to left-handed leptons the signals compatible with DM life-time are always below the experimental sensitivity.

Indirect detection signals and the 3.5 keV anomaly
DM decays provide indirect detection signals: the constraints on tree-level and loopinduced DM decay from Planck and Voyager have already been discussed in section 3.2. As these constraints can be saturated, signals from DM decays are detectable.
In particular, the unidentified X-ray emission line claimed at 3.5 keV [57,58] could be due to DM with mass M ≈ 7 keV decaying as DM → ν γ with rate Γ ≈ 1.4 10 −28 / s (see e.g. fig. 13 of [57]). In the considered models the needed DM life-time is obtained for For such values of the couplings the relic DM abundance can be reproduced via the freeze-in mechanism with a low reheating temperature T RH ∼ Λ QCD , see fig. 7. This is a non-trivial result, as other models where DM is a sterile neutrino can fit the 3.5 keV line, but its relic abundance cannot be obtained via freeze-in as non-resonant oscillations: extra mechanisms are needed, such as decays S → DM DM of an ad-hoc scalar S [51] (possibly with a partially closed phase space, if DM needs to be colder [59]). In our model on the other hand the leptoquark motivated by flavour anomalies plays a similar role, decaying into one DM particle and one SM particle. Concerning scatterings that could induce a gamma-ray signal, the only relevant interactions are those that change the DM number by one. Indeed the pair-annihilation of DM particles is extremely suppressed within the freeze-in production mechanism, as σ ann v ∝ g 4 χ and the g χ coupling is very small. We thereby focus on further indirect detection gamma-rays signals coming from DM decay and from scatterings that change DM number by one unit. A possible way to get a signal is to change the DM number by one via the scatterings with highly energetic Galactic Cosmic-Ray (CR) electrons and protons. The relevant primary channels are χe + →du and χp →ne and the corresponding charge conjugates. The spatial morphology of the induced gamma-ray signal is non-conventional and is enhanced in regions with high density of baryons and DM. As a consequence dwarf galaxies, which are usually considered as the cleanest laboratories in indirect detection, are not the best place to look for gammarays signals from DM with lepton number. The energy spectrum is similar to the one from standard pair-annihilation with the difference that the total energy is not 4M 2 but (M + E CR ) 2 with E CR the energy of Galactic CR. The cross section of highly energetic Galactic CR electrons and protons scattering off non-relativistic DM is The highest C χL(R) ∼ 10 −4 G F allowed by DM stability and consistent with the freeze-in mechanism with low reheating temperature T RH ∼ Λ QCD is obtained for very light DM. Hence we predict σ 10 −9 pb (E e,p /GeV) 2 , which gives a flux too low to be detected by any indirect detection experiments.

Flavour models with multiple DM and freeze-out
So far we considered minimal models with one DM particle, and we could not reproduce its cosmological abundance from freeze-out via thermal decays, compatibly with bounds on the DM lifetime, τ > ∼ 10 26 s, although we found viable models assuming DM freeze-in. We here show how bounds on the DM lifetime can be easily satisfied in the presence of two or more DM states χ i . 8 This extension could be motivated by flavour: for example, in the models with extra vectors such as eq. (26) one would naturally expect 3 generations of DM singlets χ 1,2,3 .
To keep the discussion simple and general, we summarize our previous results about the decay rate of one DM particle as SM ∼ 10 −18 is a typical electroweak loop factor at the DM mass M < ∼ GeV. If C ∼ G F , as needed to have thermal freeze-out at M ∼ GeV, the two suppression do not give a stable enough DM (while freeze-in needs C G F , improving DM stability). Stable enough DM is obtained in the presence of two states, χ 1 and χ 2 , with mass M 1,2 < ∼ GeV (arising, for example, as m + m π ) and effective interactions of the type Notice that we have not imposed any Z 2 symmetry that keeps the dark sector separated from the SM sector. We next assume that C DM 1 is negligibly small, and that masses M 1 < M 2 are in the range such that χ 2 → π decays are kinematically closed, while χ 2 → χ 1 π 0 decays are kinematically open (we make this latter assumption just for simplicity; we could relax it and use χ 2 → χ 1 γγ decays). Then the decay widths are where DM 12 ∼ C 2 DM 12 M 4 1,2 . The smallness of C DM 1 , at a level that saturates the decay width Γ DM 1 above, could be motivated by flavour considerations (for instance a very precise flavour alignment in the right-handed lepton sector, if χ carry lepton flavour number). Quantum effects respect this requirement. Assuming now that C DM 2 ∼ C DM 12 ∼ G F or even larger depending on the masses, we have, at the same time: 1. Very long lived χ 1 DM, that satisfies bounds on the DM life-time; 2. Cosmologically fast decays of χ 2 , so that it does not pose cosmological problems; 3. Thermal freeze-out via χ 1 ↔ χ 2 ↔ SM scatterings at finite temperature.
Although not necessary, DM stability can be improved by an extra DM 23 factor adding a third generation χ 3 and repeating the above structure: each mediation step can give one extra suppression at zero temperature, while rates at finite temperature can reach thermal equilibrium, Γ ∼ H at T ∼ M , if C DM ∼ G F . This latter condition (needed to have thermal equilibrium down to the DM mass) can be borderline or problematic in specific models.

Conclusions
We explored the possibility of matching the cosmological DM density via particle-physics processes that change the total DM number (DM particles plus their anti-particles) by one unit, dubbed thermal decays. Differently from the usual models where two DM particles annihilate, no Z 2 symmetry is introduced to distinguish the dark sector from the SM sector, and to keep DM stable.
In section 2 we wrote the relevant Boltzmann equations and provided approximate analytic solutions, finding that thermal decays can be much more efficient than DM annihilations, and thereby can need much lower rates. In the most favourable situation, thermal decays just need to reach thermal equilibrium at temperatures around the DM mass. Fig. 2 shows the needed values of DM couplings: dimension-less couplings need to be larger than 10 −9 , and dimension-6 operators can be suppressed by a scale only mildly larger than the weak scale.
The price is that DM can decay, and the bounds on its life-time are 30 orders of magnitude stronger than the DM rates needed for thermal freeze-out. So the above possibility seems largely excluded. To avoid this conclusion we considered DM in costability ranges such that main DM decays at zero temperature are kinematically blocked, while scatterings at finite temperature are open. Kinematical blocking allows to improve DM stability by one electroweak loop factor -about 18 orders of magnitude if the DM mass is around the muon or pion masses. This is not enough to allow for DM thermal freeze-out, and it is enough to allow DM freeze-in.
In section 3 we implemented this mechanism in models where DM is a right-handed neutrino: a SM singlet with lepton number one. Such a state can interact with SM fermions in presence of an extended gauge group, or additional generic mediators. This includes, for example, the TeV-scale Pati-Salam leptoquark recently considered to fit the B-physics anomalies. We find that these models can successfully reproduce the DM abundance via freeze-in, with small DM couplings possibly compatible with a Froggat-Nielsen-like picture of flavour.
We discussed possible specific signatures of these models, depending on the DM mass. Concerning indirect detection, DM decays for example allow to fit the 3.5 keV anomaly compatibly with freeze-in DM production, if M ≈ 7 keV and the reheating temperature is low. Concerning direct detection, one has unusual processes where DM can convert into SM particles. An explanation of the Xenon1T excess peaked around 3 keV is compatible with bounds from DM decays and freeze-in production, but disfavoured by LHC bounds. A dedicated experimental search is therefore needed to precisely establish order one factors. The collider signals are those of heavy vectors, and flavour signals include their effects.
In section 5 we found that extra suppression of DM decays, as needed to reproduce the DM cosmological abundance via freeze-out of thermal decays, can be obtained in models with two or more new states. Their presence could be motivated, for example, by flavour considerations.
The vector-like DM interaction in eq. (32) is similarly obtained, replacing the left-handed current with (χγ µ R ). The scalar interaction is obtained converting (χL L )(Q L u R ) into

A.1 Tree-level DM decay
If M > m , the DM can decay to a lepton and an off-shell pion (subsequently decaying to eν) through the interaction in eq. (32). The amplitude reads where the second term is absent if χ is Dirac. The squared amplitude is where p π = p e + p ν is the pion 4-momentum and the (2) Off-shell decay rates have been computed both analytically and numerically with MadGraph [32], after having implemented the interaction in eq. (32) and the relevant SM pion interactions.