Asymmetric dark matter: residual annihilations and self-interactions

Dark matter (DM) coupled to light mediators has been invoked to resolve the putative discrepancies between collisionless cold DM and galactic structure observations. However, $\gamma$-ray searches and the CMB strongly constrain such scenarios. To ease the tension, we consider asymmetric DM. We show that, contrary to the common lore, detectable annihilations occur even for large asymmetries, and derive bounds from the CMB, $\gamma$-ray, neutrino and antiproton searches. We then identify the viable space for self-interacting DM. Direct detection does not exclude this scenario, but provides a way to test it.

: Constrained regions (shaded areas in blue, pink, purple and cyan) in the parameter space of the DM mass M p D and mediator mass M V , together with areas of sizable self-interactions (shaded green), for the symmetric DM case, i.e. r ∞ = 1, and including the effect of the dark electrons. The areas of sizable self-interactions are ruled out by the CMB and FERMI dSphs constraints. This panel is to be contrasted to the asymmetric cases, r ∞ < 1, illustrated in Enter a dark asymmetry. In this paper, we entertain a slightly more complex possibility. We introduce a particle-antiparticle asymmetry in the dark sector, in analogy to and partly motivated by the ordinary baryon asymmetry [46,47]. The density of asymmetric DM (ADM) is determined, at least partly, by the excess of DM particles over antiparticles, is the DM (anti)particle-to-entropy density ratio.
The presence of an asymmetry suppresses the indirect detection signals [48][49][50], simply because there are fewer antiparticles to annihilate with the DM particles. Hence, large ranges of the DM and mediator masses that are ruled out in the symmetric DM scenario, are viable in the presence of an asymmetry. However, an accurate computation reveals that even highly asymmetric dark matter with long-range interactions can produce significant indirect detection signals, due to the residual DM annihilations being enhanced by the Sommerfeld effect. This is a non-trivial result, first pointed out in [51], which deserves further scrutiny.
The purpose of this paper is therefore twofold: (i) to compute indirect detection constraints on ADM with long-range interactions; (ii) to identify the parameter space where all constraints are evaded, while DM self-interactions are sizable.
The rest of this paper is organized as follows. In Sec. 2, we specify the details of the model. In Sec. 3 we briefly discuss the processes which set the relic density. In Sec. 4 we compute the DM self interactions, before deriving the indirect detection constraints in Sec. 5. Direct detection is discussed in Sec. 6 (providing updated constraints to [52,53]). Important caveats regarding the formation of dark atomic bound states and the reannihilation effect are explained in Sec. 7. We then briefly comment on our findings and conclude.

The Model
We consider a dark QED-like model, in which DM consists of Dirac Fermions that couple to a dark gauge U (1) D force. If DM also carries a particle-antiparticle asymmetry, then gauge invariance implies that there must be at least two dark species with compensating asymmetries, such that the total gauge charge of the universe vanishes. While this is evident in the case of an unbroken U (1) D , it remains inevitable under reasonable assumptions if U (1) D is mildly broken. Importantly, this encompasses the parameter space in which the dark vector boson is sufficiently light to mediate sizable self-interactions that can affect the galactic structure [54] 3 .
We shall thus introduce two dark species, the dark protons p D , and the dark electrons e D , that couple with opposite charges to the dark vector mediator V . For our purposes, both p D and e D are elementary point-like particles with masses M p D and m e D , and M p D m e D . The Lagrangian is where D µ = ∂ µ ± ig D V µ is the covariant derivative for p D and e D . The field strength tensor is is the dark fine-structure constant. In the case we consider here, there is a conserved dark baryon number associated with p D and a conserved dark lepton number associated with e D . A linear combination of dark baryon and dark lepton number needs to be broken for generation of the asymmetry, the inclusion of the dark electrons means this can be achieved in a gauge invariant way, in analogy with higher dimension baryon violating operators added to the SM. For concreteness we assume here that the Stückelberg mechanism generates M V . The introduction of a dark Higgs is not expected to change our conclusions, beyond minor numerical differences. It is important to note though, that the assumption of asymmetric DM implies that the Higgs mechanism does not generate a Majorana mass for at least one of the dark species, here the dark protons. This only restricts the U (1) D charge of the dark Higgs. For more complete models, see e.g. [55][56][57][58].  [59]. Bound state formation becomes important at relatively large values of α D and M p D and when the temperature drops far enough for the inverse ionization process to become slower than the relevant decay rates.

Relic Density
The required α D to obtain the observed Ω D from thermal freeze-out in the presence of an asymmetry is computed as in [51]. We take into account the Sommerfeld enhancement of the annihilation processes, as well as the formation and decay ofp D p D bound states. New here is that the dark electrons provide an additional annihilation channel, shown in Fig. 2, and an additional decay channel for the spin-1 bound state. The annihilation and bound-state formation cross-sections are where S ann is the s-wave Sommerfeld enhancement factor, and S BSF arises from the appropriate convolution of the bound-state and scattering-state wavefunctions and includes the Sommerfeld effect. S ann and S BSF depend only on the ratio α D /v rel in the Coulomb approximation [60,61], which is satisfactory during freeze-out [44] (for a caveat, see [62] and Sec. 7.2). The decay widths of the spin-0 (↑↓) and spin-1 (↑↑) bound states are and are summarised in Fig. 3. Since Γ(↑↑→ē D e D ) Γ(↑↑→ 3V ), the formation ofp D p D bound states depletes the DM density more efficiently. The e D ,ē D degrees of freedom also affect the darkto-visible sector temperature ratio after the two sectors decouple. The combined effect of the dark electrons is to suppress the required α D . The DM mass is related to the proton mass m p by 4 where Y B is the baryon asymmetry, and r ∞ ≡ (Y − /Y + ) t→∞ is the ratio of DM antiparticles to particles today. Note that a smaller r ∞ requires a larger annihilation cross-section [48,51], thus larger α D . Following Ref. [51], we employ unitarity -which sets an upper limit on the partial-wave inelastic cross-sections -to estimate the maximum mass for which DM can annihilate down to the observed density. In the present model, the direct annihilations are s-wave, while bound-state formation is pwave. Using the corresponding leading order cross-sections (2) and (3), we find that the s-wave limit is saturated for lower α D . We employ this value to estimate the maximum DM mass, taking into account the depletion of DM via both the s-wave and p-wave inelastic processes. Since smaller r ∞ requires more efficient annihilation, the maximum M p D decreases with decreasing r ∞ [51]. We note though that any computation close to the unitarity limit using (2) and (3) is likely inaccurate, since at those values of α D higher order corrections should become significant, such that the cross sections remain below the unitarity bound for any value of α D .

Self Interactions
The cross section relevant for comparison to simulations of SIDM is the momentum transfer cross section, which should take into account both attractive, p D −p D , and repulsive, p D − p D orp D −p D , interactions. The attractive and repulsive cross sections are weighted with the appropriate densities, in order to obtain the effective transfer cross section, where σ att (σ rep ) is the attractive (repulsive) transfer cross section, n ± ∞ is the DM (anti)-particle density today and n sym ∞ is the DM density today for symmetric DM of the same mass. We use the analytic approximations for σ att and σ rep given in Ref. [63]. The areas of sizable self interactions at velocities relevant for dwarf galaxies, v rel = 10 −4 , are shown as green bands in Figs. 1 and 4. For larger asymmetries (smaller r ∞ ), the importance of the attractive interaction decreases, hence the parametric resonances begin to disappear. The bands also move up to slightly higher values of M V because α D increases with decreasing r ∞ . Note σ T will be suppressed at velocities relevant for cluster collisions, v rel ≈ 3 × 10 −3 , due to the light mediator. Hence the constraints from systems with typically higher velocities [64] are naturally evaded.
In our computations, we neglect scatterings involving e D orē D . The effect of the dark electrons depends of course on their mass, and we do not attempt a detailed exploration along this parameter in the present work. We only briefly comment on the potential implications. In part of the parameter space, the dark electrons may increase the energy and momentum exchange inside halos, thereby shifting the preferred regions (green bands) towards higher M p D and M V values, which are less constrained by indirect detection. In particular, since the dark electrons are lighter than the dark protons, the p D − e D collisions are expected to be more frequent, albeit transferring less energy on average per event, than the p D − p D scatterings. On the other hand, in large portions of the parameter space where scatterings involving dark electrons are expected to be significant, dark atoms may have formed in the early universe (cf. Sec. 7.1), thus screening the DM self-interactions today [54]. Very importantly, the effect of DM self-interactions in multicomponent DM models has not been simulated. While semi-analytical estimates are possible (see e.g. [54,65,66]), they involve larger-than-usual uncertainties. Our work, in combination with the strong constraints on symmetric SIDM models [44,45,67], strongly motivate the development of multicomponent DM simulations.
Finally we note that for the range of mediator masses we consider here, M V > MeV, dissipation is expected to be negligible. For studies of dissipation in various atomic DM scenarios, see e.g. [41,42,68].

Indirect Detection Constraints
Due to the suppressed population of DM antiparticles, the effective cross section for indirect detection is [48,50,51] σ inel includes both the directp D p D annihilation, and the formation ofp D p D bound states, whose crosssections are given by Eqs. (2) and (3). For the velocities relevant to the DM halos and the CMB, the factors S ann and S BSF depend strongly on the mediator mass, and have been computed in [69]. We use σ ID to set limits from PLANCK, AMS, FERMI, and ANTARES data, and present our results in fig. 4. Before discussing the various indirect probes further, we note that we determine the decay branching fraction of the mediator into SM final states, BR(V →f f ), according to the discussion in [44], which we follow closely here. Effect of the dark electrons. In our analysis, we neglect the annihilations of the relicē D e D into dark photons. As long as m e D M p D , the dark electrons annihilate more efficiently in the early universe, and have a smaller relic population available to annihilate at late times. In the symmetric limit, r ∞ 1, the number densities scale as n ∞ (e D )/n ∞ (p D ) ∼ m e D /M p D , while for r ∞ 1, the residual symmetric population of dark electrons is more suppressed than that of the dark protons. Indeed, r ∞ decreases exponentially with the annihilation cross-section [48,51], thus r ∞ (e D ) r ∞ (p D ). A small m e D may also imply that the fluxes fromē D e D annihilations lie outside the energy ranges probed by FERMI, AMS and ANTARES. This last point does not apply to the PLANCK constraint. However, we find that the CMB constraint on α 2 D from e D annihilations is weaker by at least a factor m e D /M p D than that of dark protons 5 .
It follows that theē D e D pairs produced in thep D p D annihilation and bound-state decay processes also do not contribute to our constraints. Their density is very small and the probability that they subsequently annihilate among themselves or with relicē D , e D is entirely negligible.
The above imply that the indirect constraints weaken due to the presence of the dark electrons in the theory, even for r ∞ = 1. This is because the dark electrons contribute to the depletion of DM in the early universe, thereby reducing the predicted α D , while they do not contribute to the indirect detection signals at late times.
CMB constraints. Dark matter annihilation can increase the ionization fraction of the plasma during the CMB epoch which can affect the measured anisotropies. In our model the effect is determined by the DM mass M p D , the relative velocity during the CMB, v rel 10 −8 [44] 6 , the dark-photon branching fractions into SM final states, BR(V →f f ), which depend on M V , and the final-statedependent efficiency factor, f eff [70,71]. We apply the constraint derived from the PLANCK data [7] to obtain the bounds shown in Figs. 1 and 4, where for σ ID we have used the Hulthen approximation regularised with the prescription of [72] to respect unitarity, see [44] for more details. For M p D 100 GeV and r ∞ -dependent values of M V (e.g. M V ≈ GeV for r ∞ = 1) the region excluded by CMB is not precisely rendered in our figures because the resonances become increasingly dense. In these regions, the density of the coloured points provides a rough indication of the density of the actual resonances, and thus of the area excluded by the CMB. Constraints on ADM with contact only interactions from the CMB have been previously produced in [73]. They apply only to larger r ∞ and a limited DM mass range.
AMS antiprotons. We use the method of [44]. Constraints arise for M V 2 GeV for obvious kinematic reasons. Two parameter regions are constrained as clearly seen in Fig. 1. At M p D ∼ 10−100 GeV, any additions to the well measured flux are tightly constrained, and at M p D 100 GeV, the resonances imply signals that would overwhelm the measured spectrum.
FERMI Dwarfs. We follow closely the analysis in [44], which is based on the FERMI observations of 15 dwarf galaxies, with gamma ray energies in the range 0.5 GeV to 0.5 TeV [74]. The statistical analysis treats the J-factors as nuisance parameters. Two separate regions are excluded (clearly seen in Fig. 1). In the Sommerfeld regime, the enhancement compensates in part for the suppression of the photon flux when V decays mostly into leptonic channels.
FERMI Galactic Halo. Constraints from FERMI observations of the Galactic Halo were presented in [44], based on the analysis of [75]. We have updated the analysis considering a more recent dataset (Pass 8) and including more statistics. The constrained region disappears, except for some small areas close to the resonances. This is both due to the lower α D predicted in the presence of the dark electrons, and because the bounds become slightly less severe with Pass 8 data. We have checked that allowing the overall background level to float does not change this conclusion. We do not show these small excluded regions in our figures.
ANTARES neutrinos. We derive constraints from ANTARES [76], which were not considered in [44], and are particularly stringent for M p D 1 TeV. The ANTARES bounds come from the non-observation, from the Milky Way halo, of excesses in muon neutrino fluxes. The existing constraints are cast for DM annihilating directly into SM pairs. Therefore we cannot directly use them, because in our model DM annihilates first into V that then decay into SM particles, resulting in a one-step cascade that softens and broadens the spectra. We overcome this limitation as follows.
First, we compute the final ν µ spectra from DM annihilations directly into all SM pairs constrained by ANTARES, using the PPPC4DMID [77] with the electroweak corrections switched on [78] (which was used also by the ANTARES collaboration in [76]). We also compute neutrino spectra from one-step cascades in the simple limit M V M p D (see e.g. [79] for more details), which holds in the region where the bounds will apply. Then, we compare the spectra obtained with the above two procedures, and notice that the high energy ν µ fluxes -which drive the ANTARES limit 7 -are the least changed for quark-antiquark final states, and that the cascade yields a similar amount of high-energy neutrinos. Therefore, we exclude regions of our model in which for at least one of the two final statesf f =bb plus lighter quarks, andf f =ν e,µ,τ ν e,µ,τ . For thebb channel, we sum over all lighter quarks, as the spectra are very similar. Following ANTARES we also take into account the neutrino oscillations. Our constraints are conservative in the sense that, in our model, one would have to sum over all final statesf f , instead of constraining one final state at a given time. This procedure is impracticable here because it would require to perform the ANTARES analysis using, as signal flux, the specific one coming from our model, which is not immediately comparable with the flux from the final states considered by ANTARES. On the other hand, our constraints would be relaxed by the choice of a more cored profile than the one used in [76]. While [76] shows the impact of different profile choices for the τ + τ − final state, it does not do the same for the other ones. Therefore we choose to consistently use the only profile for which the limits is expressed for all final states.
The excluded area is shown in Figs. 1 and 4. It disappears somewhere between r ∞ = 10 −1 and r ∞ = 10 −2 . Following closely the procedure of Ref. [44] for the dark photon decays, we find another excluded region for 200 MeV M V 5 GeV and M p D 3 TeV. However, we believe that this exclusion is, at least partly, an artifact of the way dark photon decays are modeled, in the M V region where hadron final states are important (i.e. exactly in the M V region just mentioned). Indeed, in Ref. [44], the dark photon was assumed to decay with 50% probability into muons, and 50% in neutrinos, to reflect the dominant final states from π ± etc. The spectra of µ and ν µ would, in reality, be softer, because of the additional steps (the π ± ) involved. Therefore, consistently with our understanding that the limits are driven by the high-energy part of the spectra, we conservatively choose to not show on our plots the ANTARES exclusion in that region.
Other limits from high energy cosmic-rays. Our previous inclusion of ANTARES limits demonstrates that high energy cosmic-ray experiments probe otherwise unexplored regions of parameter space of this model. In this domain of energies, we have not attempted to recast recent limits from ICECUBE [80], VERITAS [81] and HAWC [82,83] to secluded models. This is motivated by the following: i) it is difficult, based on [80][81][82][83], to derive which energies of the measured cosmic-rays dominate the exclusions (see footnote 7 for ANTARES), which was a necessary ingredient for our simple recast of the ANTARES limits; ii) the limits on standard DM annihilations, given in [80][81][82][83], are much weaker than those given by ANTARES.
This second argument however does not hold for recent HESS limits [84] from 254 hours of observation of the GC. Nonetheless, to recast them for secluded DM models we could not use the same prescription we used for ANTARES limits, because of reason i) above. To recast HESS limits, a possibility would be to directly analyze HESS data. To our knowledge, however, these data are not public. To overcome this lack of public data, Ref. [85] used public data of the 2011 HESS limits on DM annihilations [86], from 112 hours of observation of the GC, analysed them for secluded models, and projected them to 254 hours. The HESS reach obtained in this way is potentially very constraining. We do not include it in our study, however, because it is partly a projection. Moreover, unlike any other bounds we included, HESS loses any sensitivity if the DM density distribution features a core at the GC as small as 300 pc, which is currently allowed from the observational point of view (see e.g. [87,88]) and not tested by current numerical simulations (see e.g. [89][90][91][92]). Finally, the HESS reach from [85] does not probe qualitatively new mass regions, the highest DM mass probed being 10 TeV.
We do encourage the HESS collaboration, as well as other high energy cosmic-ray telescopes, to cast their limits also on secluded DM models, and to present their results in a format that could easily be reinterpreted in different DM scenarios.

Direct Detection
We compute the constraints from CRESST-II [8], CDMS-LITE [9], and LUX [11], each of which in turn gives the strongest bound for increasing M p D . The constraints for different values of are shown as yellow contours in Figs. 1 and 4. As r ∞ is decreased, the required α D is increased, and the direct detection constraints become stronger. Hence, in Fig. 4 the = 10 −10 contour appears for r ∞ = 10 −4 . The method used in setting the constraints from CRESST-II and LUX is given in Appendix B of [93] (also see [52,53,94]) 8 . We use a similar method in deriving the constraint for CDMS-LITE, the details of which are given in the appendix. We checked that the constraints from XENON1T [12] and PANDAX-II [13] are similar to those arising from LUX.

Dark atomic bound states
Due to the presence of the dark electrons, DM may form Hydrogen-like bound states. The radiative capture of p D and e D into a dark atom with emission of a mediator V is possible if the binding energy is larger than the mediator mass, If dark recombination is efficient in the early universe, then the amount of dark protons available to annihilate today may be significantly reduced, thereby suppressing the expected annihilation rate, as well as the direct detection signals. A proper treatment then requires computing the residual ionization fraction of DM today. This is beyond the scope of the present work. (The cosmology of atomic DM with a massless dark photon has been computed in detail in Ref. [65].) Instead, in Figs

Resonances and Reannihilation
It has recently been pointed out that, on the parametric resonances of the Sommerfeld enhancement factor, a two stage freeze-out may occur [62]. The first stage is the standard freeze-out and happens in the Coulomb limit of the Yukawa potential, at M p D /T ∼ 20 or v rel ∼ 0.1, in which the parametric resonance is not yet manifest. On the parametric resonances, the cross section scales as ∼ 1/v 2 rel for small velocities, before saturating. The DM annihilation then comes back into equilibrium, after the DM has kinetically decoupled, resulting in a second period of significant DM density depletion [95][96][97][98]. This second stage of the freeze-out process is known as reannihilation.
As our freeze-out calculation does not take into account the period of reannihilation, the values of α D found on resonance are slightly too large and would result in a DM under-abundance. The correct procedure would be to detune α D on these points, obtain an overabundance in the first stage of freeze-out, and then reach the correct relic abundance through reannihilation. Our derived constraints on resonance are therefore somewhat too strong. However, moving significantly off resonance means suppressing the reannihilation effect, as the cross section then grows only as ∼ 1/v rel at small velocities. The detuning away from the resonance peak is therefore necessarily small. Further detailed numerical work is required to determine both the amount of detuning and the quantitative impact of reannihilations on the annihilation cross section relevant for indirect detection.
Off resonance, our indirect detection constraints are more robust. As it was shown in [62], however, by choosing α D progressively closer to the peak on lower resonances, it can still be possible to obtain the correct relic abundance from multiple values of α D for the same (M p D , M V ) parameter point. Hence, also in this case, further work must be done in order to derive the constraints at these lower α D values. However, note that the indirect detection constraints are the strongest on resonance, so the smaller α D values may not open up considerable areas of parameter space. The quantitative determination of these effects, which were also ignored in [44,45], is left for future work.

Further Constraints on the Mediator
The results here have been presented in the M p D − M V plane. To confirm that SIDM is possible in this model, we must also compare the allowed areas to those on the M V − plane in Fig. 5 of [44]. For mediator masses relevant for SIDM, between 10 −3 M V /GeV 10 −1 , there are allowed areas centred around ∼ 10 −10 and spanning one or two orders of magnitude. Lower (higher) values of are ruled out by BBN (Supernova 1987A [99]) 9 . From Fig. 4, this range of M V is also allowed by direct detection experiments and ID, provided the asymmetry is large enough. Hence SIDM is indeed viable in this model.

On the assumption of early thermal equilibrium
All our results have been derived under the assumption of early thermal equilibrium between the dark sector and the SM, T D = T SM , for temperatures larger than any of the other scales we have considered. This is justified for values of the kinetic mixing 10 −6 , see e.g. [44,100], although this condition may weaken due to the presence of dark electrons. The values of allowed by direct detection, therefore, may not be large enough to realise early equilibrium in the region of parameter space relevant for self-interactions (see Figs. 1 and 4). In that region, the assumption of early thermal equilibrium could be justified by the presence of some extra dynamics that connects the SM and the dark sector. These dynamics could lie at an energy scale large enough to not alter the phenomenological study we have carried out, and would be even welcome to give a common origin to the baryon and Dark Matter asymmetries.
In the absence of such dynamics, an early T D /T SM = 1 would result in a freeze-out abundance different with respect to the one we have computed here. This would in turn modify the values of α D in our parameter space, and thus the regions corresponding to sizeable self-interactions and to exclusions from indirect and direct detection. A more detailed study of the impact of different thermal histories go beyond the purposes of this paper.

Conclusions
The focus of this paper has been to derive constraints on a relatively simple, well-motivated example of a dark sector, featuring fermionic DM and a light massive mediator, and allowing for the possibility of a dark asymmetry.
We have demonstrated that ADM coupled to light force mediators provides a viable framework for SIDM. While the combined direct, indirect and cosmological bounds severely limit the possibility of symmetric thermal-relic SIDM [44,45,67], large portions of the ADM parameter space evade the indirect constraints, due to the suppressed annihilation rate.
Nevertheless, we have shown that ADM with a Sommerfeld enhancement can yield significant annihilation signals that are constrained by various astroparticle messengers, i.e. γ-rays, antiprotons and neutrinos. These signals remain important even for large asymmetries, in particular late-time antiparticle-to-particle ratios as low as r ∞ ∼ 10 −4 , in our model. This challenges the standard lore that ADM cannot be probed via ID. Level transitions may give rise to additional radiative signals in this class of models [101][102][103][104]. Importantly, the ID constraints derived here extend beyond the SIDM regime. We encourage high energy cosmic-ray experiments to present their results in a way amiable for reinterpretation in secluded DM models, e.g. in terms of flux as a function of energy. Future improvements to direct detection experiments will be important in testing this scenario more thoroughly.
Asymmetric DM models that feature light or massless force carriers result typically in multicomponent DM, due to various considerations, including our arguments from gauge invariance and due to possible formation of stable bound states. Since ADM is a natural framework for SIDM, there is strong incentive for the development of multicomponent DM simulations.

A CDMS-lite
In setting the limit from CDMS-LITE [9], we use a similar technique as for CRESST-II. We use Poisson statistics to calculate the 90% CL on the number of nuclear recoils, assuming the events observed by CDMS-LITE between E min and E max are due to DM induced nuclear recoils. Here E min = 56 eV ee is the threshold energy of the detector and where v E (v Esc ) is the velocity of the earth (escape velocity) with respect to the DM halo, m T is the target mass, and m DM−T is the reduced mass of the DM and target. That is, E max is the either the maximum possible recoil energy given the DM mass or 1 keV ee -depending on which is smaller. The 1 keV ee cutoff is chosen because above this the large number of counts due to the 1.3 keV ee Lshell 71 Ge electron-capture decay line will begin to contribute as background. Note the CDMS-LITE events are reported in electron equivalent recoil energies. We convert the values into nuclear recoil energies, keV nr , using the central value of the Lindhardt model as used by CDMS-LITE (k = 0.157).
We then demand the expected number of events our model will contribute at a given point in parameter space, also taking into account the efficiency of CDMS-LITE, not exceed the 90% CL on the number of signal events we have calculated.