New Dissipation Mechanisms from Multi-level Dark Matter Scattering

Multi-level dark matter with diagonal and off-diagonal interactions shows a rich phenomenology in its self-scattering. If the interactions are mediated by a particle that is less massive than the dark matter, Sommerfeld effect can lead to resonant enhancement of the scattering. For mediators lighter than the level separation, dark matter particles can upscatter to excited states and de-excite by emitting these mediators. We compute these cross-sections, both above and below the kinematic threshold, in a generic two-component dark matter model and identify the large inelastic cross-section as a result of maximal mixing between the two states. A new route for cooling of large dark matter halos and a new drag force between two colliding halos are identified and shown to arise purely from the inelastic scattering.


I. INTRODUCTION
Galaxies are observed to be surrounded by more massive halo-like structures made of a substance whose particle nature still remains unknown. The formation of these halos is predicted within the paradigm that dark matter (DM) is made of cold collisionless particles. Agreement with several other cosmological and astrophysical observations at widely different length scales has given strong support for this simple paradigm. Notwithstanding this success, it fails to provide explanation for some discrepancies between the predictions and observations of the shapes and abundances of these halos at sub-galactic length scales, viz., the small scale problems [1,2].
One of these problems is often referred to as the corecusp problem [3][4][5][6]. The central density profile of dwarf galaxies is observed to be cored, while simulations with the standard cold collisionless DM typically lead to a denser cuspy profile (∼ 1/r) near the centre. Complex astrophysical processes involving baryonic matter, e.g., tidal effects and supernova explosions remove matter from the central region, may lead to cored profiles reducing the discrepancy [7][8][9][10][11][12][13]. Indeed, these feedback processes could very well be the missing ingredient in the simulations. However, it is challenging to accurately model these processes and it is not yet established if one can obtain sufficient feedback in realistic models [14][15][16][17]. As such, this problem remains open and motivates us to consider other possibilities as well.
Self-scattering of DM particles have been shown to be effective in solving the core-cusp problem [18,19]. In this class of models, known as self-interacting DM (SIDM), the DM particles have strong interactions with each other. During the non-linear phase of structure formation, when the central density of a halo becomes large, * anirbandas@theory.tifr.res.in † bdasgupta@theory.tifr.res.in the self-scattering generates outward pressure. When this pressure equates the gravitational pull on the matter, any further accumulation of DM at the centre ceases and a stable core is formed. While these interactions may be obtained with massive mediators, the cross section in such models is velocity-independent and is strongly constrained from various observations [20]. Simulations indicate that a hard-sphere scattering cross-section per unit DM mass σ/M ∼ 0.5 − 5 cm 2 /g is required to form the cores at the centres of galaxies with DM velocity v rms = 30 − 100 km/s [21][22][23]. On the other hand, objects with larger v rms put stronger bounds on these same cross-sections [24][25][26][27][28]. For example, the non-observation of drag force between the DM components of two merging clusters puts an upper bound σ/M 1 cm 2 /g [29,30]. But a tighter constraint σ/M 0.1 cm 2 /g is given by the stellar kinematics and weak lensing data in galaxy clusters [23]. As self-interaction helps virialize DM halos efficiently, they become rounder compared to the anisotropic halos predicted by collisionless DM scenario. Therefore the observed triaxialities of halos provides an upper bound to the strength of self-interaction among the DM particles. A comparison of observations with Nbody simulation shows that a scattering cross-section as large as σ/M 1 cm 2 /g could be allowed at the larger velocity end, i.e., v ∼ 1000 km/s [31].
Therefore, one typically considers lighter mediators that lead to velocity-dependent cross-sections [32,33]. In these models, a velocity-dependent cross-section is obtained from a long-range interaction between DM particles, and have been shown to be able to satisfy observational constraints [31]. The long ranged interactions, on the other hand, have interesting phenomenology [34,35]. When the mass of the mediator particle is smaller than that of the DM, Sommerfeld effect may cause resonant scattering. In this regime, the cross-section is resonantly enhanced through virtual bound state formation [36].
A multi-component DM system brings new features in the scattering phenomenology and the dynamics of halo formation. In particular, it can solve the core-arXiv:1709.06577v2 [hep-ph] 5 Feb 2018 cusp problem by heat flow from the hotter outer region to the colder inner core. The possibility of elastic as well as inelastic scattering depending on the energy of the particles, gives a rich phenomenology that has not been explored fully. For example, the excitation and deexcitation of DM particles can give rise to observable indirect detection signals. Also, the energy dissipation from the inelastic scattering, followed by de-excitation, might lead to significant change in the shape and density profile of DM halos. Some of these features have been discussed in the context of atomic DM models in Ref. [37][38][39], two-level DM systems with purely off-diagonal interaction in Ref. [40,41] and dark bremsstrahlung process [42]. Multi-level DM models also have interesting phenomenology in the context of direct and indirect detection experiments [43][44][45].
In this work, we take a two-level SIDM model with light particles mediating both diagonal and off-diagonal interactions. In this model, the DM particles can not only elastically scatter due to the diagonal interactions, but they may also get excited to the more massive partner of DM due to the off-diagonal interactions and subsequently de-excite by emitting the light mediator particles. This leads to additional dissipation. We compute these cross sections, analytically explain their behavior in various regimes, and study the cooling of DM halos due to the new dissipation mechanism. We further identify a new dissipation-induced drag force between two colliding halos in such models.
The paper is organized as follows. In Sec. II, we discuss a formalism for two-level scattering with a description of the minimal SIDM model considered here, numerically compute the elastic and inelastic cross-sections, and explain their key features using simple analytical estimates. We then outline the key signatures and possible constraints in Sec. III, and, in Sec. IV, conclude with a summary of our results and avenues for future work.

II. TWO LEVEL DM & SCATTERING
To discuss the phenomenology of a multi-level DM model, we concentrate on a simple two-level DM system with a small mass gap ∆ between the two states, χ 1 and χ 2 , with masses M and M + ∆, respectively. We further assume a dark Z 2 symmetry under which χ 1,2 have charges ∓1, respectively. Two lighter scalars ρ 1 and ρ 2 with charges ±1, couple to the DM states as For simplicity, we have assumed the masses and couplings of ρ 1,2 to be the same and equal to m ρ and f , respectively. Most of our discussion is insensitive to these simplifying assumptions, and we will outline how the results would change in a more general model, wherever necessary.
Two colliding DM particles that are initially in the ground state can either stay in the ground state (elastic) or upscatter to the excited state (inelastic). For up- FIG. 1. Typical Feynman diagrams for elastic self-scattering of DM in the ground state (left) and for upscattering induced decay (right). The intermediate vertical lines represent multiple exchanges of scalar ρ particles in the nonrelativistic limit of the incoming DM particles.
scattering to occur, the incoming particles need to have enough kinetic energy to overcome the mass gap 2∆ between the two 2-body states. In addition, even if there is not enough kinetic energy, the excited state can be produced as virtual particles in the intermediate steps of the collisions. Moreover, the scattering cross-section between nonrelativistic DM particles is enhanced due to multiple exchanges of the light ρ particle, influencing the Sommerfeld effect. Schematic Feynman diagrams for the possible elastic and inelastic scatterings are shown in Fig. 1 where the vertical lines represent exchange of many ρ 1,2 particles in the nonrelativistic regime. In the case of inelastic scattering, the final particles decay to the ground state by emitting two light particles. This process is essentially an energy loss process, and is expected to have interesting phenomenology.
The scattering cross-sections are computed by calculating the transition amplitude between an allowed initial state |i and final state |f . The possible 2-body states are |χ 1 χ 1 , |χ 2 χ 2 and |χ 1 χ 2 . However, it is easy to see from Eq.(1) that |χ 1 χ 1 and |χ 2 χ 2 are decoupled from |χ 1 χ 2 , due to the Z 2 symmetry. Therefore, assuming that the DM particles are initially in the ground state, it suffices to work in a Hilbert space spanned by |χ 1 χ 1 and |χ 2 χ 2 only. One can remove this restriction, but it makes the calculation more difficult without yielding any qualitatively new features. This is our motivation for using two oppositely charged scalars, as opposed to a single scalar that one may naively think to be the simpler case. We neglect the scattering between χ 1 and χ 2 , because χ 2 decays soon after freeze-out and its abundance is rapidly depleted. For the same reason, the scattering between two χ 2 particles is negligible. Therefore, we have two channels, i.e., |χ 1 χ 1 → |χ 1 χ 1 (elastic) and |χ 1 χ 1 → |χ 2 χ 2 (inelastic).
The overlap between two 2-body states is defined as Ψ ab ≡ χ a χ a |χ b χ b with a, b = 1, 2 and satisfies the Schrödinger equation, where is the orbital angular momentum, and k and µ are two diagonal matrices with the momentum and reduced mass of the incoming 2-body state, respectively, k = k a δ ab , and µ = µ a δ ab . ( The incoming momentum k a is different for the two 2body states due to the presence of the mass gap 2∆ = 2(M 2 − M 1 ) between the states |χ 1 χ 1 and |χ 2 χ 2 . Depending on the energy of the incoming particles, two cases are possible: • Below threshold, µ 1 v 2 /2 < 2∆, when the initial energy of the incoming state is below threshold, then the heavier |χ 2 χ 2 state is kinematically closed as χ 2 s cannot be produced onshell, • Above threshold, µ 1 v 2 /2 > 2∆, when the incoming energy is above threshold, the excited state is open and DM particles can upscatter to the excited state.
The exchange of the scalar particles between the DM particles as dictated by the Lagrangian in Eq.(1) gives rise to an attractive potential between the 2-body states in the nonrelativistic limit of the theory. The potential matrix V (r) in Eq.(2) is given by where α ≡ f 2 /4π . This attractive Yukawa potential matrix, with equal entries, is a result of assuming the same interaction strength between either pair of 2-body states. The structure of the potential matrix would be different in other DM models, e.g., a broken dark gauge symmetry would provide both attractive and repulsive interactions, which will become purely repulsive for a late-time asymmetric DM population. With additional scalars one can engineer diagonal and off-diagonal potentials of different strengths. As we shall see, the qualitative nature of many results discussed in this work remain unchanged as long as the matrix has nonvanishing off-diagonal components. Therefore we shall not delve into these details here.
The set of Schrödinger equations in Eq.(2) is to be solved with appropriate boundary conditions for above and below threshold scatterings. The equations are solved for each partial wave and the large-r wavefunctions are matched with plane wave solutions to extract the elements of the transition matrix T which consists of the transition amplitudes from state |i to |f . The scattering matrix S is written as S ≡ 1 − T . Finally the scattering matrices are added upto some = max , to yield the total cross-section. Although theoretically max goes upto infinity, in practice a finite value must be chosen by ensuring numerical convergence of the sum. This value depends on the range of the potential and the momentum of the incoming particles.
The total scattering cross-section σ tot is given by This definition gives equal weight to all scattering angles which is useful in case of hard sphere scattering mediated by a heavy particle with short range. For a small mediator mass, the cross-section is peaked in the forward and backward directions. While this leads to an overall large value of cross-section, the effective momentum transfer in each collision is small if the particles are identical. However, in a DM halo and in N-body computer simulations, momentum transfer during a collision is the quantity that determines the virialization and shape of a halo during its evolution and the dynamics of colliding halos. In Ref. [46], it was pointed out that the transfer crosssection σ T , which removes the forward direction peak, is a more important quantity for transport phenomena. It is defined as Another quantity which is often used is the viscosity cross-section σ V [46], This definition removes the contributions of the backward scatterings, in addition. The details of the numerical calculations of these cross sections are given in appendix A. We now discuss the nature of the elastic and inelastic cross-sections in the limits when µ 1 v 2 /2 < 2∆, i.e., below threshold and when µ 1 v 2 /2 > 2∆, i.e., above threshold.

A. Below Threshold
If the total energy of the incoming particles is below the threshold, i.e., µ 1 v 2 /2 < 2∆, then the excited state cannot be produced onshell as a final state, and only the elastic channel is relevant. The elastic cross-section shown in Fig. 2 can be understood based on a simple analytical argument. In the regime µ 1 v 2 < 2∆, transitions to the heavier excited state are classically forbidden by the vanishing tunnelling probability and can be neglected. As a result, it is easiest to diagonalize Eq.(2) locally at at each value of r, and solve only for the elastic scattering cross-section in the ground state with the potentialṼ 1 (r) given byṼ The χ1χ1 → χ1χ1 cross-section in the regime µ1v 2 /2 < 2∆, i.e., below threshold for our two-level SIDM model. This is approximately the same as in a single state SIDM model but with the potential given by Eq.(9).
Since the potentials are functions of the radial distance, the rotation angle will also be a function of r, If the spatial derivatives of θ(r) are small enough, then one can assume the system to be adiabatic [47], i.e., the two instantaneous energy levels do not mix. Note that θ(r) → π/4 for r → 0 and θ(r) → 0 for r → ∞. Therefore the two states are completely unmixed at large distance but get maximally mixed as the incoming particles get closer. The radial dependence of θ(r) is shown in the green curve in the left panel of Fig. 3. The transition from zero to π/4 happens in the region where |V 1 (r)| ∆. Beyond this point towards large r, mixing ceases as the potential in the off-diagonal position becomes smaller than the diagonal mass gap term. This behaviour of the mixing angle allows us to use the rotated basis to determine the elastic scattering cross-section in the below threshold regime. The smallness of the radial dependence of the elements of the rotation matrix ensures that the system remains in one of the eigenstates during the complete scattering. In Fig. 3, we also show the behaviour of the eigenvalues of the potential matrix with r for = 0. The eigenvalues never cross each other and their separation goes as ∼ 2/r for r → 0 reaching a constant 2∆ for large r. This remains true for all higher partial waves and signals the fact that the evolution is always adiabatic and the elastic cross section is approximately given by the effective potential in Eq.(9).

B. Above Threshold
If the energy of the incoming particles is sufficiently large, i.e., µ 1 v 2 /2 > 2∆, then the excited state can be produced onshell. In this case, the DM particles in the ground state can upscatter inelastically to the excited state. A measure of the inelasticity in the system is given by the magnitude of the off-diagonal elements of the scattering matrix S . We study the behaviour of the scattering matrix using the variable phase method following Ref. [48]. In this method the wavefunction is first written in an integral form, with the Riccati-Bessel functions J (kr) and N (kr) as defined in Eq.(A4). In order to isolate the part of the wavefunction arising from the interaction potential V (r), it helps to define another function F (r) as On substitution of this in Eq.(11), one gets where one can clearly identify H (1),(2) ≡ (µ/k) 1/2 (N ± iJ ) as the free wave solutions and the unknown functions F and F * as the distortions to the plane-wave solution due to the potential. The scattering matrix function The significance of the function S (r) lies in the fact that its asymptotic value at large r yields the scattering matrix S . The differential equation for S is easily obtained, by taking a derivative of the previous equation and using Eq. (11,12), The initial condition is S (0) = 1 , because at r → 0 the function S (r), as given by Eqs. (14) and (12), has zero off-diagonal entries because the integral in Eq. (12) vanishes.  9). The mixing angle θ(r) is shown as the green dotted line. The unit on the vertical axis is arbitrary. Note that the eigenvalues remain well-separated and there is no level-crossing, explaining why the evolution is adiabatic. The parameters chosen here correspond to an above threshold scenario, but the behaviour of the eigenvalues does not depend on that. (Right) Real and imaginary parts of the off-diagonal component of the scattering matrix function, Re(S12(r)) (blue solid) and Im(S12(r)) (red dashed), as well as the instantaneous mixing angle θ(r) (green dotted). Note that θ(r) reaches π/4 at r → 0, shown as a grey dot-dashed line. The real and imaginary parts of S12(r) are shown multiplied by factors of 100 and 10, respectively, for visual clarity. Note that the off-diagonal component S12(r) deviates from its value at large r only at r 50, where the eigenvalues are well-separated but the mixing is maximal. This shows that any nonzero inelastic scattering, that comes from a nonzero S12(r), is a result of maximal mixing and not level-mixing or lack of adiabaticity. The velocity was chosen to be 100 km/s (above threshold).
The off-diagonal components of S (r) track the behaviour of the inelastic cross-section with r. In the right panel of Fig. 3, we plot the real and the imaginary parts of the off-diagonal elements of S as a function of r, along with the rotation angle θ(r) of the potential matrix V (r) in Eq.(10). The inelastic cross-section grows from zero to its asymptotic value during the region where the mixing angle is π/4 (maximal mixing). More precisely, it saturates at around r ∼ 1/m ρ (In this case, m ρ = 0.1 GeV). Beyond this point the off-diagonal potentials are exponentially screened and the off-diagonal coupling between two states vanishes. Also note that nothing special happens in the 'nonadiabatic' region, i.e., where the angle varies from π/4 to zero towards large r. While we show this for particular values of the parameters in the potential, this qualitative behaviour does not change for other values of the parameters. Therefore, one can conclude that the inelasticity is driven by the maximal mixing between two states near the origin (the adiabatic mixing) which yields the large upscattering cross-section from the ground state, and not by the nonadiabaticity in the system. As soon as this mixing goes to zero, the inelastic cross-section saturates to its asymptotic value. We also note that the asymptotic value of the off-diagonal elements of the S matrix is significantly large, which hints towards a large inelastic cross-section in the presence of an off-diagonal potential. For obtaining our numerical results, we have used the method shown in Appendix A, but the main advantage of the variable phase method described here is that it reveals that the origin of inelasticity is large mixing, not non-adiabaticity. A variant of this method was previously used to compute the crosssections in Ref. [41].
The left panel of Fig. 4 shows the velocity-dependence of the elastic cross-section for a few representative values of the mediator mass m ρ , as indicated in the figure.
In general, we see that irrespective of the value of m ρ , the σ el is larger for small DM velocity and decreases for large velocity. Therefore it is possible to enhance the self-scattering in the dwarf-sized objects and address the core-cusp problem, while simultaneously suppressing it in the larger cluster-size objects and satisfying the upper bounds coming from colliding clusters [29,30]. The values of m ρ were chosen such that they span across a resonance. The curve labelled by m ρ = 0.116 GeV corresponds to a resonance in the cross-section and hence shows large enhancement in the small velocity regime. On either side of this resonance, the cross-section decreases. These features are unaltered relative to singlelevel DM models.
In the right panel of Fig. 4, we show the behaviour of the elastic transfer and viscosity cross-sections with the mediator mass in the ∆ → 0 limit. Two distinct regions, Born (αM/m ρ 1) and resonant (αM/m ρ 1), are apparent. The dashed grey line shows an approximate analytical estimate of the cross-section in the Born limit [32]. Although the physical system in this work is different than that in single-state DM models, it is possible to get an approximate expression of the Born crosssection by a substitution α → 2α. This substitution is based on the understanding that in the limit of small ∆, the two states |χ 1 χ 1 and |χ 2 χ 2 become indistinguishable and related to each other through the relation |χ 2 χ 2 = (−1) +s |χ 1 χ 1 where s is the total spin of the state, as was explained in a previous paper by us [45]. Only one linear combination of the states survives in this limit and one can map the two-level system onto a single state with an effective potential This also explains the substitution α → 2α in our adapted estimate of the resonant condition given by [49], that explains the positions of the resonances and agrees very well with the numerical results. In all these cases the inelastic cross-section is almost equal to the elastic cross-section as shown in the right panel of Fig. 4. This strong correlation between the two cross-sections is a result of setting all components of the potential matrix in Eq.(5) to be equal. If the diagonal potentials are weaker than the off-diagonal counterparts, then the two cross-sections can be different by several orders of magnitude. An extreme example of such a case is the model of a two-component Majorana DM particles charged under a broken U(1) gauge symmetry [40,41]. The conserved currents in this model are given byχ 1 γ µ χ 2 andχ 2 γ µ χ 1 . Hence at tree level, elastic scattering χ 1,2 χ 1,2 → χ 1,2 χ 1,2 is not possible although inelastic scattering χ 1,2 χ 1,2 → χ 2,1 χ 2,1 can take place through an exchange of a gauge particle. The lowest order elastic process will involve at least one loop. Therefore in the Born limit (large m ρ ), the elastic cross-section will be suppressed. In the resonant and classical regime these interactions would give rise to purely off-diagonal attractive Yukawa potentials and both elastic and inelastic cross-sections will be comparable (see Fig. 5). However, depending on the parameter values, one may be larger than the other.

III. SIGNATURES OF INELASTIC SCATTERING
The new effects of this dissipation mechanism is driven by the large inelastic scattering rate given by Γ up ≡ n χ σ in v. An order of magnitude estimate of the timescale associated with the upscattering rate is given by (18) The DM velocity was chosen to be O(1000) km/s so that the upscattering and decay processes are kinematically allowed. Clearly this typical timescale is 1-2 orders of magnitude larger than the age of the Universe, whereas large DM densities required for upscattering to take place have only been present for a much shorter time (only since nonlinear structures have formed). Therefore the effects of these upscatterings cannot be too large. We now discuss two possible effects due to this inelastic scattering.

A. Halo Cooling
If the upscattering rate n χ σ in v is not too small, χ 1 can upscatter to the excited state χ 2 and thus produced χ 2 will promptly decay into the light mediator particle ρ and χ 1 . If the χ 1 − ρ scattering cross-section is small in the given DM halo, then these light particles may escape the halo, thereby cooling the halo at a considerable rate. Large upscattering requires the colliding DM particles to be energetic enough so that sufficient phase space is available for the excited state. For example, a 10 GeV mass DM with ∆ = 1 MeV has a velocity threshold of ∼ 1000 km/s. Thus, this phenomenon will mostly be important in objects with large DM velocity dispersions, e.g., in large galaxies and galaxy clusters.
A thorough analysis of the effect of this dissipation mechanism on the DM halo structure does not yet exist in the literature. We will not attempt a full treatment here. However, a qualitative understanding can be gained from the response of DM halos for similar cooling processes present in the baryonic matter, as we recap below. After falling towards the centre of a halo, the baryons interact with each other and condense into lower energy states. In the process, the particles dissipate away a considerable amount of energy in the form of radiation which may escape the halo. The less energetic baryons then condense and undergo further infall towards the centre. The changing shape of the baryon density profile affects the DM profile by increasing density near the centre. The analytical estimations of this effect have been worked out in the adiabatic contraction regime [50]. In this regime, the DM particle orbits are assumed to be circular or nearly circular and the total mass enclosed by the orbit is assumed to be changing very slowly compared to the orbital time period of the DM particle. In this adiabatic regime, the invariance of pdq implies Here M (r) is the total mass enclosed inside the orbit of radius r. Using this invariance, an analytical estimate has been obtained which fairly matches with the numerical N-body simulation results [50,51]. The main effect is the steepening of the DM density profile near the centre and forming a denser core. As more baryons fall towards the centre, the gravitational well becomes deeper and more DM particles are attracted inward. This increases the slope of the central density profile [51]. In our case, the DM component itself will have a dissipation or cooling mechanism through an upscattering and a subsequent decay of the excited state. This process is independent of and in addition to the baryonic cooling.
Hence the effect of halo cooling will presumably be more prominent in this scenario and one would expect more complexity and richness in the small-scale structure of DM. As a result larger portion of the parameter space can be constrained.
The rate of this new dissipation mechanism will mainly be given by the upscattering rate as the the decay is very fast and can be assumed to be prompt. Here we will give a rough estimate of the rate of energy loss due to the upscattering and decay from the excited state. In the limit of nonrelativistic DM and ∆ M , the net kinetic energy lost per particle is approximately equal to ∆ itself. The upscattered χ 2 particles will decay and produce lighter particles with some amount of kinetic energy from the phase space available. One can estimate the leading order contribution to this energy gain to be O(∆ 2 /M 2 ) and O(v 2 2 ∆/M ) where v 2 is velocity of the upscattered χ 2 particles prior to decay. Therefore for all relevant parameter choices, the gain in the kinetic energy from the decay is negligible relative to the energy loss from the upscattering. The requirement for the upscattering and the decay to happen constrains the parameter space as µ 1 v 2 /4 > ∆ > m ρ . We shall assume that all light particles generated from the decays leave the halo.
In a halo, the average rate of energy loss in a DM shell of radius r and width dr, is estimated by The radial dependence of DM velocity could be estimated from simple Newtonian dynamics. It peaks around the where ρ s and r s are the scale density and radius, respectively. The individual DM velocities at a given radius will follow a thermalized Maxwell-Boltzmann (MB) distribution characterized by a virial velocity dispersion v(r). Essentially, in a fully virialized halo, the high energetic DM particles will most often occupy the outer edges of the halo. The halo cooling rate will be given by a convolution over the DM velocity distribution function where we take f (v) to be approximated by a Maxwell 3 . Note that the velocity distribution f (v) also depends on radial distance r throughv(r).
An approximate radial dependence of the cooling rate dE/dt for a halo of the size of that of the Virgo cluster is shown in Fig. 6. The profile was taken to be an NFW with a scale radius r s = 560 kpc and density ρ s = 3.2 × 10 5 M /kpc 3 . For simplicity the inelastic cross-section was taken to be velocity-independent constant σ in /M = 1 cm 2 /g. The resulting cooling rate shows a strong radial dependence and is largest near the virial radius. This cooling rate is to be compared with the energy inflow from the gravo-thermal collapse of the DM particles and due to the heat diffusion through selfscattering. The gravitational collapse brings faster (hotter) particles from the outer region of the halo to cooler inner part. And the scattering between the particles help diffuse the kinetic energy from the hotter particles to the relatively colder ones. The process of gravo-thermal collapse can be modelled following the Refs. [52,53]. The negative specific heat of a halo after virialization leads to this collapse. If we treat the DM particles as a fluid, the heat radiated inward at some radius r is given by Here the two terms within parenthesis on the RHS correspond to two different mean free path regimes. The first term describes the hard sphere scattering with the dimensionless coefficient a = 16/π. The second term describes the short mean free path regime which is proportional to the gravitational constant and the numbers b = 25 √ π/32 and C ≈ 0.75. Typical values of this heat inflow rate are 2-3 orders of magnitude larger than the cooling rate discussed above. Nevertheless, in the resonant region of the parameter space this halo cooling is expected to be efficient enough to distort the halo.
Upscattering and decay do not start abruptly but rather are continuous processes which will be present during the virialization process of the halo. At the initial epoch of structure formation the DM particles are highly nonrelativistic and there will be no dissipation. After DM falls towards the centres of the potential wells, acquires more energy and inelastic collisions become possible, it leads to cooling. From Eq. (18) it is clear that the inelastic scattering is a rather slow process and the halo will virialize at a faster rate than the dissipation. As a result, subsequent changes in the halo shapes is expected to be continuous and not episodic. A more detailed study of the effect of this new cooling mechanism will require an N-body simulation with this extra energy loss implemented in the dark sector [54].
A similar halo cooling mechanism was considered in Ref. [39] in the context of an atomic DM model. There neutral atomic dark hydrogen makes the DM abundance in the present Universe. The hyperfine splitting in the ground state of the dark atom leads to inelasticity in the system and the excited state decays to the ground state emitting massless dark photons. The masslessness of the dark photon implies that this cooling mechanism is more important for smaller halos because of their lower gravitational binding energy. On the contrary, in our case the particle ρ is massive. Hence the cooling mechanism shuts off for small mass halos where the DM particles do not have enough energy to upscatter, and the dissipation arises mainly in large galaxies or clusters. Note that the details of the particle physics model do not affect the radial dependence shown in Fig. 6, and all such details are encapsulated into the velocity dependence of the crosssection that determines this feature.

B. Drag and Evaporation from Inelastic Scattering
Independent bounds on DM scattering can be obtained from particle evaporation during collision of clusters and the movement of smaller dwarf-size halos through larger halos [55,56]. As was pointed out in Ref. [55], the SIDM particles will experience collisions in colliding clusters, whereas the stellar components of the objects will move freely without any appreciable friction. If the momentum transfer in a DM-DM collision is such that the final velocity is larger than the escape velocity of the parent halo then it would leave the halo and would result in DM evaporation from the halo. The existing observations from colliding clusters put strong constraint on this process yielding an upper bound on the DM self-scattering crosssection. An estimate of the of rate such collisions can be obtained following Ref. [55], in the limit of long-range interaction (as the hierarchy µ 1 v 2 /4 > ∆ > m ρ is easy to satisfy with smaller value of m ρ even at cluster size scale). In Ref. [55], the cumulative evaporation rate was shown to be more important than the immediate evaporation when DM has long range self-interaction. This rate is given by Here v 0 is the relative velocity between the two colliding clusters and ρ DM is the DM density in the bigger halo. The parameter θ min encodes the screening length and regulates the forward divergence. Because of this evaporation rate, the clusters will feel a drag force given by where η is an O(1) numerical factor depending on the nature of the mediator. In the last equality, following Ref. [55], we have defined the cross-sectionσ as The existing bound onσ from the abundance of dwarfs in our MW halo is very strong,σ/M 10 −11 cm 2 /g [55]. For two-level DM, two distinct cases may arise. Firstly the usual evaporation of DM particles is still feasible in this model, and has contributions from both elastic and inelastic scatterings. If the velocities of the scattered particles are larger than the escape velocities then they can escape the halo and would collectively cause dynamical friction between the halos. Second, inelastic scattering and subsequent decay provides an additional way for energy dissipation and gives an additional contribution to the drag force. For simplicity, if we assume that all DM particles are moving at the same velocity v 0 as the halo then n χ σ in v 0 is the upscattering rate per unit time. After each upscattering and decay event, two light particles escape the halo taking away an amount of energy which is roughly E decay 2∆. Therefore, the halo loses energy at a rate dE/dt, The resulting drag force per unit DM mass (or deceleration) due to this energy loss is given by Then the net drag force acting between the halos is given by The first term on the r.h.s above represents the cumulative evaporation rate, due to elastic and inelastic processes that are approximately equal across a large portion of the parameter space. We neglect the tiny velocity gain of χ 1 from the decay as we have seen it to be of even smaller order of magnitude in the previous subsection. The second term corresponds to the new dissipation mechanism from upscattering and decay. The quantity E decay denotes the energy loss rate averaged over the phase space of the final particles which, in the last equality, has been approximated to E decay 2∆. For simplicity here we have assumed that all DM particles in the incident halo have velocity v 0 . Of course a more careful analysis would require averaging over a Maxwellian distribution characterized by a velocity dispersion v 0 .
The relative size of the new term in Eq.(30) compared to the old term is given by ∼ 4v 2 0 ∆/M 10 −8 for M = 10 GeV, ∆ = 1 MeV and v 0 = 1000 km/s and assuming σ in σ in . The parametric smallness of the new drag force term may be traced back to the smallness of the mediator mass. A light particle-mediated interaction has a negative power dependence on velocity and is enhanced at small velocities, whereas the new term is virtually velocity independent. This velocity dependence may be useful to extract the impact of the second term, relative to the larger first term. We leave this investigation to a more detailed study.
There may be other signatures of this energy loss process. For example, just as the baryonic energy loss processes like Compton scattering and bremsstrahlung are responsible for the collapse of the ordinary matter into disk-like structures forming the galaxies, for two-level DM, upscattering and subsequent decay processes help DM lose energy and can lead to the formation of a rotating dark disk in DM halo [57][58][59][60][61][62]. As another signature, the authors in Ref. [63], observed a discrepancy between the predicted positions of the splashback radii (see [64][65][66]) of cluster-size halos in simulation and the observational data [67,68]. This mismatch could in principle be addressed by this energy dissipation mechanism through DM inelastic scattering.

IV. SUMMARY & OUTLOOK
In this work, we have studied the self-scattering of a two-level DM model. The off-diagonal interaction leads to inelastic scattering of a pair of DM particles from the ground state to the excited state, in addition to the ordinary elastic scattering.
If the incoming energy of the particles is below threshold, the excited state is not produced as physical states. Nevertheless, those states are produced offshell in the intermediate steps of the scattering and can affect even the elastic scattering cross-section. It was shown that the equations in this case, can be rotated to a new basis where the potential matrix becomes diagonal, and because of adiabaticity can be solved as a single state system with an appropriate potential.
When the incoming particles are above threshold, inelastic scattering may also take place. We have shown that in a large part of the parameter space, the inelastic cross-section is comparable to its elastic counterpart. This large inelasticity is a result of the maximal adiabatic mixing between the two states. We have also identified the Born and resonant regions in the relevant parameter space, and an estimate for the resonance condition has been given using a mapping of the two-level system to an equivalent one-level equation.
The off-diagonal interaction between the DM states allows the heavier state decay to the lighter one and the mediator. The upscattering and subsequent decay thus provides a mechanism for energy dissipation in DM halos. Assuming the decay to be prompt, the rate of the upscattering induced decay is given by the inelastic scattering rate which we computed to be 1-2 orders of magnitudes larger than the age of the Universe. Therefore, the DM halos can not condense into smaller halos via this mechanism. Rather, the inelastic process takes place only in larger objects and is effective only after the DM density becomes large enough at the centres of those objects. We compared this cooling rate with the heating due to ordinary elastic scattering and found that in some regions of the parameter space, the cooling rate could be a large fraction of the heating rate. We expect that this will leave an observable imprint on DM halo formation and evolution which can be only be probed by an N-body simulation incorporating this dissipative feature.
The same dissipation gives rise to an additional drag force between two colliding halos or for a small halos drifting through a larger one. When two halos collide with each other, the self-interacting DM particles scatter with each other and lose energy by emitting the light scalars. This energy loss can be interpreted as a new drag force acting between the halos. We calculated an analytical expression for this new drag force and found that it is small relative to the other contribution from ordinary scattering, but has a distinctive velocity independence unlike the usual drag force.
This work was partially funded through a Ramanujan Fellowship of the Dept. of Science and Technology, Government of India, and the Max-Planck-Partnergroup "Astroparticle Physics" of the Max-Planck-Gesellschaft awarded to B.D and has received partial support from the European Union's Horizon 2020 research and innovation programme under the Marie-Sklodowska-Curie grant agreement Nos. 674896 and 690575. The work of A.D. was supported by the Dept. of Atomic Energy, Government of India.
Here K is the reaction matrix and [J (kr)] ab = +k a r j (k a r)δ ab , above threshold, = +k a r ι (k a r)δ ab , below threshold, [N (kr)] ab = −k a r n (k a r)δ ab , above threshold, = −k a r κ (k a r)δ ab , below threshold. (A4) Here j (x) and n (x) denote spherical Bessel functions of first and second kinds, and ι (x) and κ (x) are the modified spherical Bessel functions of first and second kinds, respectively. These two types of functions serve as the asymptotic forms of the wavefunction as indicated above. In the below threshold case the boundary conditions need to be changed for the excited state. In Ref. [70,71], the author has shown that in the below threshold case, only the open-open part (the part which consists of only the open channels) of the K matrix contributes to the final scattering matrix though one has to solve the full system of Schrödinger equation. In this case the asymptotic wavefunctions are either exponentially growing or decaying which may cause trouble in the numerical computation (see the second line in Eq.(A4)). It is solved by normalizing the closed channel wavefunctions and their derivatives by J and N respectively such that the new asymptotic wavefunctions become J (kr) → 1, J (kr) → J (kr)/J (kr) and similarly for N (kr).
We solve for K from Eq.(A3) by taking logarithmic derivative of the equation, This S is computed for all partial waves starting from = 0 to max . As stated in the text, the value of max depends on the initial momentum of the particles and the range of the potential. A useful lower bound on its value can be given as max ≥ k/m ρ for the case discussed in this paper. The final total cross-section is given by (2 + 1)|(T ) ab | 2 .
(A7) where (T ) ab = (T ) ab + (−1) (T ) a b . Here the prime on a denotes an exchange of particles in the final 2-body DM state. Note that the last term in Eq.(A7) is present only when the final state particles are identical. In case of distinguishable particles, this term will be absent and so will be the extra factor of 1/2. The other two quantities of interest are the transfer and viscosity cross-sections. The definition of the transfer cross-section σ T is given in Eq. (7). Expanding the differential cross section in the partial wave basis gives (A9)