Late-Time Dark Matter Oscillations and the Core-Cusp Problem

The core-cusp problem persists as an unresolved tension between the predictions of $\Lambda$CDM cosmology and observations of dark matter (DM) profiles in dwarf spheroidal and other galaxies. We present a novel scenario for converting cusps into cores through reactivation of DM annihilation in galaxies at late times. This can happen in asymmetric DM models when there is a very small DM-number violating mass term that causes oscillations between DM and its antiparticle. Using analytic methods as well as gravitational N-body simulations, we show that this mechanism can robustly eliminate cusps from galactic DM profiles for light fermionic DM of mass $m_\chi\sim (0.1-1)$ GeV and a lighter mediator into which the DM can annihilate. We identify regions of parameter space where annihilation of DM particles is more efficient than elastic scattering at reducing the inner density of the DM profile. Dark matter annihilation is therefore a qualitatively distinct alternative to the mechanism of elastic self-interacting dark matter for addressing the cusp-core problem.

The core-cusp problem persists as an unresolved tension between the predictions of ΛCDM cosmology and observations of dark matter (DM) profiles in dwarf spheroidal and other galaxies. We present a novel scenario for converting cusps into cores through reactivation of DM annihilation in galaxies at late times. This can happen in asymmetric DM models when there is a very small DM-number violating mass term that causes oscillations between DM and its antiparticle. Using analytic methods as well as gravitational N-body simulations, we show that this mechanism can robustly eliminate cusps from galactic DM profiles for light fermionic DM of mass mχ ∼ (0.1 − 1) GeV and a lighter mediator into which the DM can annihilate. We identify regions of parameter space where annihilation of DM particles is more efficient than elastic scattering at reducing the inner density of the DM profile. Dark matter annihilation is therefore a qualitatively distinct alternative to the mechanism of elastic self-interacting dark matter for addressing the cusp-core problem.

I. INTRODUCTION
In many respects the standard ΛCDM paradigm of cosmology gives an extremely good description of the observed universe. But it has long been recognized that simulations of structure formation that neglect the presence of baryons predict singular (cuspy) density profiles of the dark matter (DM) toward the centers of galaxies, whereas observations suggest flatter (cored) distributions. More recent simulations include the effects of baryonic feedback, which can expel material from denser regions and help to ameliorate this discrepancy, but there is not yet any consensus that this provides a complete solution. Moreover in systems like dwarf spheroidals, where baryons are relatively scarce, one does not expect baryons to have a significant impact on the small scale structure. These issues have been reviewed in ref. [1].
Another proposed solution is self-interacting dark matter (SIDM) [2], with a scattering cross section at the level of that is close to upper bounds from colliding galaxy clusters, such as the Bullet Cluster [3], even though it may be not so robust [4]. N-body simulations incorporating such interactions have shown that cross sections consistent with eq. (1) can produce cored DM profiles in a wide range of systems [5,6]. However, more recent studies indicate that a constant cross section is not the ideal solution, since then σv increases with the size of the system, contrary to the observation that cores are less pronounced on the scales of galactic clusters. A weak velocity dependence of the form σ ∼ 1/v is found to give a better fit to the full range of structures [7]. The most common assumption is that the DM selfinteraction is in the form of elastic scattering, but a more exotic possibility was proposed in ref. [8], in which fusion of DM particles into bound objects is the interaction leading to cored profiles. Like other exothermic processes, this has the advantage of predicting a cross section with σv remaining constant at low velocities, as desired for fitting the DM profiles of both large and small galactic structures. Other interesting possibilities to achieve the correct velocity-dependence of SIDM have been studied, for instance in the context of resonant SIDM [9][10][11], puffy DM [12], self-heating DM [13], maximally SIDM [14] and DM bound states produced in the early universe by threebody recombination [15,16].
Here we explore a different alternative, motivated by the fact that DM annihilation is also an exothermic process with σv becoming constant as v → 0. The challenge for such a scenario is to explain how annihilations could go out of equilibrium in the early universe, but then come back at late times [17]. In fact, a mechanism to do this is well known in the context of asymmetric dark matter, where there is an asymmetry between the DM χ and its antiparticleχ. By allowing for a small mass term that violates the conservation of DM number, oscillations between χ andχ can reactivate the annihilations at late times [18][19][20][21].
The reactivation of DM annihilation at late times is usually seen as a danger to be avoided, since it is known that the DM density should not change appreciably between the era of the CMB (redshift z ∼ 1100) and structure formation [22], but in the present work we demonstrate that this mechanism can efficiently produce cored profiles in galaxies without changing the total DM density significantly. The reason is that the efficiency of oscillations leading to regeneration of the anti-DM component can depend strongly on density, so that annihilations are effective in the centers of galaxies but not in the outer regions.
For a more quantitative investigation, one should integrate quantum Boltzmann equations for the density matrix, that account for the coherence of states undergoing oscillations, analogous to those used for the study of neutrino oscillations in a medium. This formalism was initially worked out for DM in ref. [20], and some important corrections were realized in ref. [21], which we follow closely in the present work.
We consider two models of quasi-Dirac fermionic DM χ of mass m χ . In the first, the dark matter couples to a lighter vector boson V µ (Model 1), with effective Lagrangian In Model 1, the dark matter freeze-out and the late-time depletion are both allowed by the annihilation process χχ → V V . In the second model, dark matter couples to a complex scalar Φ = φ + ia (Model 2), Model 2 allows freeze-out and late-time depletion from χχ → φa (which unlike χχ → φφ or χχ → aa is s-wave, hence not suppressed at low velocities). The coupling between χ and either kind of boson is denoted as g , and its associated fine-structure constant is α = g 2 /4π. The DM-violating mass term is The parameter δm violates not only dark matter number, but also the gauge symmetry of Model 1, which is additionally broken by the Stueckelberg mass term for the vector. It would be possible to replace both of these explicit breakings by a Higgs mechanism, but for simplicity we adopt the simpler effective theory. We begin by making preliminary estimates to identify viable regions of the parameter space, in section II. The essential details of the density matrix Boltzmann equation formalism are reviewed in section III. In section IV we will show that, for appropriate choices of the model parameters, integration of the Boltzmann equations in the early universe leads to the conventional freeze-out of DM annihilations, leaving only the asymmetric component of the DM. This justifies the initial conditions for the second step, described in section V, where we re-solve the analogous Boltzmann equations in the background of an already-formed galaxy and show how an initial cusp gets erased by reactivated annihilation following χ-χ oscillations. This is a somewhat crude approach since it considers formation of the galaxy to happen suddenly and neglects the role of gravity in shaping the DM halo. In section VI, we improve on this by carrying out a gravitational N-body simulation of galaxy evolution, in a code adapted to properly account for the new physics effects. Section VII is devoted to a discussion of the possible constraints that arise from cosmological observables, such as the CMB and Big Bang Nucleosynthesis (BBN). Conclusions are given in section VIII, and details of the quantum Boltzmann and N-body simulation methods are presented in the appendices.

II. ANALYTIC ESTIMATES
Before embarking on a detailed analysis, we analytically estimate the regions of parameter space that are of interest for our mechanism. First, the annihilation cross sections at threshold for the two models are σv a = π α 2 m 2 χ × (1 − r 2 m ) 3/2 /(1 − r 2 m /2) 2 , Model 1 (1 − r 2 m /4), Model 2 (5) where r m is the ratio of the mediator to the DM mass, r m = m V /m χ for χχ → V V (Model 1) or r m = m φ /m χ for χχ → φa (Model 2, where we assume m a m φ ). We neglect the p-wave suppressed channels χχ → φφ, aa.
To compare this to the desired cross section (1), consider a reference velocity v 0 = 100 km/s characteristic of DM in a Milky-Way-like galaxy, and the upper value in the range (1), giving σv/m χ ∼ 100 cm 2 km/s/g. Equating this to σv a suggests the constraint For example with m χ = 100 MeV, α ∼ = 0.02; we will adopt these as approximate benchmark values. However nothing prevents us from taking somewhat heavier DM, up to m χ ∼ 1 GeV; above this, the theory starts to be strongly coupled.
It is impossible to avoid χχ elastic scattering mediated by the annihilation products, and we choose to constrain these cross sections so that they are below the level that would change the DM density profile independently of the annihilation effect, which is the focus of this work. The elastic scattering cross sections at low velocities are where v = v rel /2 is the center-of-mass velocity. To avoid that the scattering self-interactions play a leading role in the galactic dynamics, we require that σ s v rel σ a v . This implies (m V,φ /m χ ) 4 4v rel , which is most stringent for large systems, galaxy clusters, that have the highest DM velocities. For example, the cluster A2537 has velocity dispersion σ v ∼ 1000 km/s [23], with v rel = (4/ √ π)σ v (assuming the velocity is Maxwell-Boltzmann distributed), and demanding that σ s v rel < 0.3 σ a v gives the constraints 0.6 < r m < 0.94, Model 1 0.6 < r m < 1.99, Model 2 (8) where the upper limits come about because of phase space suppression of the annihilation. 1 It is consistent to take m a /m φ nearly as small as (v/c) ∼ 0.003 (the DM velocity in clusters) while still ignoring its d-wave suppressed contribution to scattering in (7). But to satisfy cosmological constraints we assume m a 1 MeV; see section VII. The χ-number violating mass δm must be small enough so that χ-χ oscillations have not yet started at the time of DM freeze-out, T χ,fo ∼ m χ /20, where we allow for a lower temperature in the dark sector, as discussed in more detail in Sec. VII. For annihilations to recouple during structure formation, the oscillations should start before ∼ 0.1 Gyr. In this way we identify a large range of possible values where the last estimate assumes m χ ∼ 100 MeV. The upper bound in eq. (10) comes from requiring that oscillations are not coupled at the time that annihilations of χ particles are freezing out. Since n χ σv a ∝ a −3 whereas H ∝ a −2 before the formation of the CMB, if n χ σv a < H during annihilation freeze-out then n χ σv a /H will remain smaller than 1 until after overdensities go nonlinear. Regardless of the precise value of the upper limit, it is clear from eq. (10) that there is a wide range of parameters for which annihilations can decouple in the early universe, while still being important in overdense environments at late times.

III. OSCILLATION FORMALISM
In the presence of DM oscillations, the distinction between particle and antiparticle becomes time-dependent. If we define a basis then it is straightforward to show that the time dependence of a state that is initially pure |χ is with c ϕ = cos ϕ, s ϕ = sin ϕ, ϕ = δm t. To this state we can associate the density matrix for a single-particle state, Naively it might seem like appreciable amounts ofχ appear as soon as ϕ ∼ 1 and χχ annihilations could recommence, but this need not be true if all the particles in the plasma are oscillating with the same phase. Ref. [21] showed that recoupling of annihilation depends on the nature of the interactions. Interactions of fermionic DM with vectors V are called "flavor sensitive," while interactions of χ with scalars or pseudoscalars are "flavor blind," leading to very different behaviors of the annihilation probabilities. In the collision of two particles with respective phases ϕ and ϕ , the annihilation rates are modulated by the factors In the first case, a bath starting as pure |χ and maintaining phase coherence never undergoes annihilations since ϕ − ϕ remains zero, despite the oscillations. In the second, the modulation factor averages to 1/2 for fast oscillations, and is therefore effective even when the particles stay in phase with each other. For a thermal bath, the matrix n 1 in eq. (13) is replaced by an integral over the corresponding matrix distribution function F(k) for the states of momentum k, where s = 1/2 for fermions as we consider. Then n 11 (n 22 ) represents the number density of particles (antiparticles) defined with respect to the basis {|χ , |χ } as in (11); the off-diagonal elements keep track of the coherence between these two states. The Boltzmann equation for n reduces to the usual form when we consider only the diagonal elements, but it has additional terms due to the off-diagonal elements, which depend on whether the interactions are flavorsensitive or flavor-blind.

III.1. Model 1: vector mediator
We first consider the flavor-sensitive case, applicable to Model 1, for which the Boltzmann equation iṡ where H is the Hubble parameter, the thermally averaged free Hamiltonian is σv s is the χχ or χχ scattering cross section (that coincide at low energies), σv a is the χχ → V V annihilation cross section, and n eq is the equilibrium number density. The scattering term in eq. (17) is derived in appendix A, while the other terms can be found in ref. [21]. 2 The scattering term in (17) has the effect of damping the off-diagonal elements of n, which destroys the coherence of the quantum superpositions and effectively measures the state of an oscillating system. The loss of coherence results in det n = 0, which activates the annihilations. Since σv s is proportional to the DM velocity, this makes the effect stronger in systems with large velocity dispersions. We will see in section V that this is contrary to observations, disfavoring Model 1 taken by itself.
The origin of the factor (14) can be heuristically understood from (17) by interpreting the annihilation term det n 1 as the matrix [21] det n 1 → 1 where n 1 , n 2 represent the density matrices (13) of two particles, having respective phases ϕ 1 , ϕ 2 , and σ 2 is the Pauli matrix. 3 Similarly, the effect of the scattering term can be understood by replacing n → n 1 everywhere except in the trace, where Tr n → Tr n 2 , which does not depend on ϕ 2 , and just represents the total density n of DM scattering on particle 1. Then the off-diagonal parts of the Boltzmann equation determine the damping of ϕ 1 as which has the solution c ϕ s ϕ ∼ exp(− 3 2 Γ s t)(c ϕ s ϕ ) 0 , where Γ s = n σv s is the elastic scattering rate.

III.2. Model 2: scalar mediators
For scalar interactions, the Boltzmann equation simplifies, because elastic scattering no longer has any effect on the density matrix. The form of the annihilation term is also changed, in a way that makes it lead to decoherence by itselḟ − σv a det n (Tr n)n 21 (Tr n)n 12 det n − n 2 eq 1 .
Here we define det n ≡ n 11 n 22 + n 21 n 21 . In contrast to eq. (17) for the vector model, there is no dependence on the DM velocity in (21), leading one to expect a more similar level of cusp erasure in both large and small galactic systems, independently of the differences in their velocity dispersions.
The diagonal term goes to (15) when the two phases are different from each other. The off-diagonal term leads to phase damping similarly to (20), but now with c ϕ s ϕ ∼ exp(−Γ a t)(c ϕ s ϕ ) 0 , where Γ a is the conventional annihilation rate, without the sin 2 2ϕ modulation factor.

IV. EARLY COSMOLOGY
For very small values of δm 10 −30 eV, oscillations are unimportant until the epoch when structure formation begins. For larger values of δm they can cause annihilations to temporarily recouple, further reducing the density of the asymmetric component, before structure formation begins and annihilations are reactivated once again. In this section we illustrate these possibilities by solving the Boltzmann equation at early times. This is meant to provide the initial conditions before the effects of oscillations on structure formation begin, that we will investigate in the following sections.

IV.1. Model 1
Like for conventional freeze-out, it is convenient to use x ≡ m χ /T as the independent variable, and the abundance Y ≡ n/s as the dependent variable, where s = 2π 2 m 3 χ /(45g * s x 3 ) is the entropy density and Y is now a matrix. The Boltzmann equation becomes Here is the Hubble parameter, and we have allowed for the DM temperature to differ from that of the standard model by putting the appropriate factors of ξ = T χ /T . The averaged scattering cross section is and the equilibrium abundance is in the Maxwell-Boltzmann approximation. It is the abundance of just the DM particle χ, not including the antiparticle. There is an additional possible source of decoherence that is not captured by eq. (23). The full Hamiltonian (18), before taking the non-relativistic limit, depends on the momentum k of the state, which is neglected in (23). This causes states of different momenta to oscillate at slightly different frequencies δω ∼ δm (k 2 /2m 2 χ ), giving rise to thermal decoherence even in the absence of scattering. To fully investigate this effect would require solving for the full distribution function F(k), which is numerically prohibitive. Instead we model it in an approximate way, by splitting the integral over k in (16) into two bins of small and large momenta, Y = Y s + Y l . The averaged Hamiltonians H s,l for the respective bins are shown in eq. (B7). The resulting coupled Boltzmann equations are given in eqs. (B3). They have a more complicated matrix structure than (23), but the sum of the two agrees with (23). We found this additional source of decoherence to have a negligible effect, compared to that due to the scatterings.
The left panel of Fig. 1 shows the effect of thermal decoherence in the evolution of the oscillating dark matter in the vector model at early, intermediate, and late times. We see that after oscillations commence at late times (values of x ∼ 10 7 − 10 8 in the two models), annihilations recouple briefly before freezing out again. The dark matter density, Y = Yχ + Y χ , is reduced by O(5%), which is roughly compatible with observations of density perturbations in the CMB [22], and such a change in the dark matter abundance may even have some implications for the H 0 [24] and σ 8 [25] tensions. We leave further exploration of these implications for future work.

IV.2. Model 2
The cosmological version of the Boltzmann equation (21) is where scatterings no longer play any role. Its solution is shown in the right panel of Fig. 1. The implications of the brief recoupling in dark matter annihilation are similar as in Model 1.

V. STRUCTURE FORMATION
We start with an approximate treatment of the effect of χ-χ oscillations on galactic dynamics, by imagining that an NFW-shaped halo with has already formed at some time t 0 , with the initial condition on the matrix density that at each position r in the collapsed system: this corresponds to a pure χ state, in which oscillations have not yet had any effect.
To apply the Boltzmann equation (17) in a galactic environment that has separated from the Hubble expansion, we drop the 3Hn term, and set n eq = 0, since the annihilation products are assumed to escape. We evolve the initial density (28) at each radial position r, up to a final time t f of order 10 Gyr. This leads to a modified density profile ρ χ = (n 11 + n 22 ) m χ , that is to be compared to present-day observations.
In addition to knowing the initial density profile, it is also necessary (for Model 1 only) to specify the DM velocity profile, since the relative velocity enters into the scattering rate through σv s (whereas the annihilation rate is insensitive to v rel ). We have adopted the analytic solution for the radial velocity dispersion σ r (r) from ref. [27] (see eq. (14) of that reference), derived by solving the Jeans equation for an NFW profile. This determines σ r (r) for given NFW parameters r s , ρ s . The latter can be related to the virial radius where ρ c is the present critical density, and g(c) = [ln(1+ c) − c/(1 + c)] −1 .

V.1. Model 1
We applied this procedure first within Model 1 for a particular dwarf spheroidal galaxy, DDO 154 [28], that has been discussed from the point of view of selfinteracting dark matter in ref. [26]. There the NFW parameters were determined using data from ref. [29] and the mass-concentration relation from ref. [30]. The resulting NFW profile is shown in Fig. 2 (top left). This profile disagrees with the observed rotation curve in the inner part of the galaxy, whereas the solid "SIDM" curve, which arises from elastic DM self-interactions adjusted to the appropriate cross section, gives a good fit.
The dot-dashed curves show the results for our model, with m χ = 100 MeV, α = 0.02, and several values of m V . The profile is significantly cored, depending on the value of m V , and has a different shape from that predicted by SIDM. The closest match between the SIDM profile and ours is produced for m V = 26 MeV, which however is inconsistent with the constraint (8). It means that we should not neglect the effects of elastic scattering by itself, which go into the usual SIDM treatment. This problem can be overcome by simultaneously increasing m V and α ; for example m V = 60 MeV and α = 0.1 gives a reasonable fit. However neither of these models are consistent with data from galaxy clusters, as we discuss next.
Ref. [7] presents evidence for the DM profiles of galaxy clusters also being cored, to a somewhat lesser extent than dwarf spheroidal galaxies. These larger systems have much higher velocity dispersions, which leads to a stronger reduction of the central density by our mechanism, using Model 1. This is shown for the cluster A2537 in Fig. 2 (top right), where for the same values m χ = 100 MeV and α = 0.02 as before, the best match to the SIDM curve is for m V ∼ = 51 MeV; the lower value of m V = 26 MeV, favored by dwarf spheroidals, leads to unacceptably large suppression of the central density to be compatible with measured stellar velocities profiles. More detailed quantitative comparisons between the theory and data will be presented in section VI, in terms of the predicted versus observed velocity profiles.
One may also wonder to what extent a given model can match the observed properties of different spheroidal dwarf galaxies, whose density profiles can be diverse. Although an exhaustive comparison is beyond the scope of the present work, we have studied a contrasting example, DDO 126, whose DM density profile (like that of DDO 154) was estimated by ref. [29]. The best fits to the circular velocity measurements for the two galaxies occur at different values of the model parameters, as shown in the left panel of Fig. 3, where we fixed m χ = 65 MeV, α = 0.015 and allowed m V to vary. (Notice we have chosen a lower value of m χ in this example; it is motivated by the discussion in section V.3.) However, an acceptable fit to both systems can be found at an intermediate value m V ∼ = 20.6 MeV, resulting in χ 2 /d.o.f. ∼ = 0.8 for either system. We have allowed for systematic uncertainty in the magnitude of the DM density profiles, reflecting an estimated ∼ 25% uncertainty in the baryonic content of the galaxies [7]. Since the baryons comprise ∼ 10% of these systems, this translates to a 2.5% uncertainty in the overall DM densities, that we have marginalized over to slightly improve the fits.

V.2. Model 2
In Model 2, the situation is the opposite, though with a smaller discrepancy. In this case nothing depends on the scalar mass m φ , as long as it satisfies the consistency condition (8). For a fixed value of m χ (here still at 100 MeV), only α matters. The SIDM profile can be approximately matched by taking α ∼ = 0.01 in the DDO 154 dwarf galaxy, while for the same parameter choices, the predicted inner profile of cluster A2537 lies somewhat above the SIDM fit, a factor of 1.7 higher as illustrated in the density profiles shown in Fig. 2.
We have found that this qualitative difference between scalar and vector mediators is generic: the velocity dependence of the decoherence mechanism in Model . Left: predicted circular velocities due to the DM component alone from the same two models and from SIDM (ref. [26]), and data from ref. [29]. In each case, one mediator dominates the coring effect of the central profile in one system, while having little effect in the other system. Right: stellar velocity dispersion along the line-of-sight for cluster A2537, with predictions based on the DM density profile from two of our models, from SIDM (ref. [7]) and data from ref. [23].
1 makes it more effective for cusp suppression in highvelocity systems (clusters), whereas the lack of such dependence in Model 2 leads it to be more efficient in higher density systems (dwarfs).
The mild tension in simultaneously explaining the density profiles of different spheroidal dwarf galaxies, described for Model 1, is also present in model 2, as illustrated in the right panel of Fig. 3 for the case of m χ = 65 MeV: the best fits occur at different values of α for the DDO 154 and DDO 126 galaxies. Like for Model 1, it is not a serious difficulty since an intermediate choice α ∼ = 0.0053 results in an acceptable χ 2 /d.o.f. = 0.72 for both systems. We leave a more exhaustive study, both of the allowed parameter space and including more galaxies, for future work.

V.3. Hybrid models
The previous results suggest that the challenges for Models 1 or 2 to simultaneously fit the rotation curves of both dwarf galaxies and clusters could be overcome in a model with both mediators present. Here we present an example that supports this hypothesis, leaving for future work a more rigorous or detailed analysis.
Since it is technically difficult to implement both kinds of mediators simultaneously, we will be content here to give an example in which a vector mediator gives a good fit to a cluster, while leaving a dwarf galaxy relatively unaffected, and at the same time a scalar mediator that achieves the opposite. Since each model has a relatively small effect on one of the systems, it seems likely that by combining them, one can add the coring effects to both systems in a roughly linear fashion.
For example we find that for lighter DM with m χ = 65 MeV, and Model 1 parameters α = 0.015, m V = 44 MeV, we fit the observed stellar line-of-sight velocity dispersion profile (described in more detail in the next section) for A2537 extremely well, while leaving the predicted circular velocity V circ (r) in DDO 154 too high from not sufficiently reducing the central density in the dwarf system. On the other hand, choosing a lower coupling α = 0.004 in Model 2 gives an excellent fit to the DDO 154 rotation curve, while having a small impact on the inner profile of A2537. These outcomes are shown in Fig. 4, indicating that by combining the two mediators, it is possible to get as good a fit as an elastic SIDM model with a velocity-dependent cross section that is tuned to fit both systems. In elastic SIDM, a cross section of σ/m ∼ = 3 cm 2 /g [26] is needed to agree with dwarf spheroidals, whereas a smaller value ∼ 0.1 cm 2 /g is used to explain clusters [7].

VI. N-BODY SIMULATIONS
To more quantitatively predict the evolution of galactic structures in our scenario, we have performed N-body simulations that take into account the peculiar interactions described by the Boltzmann approach of the previous sections. The two approaches should be viewed as complementary since each has its own limitations. The challenge for N-body simulations, even if modified to account for self-interactions, is that they treat test particles classically, with scatterings occurring probabilistically rather than quantum mechanically. For conventional self-interactions this is not a serious limitation, but in the present context, the overall coherence of the DM ensemble is of primary importance.
To address this, we have modified the public version of the GADGET-2 code [31,32], which is widely used to generate N-body cosmological simulations. 4 The novel feature, apart from including DM scattering and annihilation (see appendix C for implementation details and code tests), is to keep track of the phase ϕ of each test particle, that describes the oscillations as in eq. (12). We assume that all particles in the halo are initially in phase with each other. Depending on whether the model is flavor-sensitive (Model 1) or flavor-blind (Model 2), this phase plays different roles, and evolves differently. In the absence of interactions, the phase of each particle would evolve trivially as ϕ = δm t. To mock up the behavior predicted by the quantum Boltzmann equations while still treating the particles classically, we implement scattering as follows.
Model 1. Elastic scatterings damp the quantum coher-4 https://wwwmpa.mpa-garching.mpg.de/gadget/ ence, as described by the off-diagonal elements of the collision term in (17). Integrating the off-diagonal elements over a collision time ∆t = 1/Γ s = (n σv s ) −1 leads to a phase change ∆ ln(c ϕ s ϕ ) = −3/2, as shown in eq. (20). For strongly damped systems such that Γ s > δm, this can be modeled by replacing the phase of each particle undergoing elastic scattering by leading to decoherence of the ensemble, that allows annihilations to occur. The annihilation probability of two particles with respective phases ϕ 1 and ϕ 2 is reduced relative to its usual value by the factor sin 2 (ϕ 1 − ϕ 2 ), as derived in eq. (19).
Model 2. In this case, the scattering self-interactions have no effect on the phases, and they play exactly the same role as in conventional SIDM. Instead, decoherence is caused by the annihilation interactions themselves. The phase reduction described above now becomes a fac-tor of e −1 each time an annihilation would have occurred, for a fully decoherent mixture of χ andχ. The annihilation probability is modulated by the different factor sin 2 (ϕ 1 + ϕ 2 ) as was explained below eq. (22).
As initial conditions for DM halos corresponding to the dwarf galaxy DDO 154 and the galaxy cluster A2537, we took Hernquist profiles [33], which are described by the total mass M and the scale radius a. Unlike NFW, these profiles have finite mass without any need of truncation and they are perfectly stable in time when evolved with collisionless DM [34,35], as we show in appendix C.
To match the initial N-body profiles to the ones assumed in section V, we used the procedure described in ref. [34]. It consists of choosing the value of the Hernquist mass M as the virial mass M 200 of the NFW profile and requiring the two density profiles to coincide in the inner region where r r 200 . The latter condition gives a relation between the Hernquist scale radius a and the NFW one r s , namely  Fig. 2. The overall agreement observed for both DM models suggests that the Nbody simulations model reasonably well the physics encapsulated in the quantum Boltzmann equation, where the coherence of DM particles plays a decisive role. Differences between the simulation and the approximate approaches are perhaps more evident in the dwarf galaxy because gravitational effects and DM dynamics have a relatively larger effect in small systems than in large ones.
To compare our model predictions with existing data, we converted our results for the DM density into observed quantities, namely the circular velocity for dwarf systems and the projected stellar velocity along the lineof-sight for galaxy clusters. The former is defined as V circ (r) = [G M (r)/r] 1/2 , where M (r) is the enclosed mass at radius r. The left panels of Fig. 6 show our DM predictions for the rotation curve of DDO 154 dwarf galaxy within the two classes of models considered in this paper compared to current data. The grey points show the total circular velocity of the dwarf as observed by the LITTLE THINGS survey [29], whereas the white points represent just the DM contribution to V circ (r), obtained by subtracting the gas and star components after carefully modelling their distribution within the galaxy [29].
The vector model with m V = 34 MeV provides the best fit to data among the models displayed in the top panel of Fig. 6, with a χ 2 /d.o.f. < 1, comparable to the SIDM curve found in ref. [26]. For the scalar case, a choice of α somewhat smaller than 0.01 would provide good agreement between our model and observations as shown in the bottom left panel of the same figure.
For relaxed clusters dominated by a central early-type galaxy, such as in A2537, it is possible to measure the stellar line-of-sight velocity dispersion profiles σ LOS (r) with spatially-resolved spectroscopy [23,36]. In order to convert our model predictions into σ LOS (r), we used the procedure outlined in appendix A of ref. [37] combined with the information for A2537 cluster contained in ref. [23]. In particular, as done in the latter reference, we modeled the stellar luminosity density ν (r) with a dual pseudo isothermal elliptical profile (dPIE) [38] and converted it into a baryonic density via the relation Here Υ V ≡ M /L V is the stellar mass-to-light ratio in the V-luminosity band, which is usually assumed to be spatially-independent across the cluster [23,37]. The value of Υ V could be inferred from the stellar population synthesis (SPS) up to an unknown initial mass function (IMF) and therefore one usually parametrizes this ignorance with the free parameter log (Υ V /Υ SPS V ), where Υ SPS V is the SPS predicted mass-to-light ratio for a given IMF. We considered a Chabrier IMF [39] as done in ref. [23] and fixed the value of log (Υ V /Υ SPS V ) for the A2537 cluster by matching the baryonic density computed by eq. (32) with that obtained in ref. [7]. The results of this procedure are shown in the right panels of Fig. 6 for both the vector and scalar models. Good agreement between them and the existing data is obtained for a wide set of parameters in both classes of models because of the large error bars in the observational data.
The N-body approach allows us to distinguish between the complementary effects of ordinary self-interactions by scattering, versus the novel one from annihilations, which we have investigated in both Models 1 and 2. To estimate the annihilation contribution to the total profile, we turned off the elastic scattering processes. Similarly, the scattering contribution can be estimated by turning off the annihilations. Fig. 7 shows that the major effect in shaping the halo density profile is given by DM annihilation for the choices of parameters both in Model 1 and Model 2 considered in this paper. This verifies that the annihilation mechanism, investigated here for the first time in quantitative detail, has the capacity to alleviate the small-scale structure problems of CDM in the way originally suggested by [17].
Using dark matter annihilation to solve the core-cusp problem naturally gives a roughly constant value of the rate of core formation [8], as is suggested by fits to astrophysical objects spanning five decades in mass [7]. Relying on dark matter dynamics to resolve these issues is potentially under some tension from the measurement of  [29]. In particular, the grey dots show the total effect of DM, gas and stars on the rotation curve, whereas the white dots show just the DM contribution obtained after a careful modelling of stars and gas components (see ref. [29] for details). Right: projected stellar velocity dispersion along the line-of-sight as a function of radial distance for the cluster A2537. The data points and the error bars are taken from ref. [23]. In all panels, N-body simulation results are shown as solid lines surrounded by the 1σ uncertainty band, obtained by assuming that the number of particles in each bin is Poisson-distributed. The black dotted curve corresponds to the original NFW profile, whereas the matched Hernquist profile is shown with the red dashed line. The other dot-dashed curves are the results of Fig. 2. The orange solid line is the SIDM prediction from ref. [26] for DDO 154 and from ref. [7] for A2537. cusps in the centers of classical dwarf spheroidal galaxies [40,41] and from recently discovered ultrafaint galaxies [42], although out-of-equilibrium dynamics like tidal effects of the host galaxy may play a role in contributing to the diversity of these systems [43][44][45][46][47]. Doing selfconsistent fits to the observational data across many different systems will be critical for determining if the mechanism we investigate in this paper is as quantitatively successful as the elastic SIDM mechanism has been. Exploring this model in cosmological N-body simulations to compare against the subhalo abundance, for instance, will also be an important route for future work.

VII. COSMOLOGICAL CONSTRAINTS
The framework we propose necessarily introduces a hidden sector having several light degrees of freedom, in particular bosonic mediators like the vector V of Model 1 or the scalars φ and a of Model 2. If these coupled with appreciable strength to the Standard Model, this model would produce an indirect detection signal approximately nine orders of magnitude larger than current limits from X-ray and gamma-ray telescopes. Thus, the hidden sector must be largely secluded from the Standard Model, and we must clarify how the entropy of the dark sector dilutes.
If the V, φ, and a force mediators, which are the products of χχ annihilation, were exactly stable, they would introduce a significant matter component to the universe well before the time of the CMB formation. Therefore, they must decay into radiation of some sort to avoid dominating the energy density of the universe at low temperatures, for example at the time of big bang nucleosynthesis (BBN). The simplest solution is to introduce a dark radiation species, such as massless sterile neutrinos ν , that couple to the mediators, and allow for the decays V, φ, a → ν ν . As long as the dark neutrinos have a mass smaller than ∼ eV, this ensures they do not come to matter-dominate the universe before the formation of the CMB. Even this single light degree of freedom might be detected by precise probes of the energy content of the early universe at BBN [48,49] and CMB [50], which constrain the number of new relativistic degrees of freedom. For single-parameter extensions of ΛCDM, the constraints are of order ∆N eff 0.2, but known parameter degeneracy with the helium abundance Y p and possible hints of beyond-ΛCDM physics such as neutrino masses or the H 0 tension can partially relax these constraints to the level of ∆N eff 0.5 [48][49][50]. 5 In the near future, highresolution studies of the CMB damping tail will improve these bounds by an order of magnitude [51]. These constraints are all substantially stronger than the 1.75 (3.5) degrees of freedom contributed by a fully thermalized Majorana (Dirac) fermion. 5 see Eqs. 68-69 and 81 of ref. [50] or Fig. 11 and Table 5 of ref. [49].
However, as specified above, our hidden sector must be rather well secluded from the Standard Model bath. This can reduce the entropy carried by each dark-sector degree of freedom in comparison to that carried by Standard Model particles. Conservation of comoving entropy indicates that for any given thermal sector, the temperature evolves as T = (g * S (T i )/g * S (T )) 1/3 T i , where g * S is the number of degrees of freedom in entropy, s ∝ g * S T 3 . If we allow for an initial discrepancy between the Standard Model and dark sector temperatures, T d,0 = ξ 0 T γ,0 , the contribution to the effective number of neutrinos at any temperature is given by where g d * , g * S , and g d * S are all evaluated at T γ . In going from the first to the second line, the Standard Model degrees of freedom in entropy before and after e + e − freezeout cancels the change in neutrino temperature, as expected. We normalize to the values appropriate for a dark sector containing one light and one heavy Dirac fermion, a vector, and a charged scalar (as needed to give vectors a mass), reflecting the field content of Model 1. If instead we have a dark sector with two heavy and one light Majorana fermion plus a complex scalar, as in the minimal case to generate the phenomenology of Model 2, we would have ∆N eff 0.31 ξ 4 0 . In either case, the CMB and BBN limits in single-parameter extensions of ΛCDM are in tension with the model for ξ 0 = 1, but the model may be reconciled to the most stringent BBN and CMB limits for ξ 0 0.9, which can be arranged with only mild discrepancies in inflationary reheating widths to the two sectors [52,53]. Regardless of this caveat, either model is compatible with CMB and BBN limits once uncertainties in Y p or other quantities are taken into account [48][49][50].

VIII. CONCLUSIONS
The long-standing discrepancies between gravitational N-body simulations of structure formation in the ΛCDM paradigm and observations of cored density profiles continue to motivate exploration of alternative dark matter models and mechanisms. In this work we have revived one of the earliest such proposals [17] by showing that dark matter annihilations in galactic structures can be responsible for erasure of the cusps, using distinctive properties of asymmetric dark matter (ADM). The key idea is that very strong annihilations would freeze out early in cosmic history, solving the problem of removing the "symmetric" ADM relic density, and are reactivated at late times relevant for structure formation by oscillations of DM into its antiparticle. The preferred annihilation rate per unit mass σv/m χ ∼ 100 cm 2 /g km/s can be fit in our model by dark matter and mediator masses of order 30 MeV m V,φ,a m χ 100 MeV, a perturbative self-coupling as given in Eq. (6), and a Majorana mass term δm within the range (10 −31 − 10 −14 ) eV.
To obtain a large-enough annihilation cross section while respecting perturbativity of couplings constrains the DM and the mediator of the strong hidden force to be light, typically below 100 MeV. We have illustrated the mechanism in two representative models, with vector or scalar mediators respectively, and using two complementary approaches to model the structure formation dynamics. A fully consistent simulation is challenging because it must incorporate the quantum coherence of the oscillating dark matter while tracking the spatiallydependent annihilation rates within a DM halo. Our Nbody simulations, which treat the coherence in an ap-proximate way, give relatively close results to a quantum Boltzmann equation approach, which models the structure formation less rigorously. We have tested the scenario on two representative dwarf spheroidal galaxies, as well as a galactic cluster. Both methods lead to significant coring of the density profile, qualitatively similar to the effects of elastic SIDM scattering that have been widely used to address the cusp-core problem.
Like the elastic SIDM paradigm, the new mechanism we propose here does not, in its simplest forms, address the diversity of halo profiles on all scales. In elastic SIDM this is accomplished by assuming velocitydependent scattering, with a cross section that goes down at larger DM speeds. Within our mechanism, scalar mediators generically have a relatively stronger coring effect on small halos than larger (less dense) ones, while vector mediators have the opposite behavior. We presented evidence that the combination of both mediators could provide a good universal fit, leaving a detailed investigation for future study.
Acknowledgments. We thank A. Benson, S. Tulin and A. Robertson   In this appendix we derive the collision term for elastic scattering of χχ or χχ through exchange of the vector boson, needed for the quantum Boltzmann equations. The diagrams in Fig. 8 are the analog of Fig. 4b in ref. [21]. We can read off the imaginary part of the self-energies Σ >,< , in analogy to their eq. (A26) of ref. [21], where dp = d 4 p/(2π) 4 and the Green's functions are given by HereF is the matrix with the diagonal entries interchanged, as in (11) of [21], and O − = diag(1, −1). The trace is over both Dirac and flavor indices. We also definẽ This has the effect of reversing the signs of the off-diagonal elements.
Considering the relevant physical processes, it is not necessary to take account of all eight of the terms that arise from each diagram, from the products of the S >,< functions. First, since annihilation diagrams are suppressed while k 0 > 0, we can ignore k 0 < 0, which would give the s-channel diagram. Second, by energy conservation, we must have either p 0 > 0 and p 0 > 0, representing XX scattering, or p 0 < 0 and p 0 < 0, representing XX scattering. Let us first write the terms that arise from the middle line of (A1), apart from the factors of 2πδ(. . . ) Similarly the last line of (A1) contributes The collision term comes from Σ >,< by where unlike Tr above, tr denotes only the trace over Dirac matrices. Since we are interested in low densities, we can neglect terms of order F 3 , which means that we need only keep terms of order After carrying out the Dirac traces and combining like terms, we find that the respective contributions from the two diagrams are where dΠ p = d 4 p δ(p 2 − m 2 χ )/(2π) 3 . In the non-relativistic limit, it further simplifies since the squared matrix element in brackets is equal to 2 m 4 χ for C 1 , while for C 2 it depends on which of the theta functions are taken: [. . . ] = −m 4 χ for positive energies and +m 4 χ for negative energies. The resulting collision term is Here, we used the identities TrF p = TrF p = TrF p = Tr F p , as well as the fact that any terms with negative energies can be transformed to the corresponding phase space integrals with positive energy by changing p ↔ p . The next step is to make the ansatz where and n eq is the equilibrium number density. Then the momentum integrals can be carried out to get collision terms as a function of the density matrix n The integrals are all equal to (m χ T ) 9/2 /T times a dimensionless number, and there are only two different possibilities, depending upon whether the two energies in the Boltzmann factors are both initial/final state, or one initial and one final. We get where the two dimensionless integrals are and the zeroth component of the delta function is in terms of the non-relativistic dimensionless energies. The second forms of the integrals, in which the delta function of energies simplifies, are obtained by shifting p → p + p and k → k + p . We evaluated I s numerically. The divergent integral is inconsequential because it multipliesñn − 1 2 {n, n} ≡ 0, which vanishes identically. In retrospect we understand that this term is unphysical, since it corresponds to the interference of the t-and u-channel scattering diagrams, which vanishes for scattering of X withX. The relevant matrix evaluates to 4(ñ − n) Tr n −ñ 2 + n 2 = −6(n 11 + n 22 ) 0 n 21 so the collision term from scattering is which would appear in eq. (20) of ref. [21]. The normalization of σv s is chosen to agree with the usual definition, in which the low-energy cross section is σ ≈ g 4 m 2 χ /(4πm 4 v ), and the thermal averaging is done as in ref. [54].

Appendix B: Thermal decoherence in the Boltzmann equation
For the vector model, we have simulated the effect of thermal decoherence due to the oscillation rate depending on the momentum in the quantum Boltzmann equation for the density matrix where H is the Hubble rate and ω k = k 2 + m 2 χ . The k-dependence in the second term of the Hamiltonian implies that high-k parts of the distribution oscillate with slightly lower frequency than low-k parts, which is an additional source of decoherence that is neglected by integrating over momenta to reduce eq. (B1) to an equation for the number density matrix n. Our goal is to verify that this neglect is justified. For the scalar model, this issue is less important since decoherence is not a requirement for annihilations to occur.
To model the effect one would like to divide the particle distribution into several momentum bins. We will be content to take just two, labeled by s, l for small and large momenta, relative to the midpoint of the distribution. Accordingly, we split the density matrix into n t = n s + n l (B2) and one finds separate Boltzmann equations for each component, that are coupled to each other through the collision terms. The Boltzmann equations take the forṁ where σv s,a are the scattering and annihilation cross sections, the matrices S i , S, A i are defined as and we defined n ij = n s,ij + n l,ij , n t = n 11 + n 22 , and n i,t = n i,11 + n i, 22 . These expressions can be read from the form of the collision and annihilation terms in terms of the F matrices, before doing the final integral over the momentum k of the particle whose distribution is being tracked in the Boltzmann equation. If one adds the two equations together, they revert to the standard equation in terms of n t alone. The decoherence effect comes from the fact that the free Hamiltonians H s,l are slightly different for the two components, which for non-relativistic particles is The important feature is the difference between k 2 l and k 2 s , so for simplicity one could take, for example, k 2 s = 1 2 k 2 and k 2 l = 3 2 k 2 , which is a temperature-dependent split. For temperatures such that scattering is still in equilibrium, we can estimate k 2 ∼ 3m 2 χ /x, where x = m χ /T . After scatterings freeze out, the wavenumber redshifts as 1/a, so k 2 ∼ 3m 2 χ x f /x 2 . This effect can be important only in the early universe when the momenta are sufficiently large that k 2 /m 2 χ is not negligible. We have applied this formalism to check the early-universe solutions of section IV. We found no appreciable effect from this extra source of decoherence.

Appendix C: N-body simulation details
In this appendix, we describe the scattering and annihilation algorithms used in the simulations presented in this paper. Validation tests of the code are also presented and discussed.

C.1. Scattering
Elastic scattering between DM particles has been implemented stochastically on top of GADGET-2 in the same way as done by ref. [35], which was derived directly from the classical Boltzmann equation [5]. We summarize below the main relevant information and refer the reader to ref. [35] for a more detailed description.
The scattering rate for a DM particle of mass m χ at position r i and velocity v i , moving in an equal-mass particle background characterized by a normalized velocity distribution f v ( r, v) and a local density ρ( r), is Here σ/m χ is the DM scattering cross section per unit particle mass, which can be velocity-dependent. In simulations, individual physical particles cannot be resolved and all the properties of the latter should be translated to those of the simulation particles. For instance σ/m χ should be replaced by σ p /m p , where σ p and m p are the scattering cross section and the mass of the simulation particles, respectively. In a similar manner, the quantities f v ( r, v) and ρ( r) should be estimated from the volume within a sphere of radius h S , called scatter search radius, centered at the DM particle position r i . Assuming all the simulation particles within the scattering volume contribute equally, independently of their location (tophat kernel), eq. (C1) can be written as a sum over the N p neighboring particles [35] Hence the probability for particles i and j, separated by a distance smaller than h S , to scatter within the next time step of size ∆t, is given by wherer ≡ | r i − r j | and ρ ij is the target density, which is constant in this case because it is estimated using a top-hat kernel function. This is the simplest choice, but not the most common one used to implement DM self-interaction in N-body simulations. For instance, refs. [5,6,55] used a cubic spline kernel W (r, h S ) like the one already used in GADGET to compute the gravitational force in the context of smoothed-particle hydrodynamics (SPH) [56]. Such a kernel allows nearby particles separated by a distance less than h S to have higher scatter probability than those further apart because W is a smoothing function peaked atr = 0. For appropriate choices of h S , both the top-hat and cubic spline kernels provide results in agreement with analytical expectations [5,6,35,55] and with each other, within numerical uncertainties [35]. Therefore, we preferred using the simple and intuitive top-hat kernel with ρ ij given by eq. (C3) and with a fixed value of h S of the same order of the gravitational softening length . In particular, we consider h S = as chosen by ref. [35] because it gives the expected scattering rate, as discussed in section C.3. To see which particles do actually scatter at each time step, the probability in eq. (C3) is computed for each pair of nearby particles and compared with a random number drawn from a uniform distribution. For isotropic scattering and equal-mass particles, the post-scatter velocities are computed as v i,j = v CM ± (v rel /2)ê, where v i and v j are the initial velocities, v CM is the center-of-mass velocity between particle i and j, v rel is the magnitude of their relative velocity andê is a randomly oriented unit vector.
This scattering algorithm is very similar to that used in several SIDM cosmological simulations [4,6,[57][58][59][60][61][62], which mainly differ in the number of neighbors within the scattering volume. The majority of these simulations have treated the DM scattering as isotropic, which usually results from a short-range interaction mediated by a massive particle. In particular the mediator mass m V,φ should be much heavier than the DM particle momenta, which translates into m V,φ 10 −3 m χ for DM particles moving in Milky-Way-like galaxies today. If this is not the case, the cross-section will depend on the momentum exchange, which increases with the collision velocity or the scattering angle, leading typically to velocitydependent anisotropic scatterings. The latter are common in long-range interactions via light or massless mediators, which occur in several motivated DM scenarios such as mirror [63][64][65], atomic [66][67][68][69] and hidden sector DM models [70][71][72][73][74]. The simplest way to take them into account is to assume the scattering is still isotropic but the cross-section σ in eq. (C3) is replaced by the momentum transfer cross section σ T , defined as [75][76][77] where θ and dΩ are the scattering and solid angles, respectively. Here, the differential cross section dσ/dΩ can be derived within the Born approximation [78] from particle scattering mediated by a Yukawa interaction [2,79] where σ s is given by eq. (7) and ω = r m c with r m = m V /m χ or r m = m φ /m χ depending on the model under consideration. In the limit v rel ω, dσ/dΩ becomes velocity-independent, σ = (dσ/dΩ) dΩ ≈ σ s and the scattering is isotropic. In the opposite regime, eq. (C5) scales as ∝ v −4 rel like in Rutherford scattering which is mediated by the Coulomb potential. The approximation of using σ T instead of σ as the scattering cross section captures most of the effects of more complicated scattering dynamics where the DM velocity distribution is close to isotropic, which is the case of dwarf galaxies [77] 6 . This is because σ T estimates the average forward momentum lost during the collision, and several simulations used it to model DM long-range interactions [6,[82][83][84][85]. Ref. [80] proposed an alternative definition of σ T , namely to account for particle indistinguishability, which implies that dσ/dΩ is invariant under cos θ → − cos θ. The overall factor of 2 is used to give σT ≈ σ in the isotropic 6 Although in most astrophysical systems DM is expected to have an approximately isotropic velocity distribution, this is not the case of colliding galaxy clusters where there is a preferred direction along which DM particles collide [80,81].
regime, as done in ref. [86]. This new definition of the scattering cross section has been shown to provide better results than σ T in simulations of isolated DM halos [81]. An even better description of anisotropic scattering within the isotropic approximation seems to be given by the replacement of σ with the viscosity (or conductivity) cross section [77] because σ V takes into account of both forward and backward scatterings at the same time in addition to particle indistinguishability. Again, in order to match σ V ≈ σ in the isotropic regime, the original expression is multiplied by an overall factor of 3/2. In the present case, the DM particles are treated as being distinguishable since they have an associated oscillation phase ϕ, which might evolve differently during the simulation from particle to particle as discussed in section VI.
As noticed by ref. [2], σ V differs just by a O(1) factor from σ T and σT for distinguishable particles and apart from the rescaling factor, any of them can be taken as good measures of DM self-interactions, since systematic uncertainties in astrophysical observations are still too big to allow for discrimination between the different prescriptions. They become identical in the isotropic regime where v rel ω, which is well satisfied for the choice of model parameters and DM halos considered in this paper (see section V). We have chosen σ T because of its wide use in SIDM simulations, and we checked a posteriori that the use of the other prescriptions lead to the same simulation results for the DM models considered in this paper.
With eq. (C5), which is valid within the Born approximation for r α (recall α is the dark fine-structure constant), the momentum transfer cross section in eq. (C4) and its modified version in eq. (C6) turn out to be [87,88] and, analogously, the viscosity cross section in eq. (C7) becomes The tests for the code implementations of DM scattering are presented in section C.3.

C.2. Annihilation
We have implemented DM annihilation in a stochastic way similar to what was done for scattering. To the best of our knowledge, the first self-consistent implementation of DM annihilation in N-body simulations was performed in ref. [89], followed by ref. [90], where the energy released from annihilations to the surrounding gas particles was properly accounted for throughout the simulation. We followed a simplified version of this annihilation algorithm, without including the energy transfer between different particle species, because our simulation considers only DM particles χ and the annihilation products are assumed to decay quickly into lighter states that can escape galaxies without affecting their environment and therefore the observable quantities.
In detail, the annihilation rate for a DM particle of mass m χ at position r i , moving in an equal-mass particle background of density ρ( r), is where σ ann v is the velocity-averaged annihilation cross section. Depending on the particle nature of DM, the estimation of ρ( r) can include all the DM particles in the system (for Majorana fermions or real scalars) or just antiparticles (for Dirac fermions or complex scalars).
Following the same steps as in section C.1 and introducing the annihilation search radius h A , we can generalize eq. (C10) and write the probability for DM particles i and j, separated by a distance smaller than h A , to annihilate within the next time step of size ∆t as wherer ≡ | r i − r j | and ρ ij is the target density, which is estimated using a top-hat kernel function. The velocityaveraged annihilation cross section σ ann v ij depends generally on the relative velocity v rel = | v i − v j | between particles i and j. In the non-relativistic limit, valid in the context of galaxies, it is usually expanded in powers of v rel as where σ ann, s and σ ann, p are constants corresponding to the s-wave and the p-wave annihilation terms, respectively. For the models considered in this paper, only σ ann, s contributes to σ ann v ij , which is given by what we refer to as σv a in eq. (5). This makes the annihilation probability in eq. (C11) velocity-independent. As in the scattering case, the annihilation search radius h A entering the density ρ ij is in principle a free parameter that should be chosen to reproduce analytical results. As we will discuss in section C.3, h A = turns out to be the best choice (recall is the gravitational softening length). To see which particles annihilate at each time step, the probability in eq. (C11) is computed for each pair of nearby simulation particles and compared to a random number. If the latter is below than the former, annihilation happens and the two particles in the event are removed from the system.

C.3. Validation tests
The simplest test for both our scattering and annihilation algorithms is a uniform cube of N c particles moving through a background of stationary particles with constant number density n b . All the particles making up the simulated system have the same mass m p . To allow simple predictions, we impose that the cube particles move with constant speed v 0 along the same axis, they can scatter with constant cross section σ p at most once, and gravity is turned off [35].
The scattering or annihilation rate for each simulation particle in the cube is Γ = n b m p (σv 0 /m χ ), where m χ and σ are respectively the mass and the cross section of the physical particles, whose properties have been translated into those of the simulation particles via σ p = m p (σ/m χ ). Naively, we can estimate the expected number of interactions after a time t as This expression has however several limitations because not only it does not take into account that not all the particles in the cube have interacted between zero and time t but also that those that have scattered or annihilated could do it just once. These effects are properly captured by considering how the number of cube particles changes with time, which can be described by With this in mind, the expected number of scattering or annihilating particles in the simulation after a time t can be better expressed by where N c (0) is the initial number of particles in the cube. The comparison between N exp and the number N of cube particles that have scattered (annihilated) in our test simulations is plotted in Fig. 9 as a function of the scatter (annihilation) search radius h. Here, each point corresponds to a simulation where only scattering was turned on, while stars are used for simulations in which only annihilation was active. As already noticed in refs. [5,35], N falls below than expected only for h smaller than 20% of the mean background interparticle separation, but this minimum value depends heavily on the chosen time step ∆t [35], as confirmed by Fig. 9. The h-dependence of the number N of scattering or annihilating particles in simulations arises in the "probabilitysaturated" regime, where the probability for a particle to scatter or annihilate within a time step, given by eqs. (C3) or (C11), becomes greater than unity. The latter probability enters the definition of N , which depends on it (∝ h −3 ) and on the number of neighbouring particles that a particle finds at each time step (∝ h 3 ). This makes N generally insensitive to h, except in the probability-saturated regime, where N ∝ h 3 as shown by the solid and dashed lines in the right panel of Fig. 9. In order to avoid probability saturation, a shorter time step should be used when using smaller h, since the probability for a pair of particles to scatter or annihilate is proportional to ∆t/h 3 .
Although ref. [35] suggested that probabilities exceeding unity within each time step should not appear in SIDM simulations powered by GADGET for reasonable values of σ/m χ , we decided to implement a time-step criterion. In particular, similar to what was done in ref. [60], an individual particle time step ∆t is modified by rearranging eq. (C3) as if the probability of interaction for any pair involving such a particle was greater than P max = 0.1 during the last tree-walk. This restriction is important only for scattering particles, because annihilating particles are removed from the system as soon as the interaction takes place and therefore their time step becomes meaningless. Although limiting the individual time step makes the Monte Carlo method more computationally costly, it allows for suppression of the probability of having multiple scatterings in the same ∆t, which is an important issue in SIDM simulations.
For simulations with only scattering, the directions and velocities of the scattered particles can also be compared to the expected normalized distributions. The latter can be obtained by transforming the differential cross section from the centre-of-mass frame of the collision to the simulation frame. For isotropic scatterings, these distributions turn out to be the same for both background and cube particles and take a simple form, with f (θ) = sin 2θ, f (φ) = (2π) −1 and f (v) = 2v, where the latter is valid for v ≤ v 0 , and becomes zero otherwise. Fig. 10 shows the comparison between the expected distributions and those reconstructed from one of our test simulations, which agree to within 1σ uncertainty.
Another important test of our code is to check whether particle scattering and annihilation are well modeled in isolated DM halos. We focus on halos with Hernquist profile [33], whose density distribution can be written as where ρ H ≡ M/(2πa 3 ), M is the total halo mass and a the scale radius. We have chosen this type of halo because it has a finite mass without need of truncation and the phase-space distribution function f (E), with E being the particle energy, has an analytic form [33]. The latter property allows to easily generate equilibrium initial conditions for N-body simulations and compute scattering or annihilation quantities analytically. The initial conditions for particle positions and velocities making up our test halos have been generated randomly from the Hernquist distribution function f (E) using the von Neumann rejection method [91], as originally done in ref. [92]. To prevent centroid motion of the generated DM halo during the simulation, we set its initial centre-of-mass position and velocity to zero by an overall boost. Since we consider a large range of halo masses in this paper, ranging from 10 10 to 10 15 M , we tested our code with simulated Hernquist halos having different values of M and a, finding good agreement between the analytical expectations and the simulation outcomes for all of them. As a prototype example, we focus here just on an isolated Hernquist halo with total mass M = 10 14 M and scale radius a = 225 kpc. The simulations for such a halo were run with N = 128 3 particles, each having mass m p 4.8 × 10 8 M , and the Plummer-equivalent gravitational softening length was set to = 4.4 kpc. We evolved the generated halos with collisionless DM first to study their stability and to check whether cores can form as numerical artifacts. The results for our prototype halo as a function of the simulation time are shown in the top left plot of Fig. 11. For a suitable choice of time-integration and tree-force accuracy parameters, η = 0.005 and α = 0.0012 respectively, the density and velocity distributions remained unchanged except for the formation of a small constant core with size similar to the gravitational softening length , as observed by ref. [35].
When DM scattering is turned on, these cores quickly become larger because particles scatter mostly in high density regions until the cores settle to a size that is independent of the value of σ/m χ [57]. The right panel of Fig. 11 shows the results for the density profile of our prototype halo at different simulation times for a constant scattering cross section per unit DM mass of 1.0 cm 2 /g. The resulting profiles are well fitted by a cored-Hernquist profile of the form [81] where x ≡ r/a, r c is the core-radius and β is an index controlling the sharpness of the transition from ρ ∝ x −1 to a constant density. Leaving r c and β as free parameters in the fit, we found that their best-fit values across a time evolution of 10 Gyr are r c /a ∈ (0.06, 0.14) and β ∈ (3.5, 62.0), in agreement with refs. [57,81]. The results shown in the right panel of Fig. 11 have been obtained by choosing the scatter search radius equal to the gravitational softening length in the computation The solid lines correspond to the best-fit cored-Hernquist profiles, given by eq. (C18), where rc and β are left as free parameters. The 1σ error bar for each data point is computed assuming that the number of particles in each bin is Poisson distributed. Bottom: Extracted scattering rate per particle for the same Hernquist profile DM halo after 3 Gyr, for different choices of the scatter search radius hS, and its comparison with the theoretical expectation. Although the simulations were run with σ/mχ = 1.0 cm 2 /g, the result is independent of the scattering cross section since a ratio is considered. The colored crosses along the analytical curve correspond to the radius equal to hS. The 1σ uncertainty for each colored line is computed assuming that the number of particles in each bin is Poisson distributed and displayed with a same-color shaded region. In all three panels, the gray dot-dashed vertical line shows the position of the gravitational softening length used in all simulations for the considered halo. We used a time-integration parameter of η = 0.005 and tree-force accuracy of α = 0.0012 in each simulation run.
of the scattering probability, given by eq. (C3). This was done in light of the result shown in the bottom panel of Fig. 11, where the scattering rate extracted from simulations is compared to the analytical prediction. In particular, the latter can be computed by integrating eq. (C1) over the velocity distribution function, obtaining the average rate for particles at position r [35] Γ scatt ( r) = σ v pair ( r) ρ( r) m χ , where ρ( r) is given by eq. (C17) and σ v pair = σ v pair if the scattering cross section is velocity-independent. The mean pairwise particle velocity v pair can be computed from the one-dimensional velocity dispersion σ 1D for Hernquist halos [33] as v pair = (4/ √ π) σ 1D , where we have assumed the velocities are isotropic and follow a Maxwell-Boltzmann distribution. From simulations, the scattering rate per particle as a function of the radius r can be estimated as the ratio between the number of particles that have scattered within the radial bin including r, and the time averaged number of particles within the same bin [35]. To avoid any modification of the density profile and velocity distribution due to DM self-interactions that would lead the scattering rate not to follow the analytic prediction, we turned off the change in the momenta of the scattered particles in the simulations used in the bottom panel of Fig. 11.
Such a plot clearly shows that our code accurately reproduces the scattering rate within the halo except for distances less than h S , where the extracted rate falls below the analytic result. This can be explained by the fact that the scatter search radius acts as a scale below which the particle density entering the scattering probability becomes smooth, preventing a faithful reconstruction of the real density. Although it suggests choosing a small value of h S in order to correctly capture the scattering dynamics in small high-density regions, there is a natural lower bound for h S set by the gravitational softening length . As we have found in simulation runs for collisionless DM, the density gets affected by a small core of size of the order of the softening length, which arises from smoothing the gravitational potential at r < and thus the particle distribution at these scales. This explains why Γ scatt ceases to change for scales smaller than in simulations where h S < , as displayed in the bottom panel of Fig. 11. Therefore, reducing h S to values below cannot improve the agreement between the extracted rate and the analytic result, but rather it tends to create probability saturation problems that can be solved by choosing smaller time steps or a time-step delimiter, as already discussed above. Considering all these factors and the results in the bottom panel of Fig. 11, we find that setting h S = provides the best reconstruction of scattering dynamics at low scales, and it limits the occurrence of events where the scattering probability becomes greater than unity.
With the same Hernquist halos it is possible to also test the goodness of our annihilation algorithm. For instance, the annihilation rate extracted from simulations can be compared to the analytic result, which is given by eq. (C19) with the replacement of σ v pair with σ ann v . The extracted annihilation rate as a function of the radius r has been obtained in the same way as done for scattering events, namely as the ratio between the number of annihilated particles located in the radial bin including r and the time averaged number of particles within the same bin. To allow for a direct comparison to the analytical result we did not remove the annihilated particles from the system, but the annihilation algorithm was used to count and localize the particles that actually annihilate. The results from this comparison are shown in the left panel of Fig. 12, where we have chosen a constant (s-wave) velocity-averaged annihilation cross section per unit DM mass of 100 cm 2 /g km/s. Similarly to what was observed in the scattering case, choosing an annihilation search radius h A equal to the gravitational softening length provides the best reconstruction of the annihilation rate. This is because the particle distribution at scales below h A and are smoothed, leading to a decreased annihilation rate compared to the true unsmoothed one.
The effect that DM annihilation has on the time evolution of the halo density profile was first studied by ref. [17], which proposed the following analytic formula for ρ(r, t) d dt where ρ A ≡ m χ /( σ ann v t 0 ) and t 0 is the age of the universe today. This equation comes directly from the definition of the DM annihilation rate and can be easily adapted to our simulation setup, where the annihilation probability is given by eq. (C11), by replacing ρ A with ρ A /2. The factor 1/2 arises because the halo density is reduced by two units of the simulation particle mass in each annihilation event. The general solution of eq. (C21) is given by ρ(r, t) = 1 ρ core (t) where the latter definition is valid for our simulations and t ini is the initial time, which can be taken as the time when the simulation starts and the density profile is ρ(r, t ini ). One observes that the core density ρ core (t) is generally greater than ρ A and the equality is reached when (t − t ini ) = t 0 . Focusing on just DM halos with initial Hernquist profile and taking t ini = 0, we find that the density profile of halos undergoing DM annihilation should be described after a time t by ρ(x, t) = ρ H x (1 + x) 3 + ρ H /ρ core (t) (C23) with x ≡ r/a, which has been obtained by substituting eq. (C17) into eq. (C22). The halo density is therefore characterized by a core of constant density ρ core (t) at small scales. The bottom panel of Fig. 12 shows the evolution of the density profile of the prototype Hernquist halo in the simulation at different times t and its comparison with the theoretical prediction. Although the data points and the solid lines given by eq. (C23) match very well at large radii, there is some disagreement at small distances. These differences are not eliminated by changing the annihilation search radius, as shown by the right panel of Fig. 12, where the analytical prediction is compared to the simulation outcome at t = 10 Gyr for different values of h A . A graphical comparison between the analytical curve and the colored lines  Fig. 11, but for the extracted annihilation rate per particle and different choices of the annihilation search radius hA. The simulations for the prototype Hernquist halo were run with σannv /mχ = 100 cm 2 /g km/s. Top right: Radial density profile of the prototype Hernquist halo undergoing DM annihilation for 10 Gyr, for different choices of hA. The black dashed line displays the original density profile in eq. (C17). The violet solid line corresponds to the theoretical prediction given by eq. (C23) at (t − tini) = 10 Gyr. Bottom: Similar to the right panel of Fig. 11, but for particle annihilation with the same constant velocity-averaged cross section considered above. Here we chose the annihilation search radius as hA = . Each colored solid curve shows the analytical expectation given by eq. (C23) at the time of the corresponding same-color data points. The dotted lines represent the best-fit cored-Hernquist profiles, given by eq. (C18), where rc and β are left as free parameters.
suggests that eq. (C23) does not provide a good fit to the simulation data at low scales, independently of the choice of h A , although it accurately captures the physics at large radii. We confirmed this observation by fitting the simulation data with eq. (C23), in which ρ core (t) was left as a free parameter, leading to agreement within 1σ between the best fit value and that computed by eq. (C22) at different times t, independently of the value of h A . In analogy to the case of scattering, as shown in the top right panel of Fig. 11, we investigated the cored-Hernquist profile function given by eq. (C18) for fitting the annihilation data displayed in the bottom panel of Fig. 12, leaving r c and β as free parameters. The results of this fit are shown with colored dotted lines in the same panel, with best-fit values r c /a ∈ (0.06, 0.48) and β ∈ (1.2, 2.0). They do not give significant improvement with respect to eq. (C23).