Astrophysical probes of inelastic dark matter with a light mediator

We consider an inelastic dark matter model, where a fermion is charged under a broken U(1) gauge symmetry, and introduce a tiny Majorana mass term to split the fermion into two states with the light one being a dark matter candidate. If the gauge boson is light, it can mediate both elastic and inelastic dark matter self-interactions in dark halos, leading to observational consequences. Using a numerical technique based on partial wave analysis, we accurately calculate the elastic and inelastic self-scattering cross sections. We assume a thermal freeze-out scenario and fix the gauge coupling constant using the relic density constraint. Then, we focus on six benchmark masses of dark matter, covering a wide range from $10~{\rm MeV}$ to $160~{\rm GeV}$ and map parameter regions where the elastic scattering cross section per unit mass is within $1~{\rm cm^2/g} - 5~{\rm cm^2/g}$, favored to solve small-scale issues of cold dark matter. If the heavy state can decay to the light state and a massless species, the inelastic up-scattering process can cool the halo and lead to core collapse. Taking galaxies with evidence of dark matter density cores, we further derive constraints on the parameter space. For dark matter masses below $10~{\rm GeV}$, the mass splitting must be large enough to forbid up scattering in the dwarf halo for evading the core-collapse constraint; while for higher masses, the up-scattering process can still be allowed. Our results show astrophysical observations can provide powerful tests for dark matter models with large elastic and inelastic self-interactions.


I. INTRODUCTION
It is well established there is a non-luminous component in the universe, called dark matter. Since its influence has, so far, only been observed to be gravitational in nature, the particle properties of dark matter are still largely elusive. Self-interacting dark matter (SIDM) is a class of particle physics models in which dark matter particles are assumed to have a large self-scattering cross section [1,2], see [3] for a recent review. N-body simulations of structure formation demonstrate strong dark matter selfinteractions can lead to heat transfer in the halo and allow its inner region to thermalize [4][5][6][7][8][9][10]. Recently, it has been shown that SIDM can explain a number of long-standing puzzles in astrophysics presented in the prevailing cold dark matter theory, such as diverse galaxy rotation curves of spiral galaxies in the field [11][12][13][14], dark matter distributions in satellite galaxies in the Milky Way [5,[15][16][17][18][19], and shallow inner dark matter density profiles in galaxy clusters [2,20]. SIDM also inherits all the success of cold dark matter in explaining large-scale structure of the universe [21] and many important aspects of galaxies [22,23].
Many SIDM models assume there is only one dark matter state and a light force carrier mediates elastic dark matter selfscattering in the halo, see, e.g., [3]. In this case, the initial and final states in the collisions are the same, and they only redistribute energy of dark matter particles as the halo as a whole does not lose energy. More recently, there is growing interest in considering particle physics realizations of SIDM with multiple states. For example, to avoid strong bounds from direct detection experiments [24][25][26][27], Refs. [28,29] propose an inelastic SIDM model, where there are two dark states and they differ by a small mass splitting [30,31]. One can adjust the splitting to kinematically forbid transitional up-scattering in nuclear recoils, but still allow strong elastic self-interactions between two light states in the halo. In addition, Ref. [32] studies dark matter self-interactions in the exciting dark matter model [33][34][35][36], where dark matter collisions can produce a heavy state that subsequently decays back to the light one and a standard model particle. More generally, if the SIDM candidate is made of composite states, such as dark atoms [37][38][39][40][41][42][43][44][45][46][47] and strongly-coupled particles [48][49][50][51], it is natural to expect inelastic excitations during dark matter collisions in dark halos.
In this paper, we consider an inelastic SIDM model and study its astrophysical implications. It assumes that a Majorana mass term induces a small mass splitting between two fermionic dark matter states and they interact with a U(1) gauge boson. We assume the gauge boson is light and develop a numerical method to calculate both elastic and inelastic dark matter selfscattering cross sections. After imposing the relic abundance constraint on the gauge coupling constant, we focus on benchmark dark matter masses, which cover a wide range from 10 MeV to 160 GeV, and search for parameter regions where the elastic self-scattering cross section per unit mass (σ V /m χ ) satisfies 1 cm 2 /g ≤ σ V /m χ ≤ 5 cm 2 /g, as favored by observations on galactic scales [3]. Our work is a natural and simple extension to the minimal SIDM model that contains only one dark matter state, but calculations of the self-scattering cross sections in the current model are much more challenging than the minimal one [52,53]. In addition, we explore a broader mass range for both dark matter and mediator particles, compared to earlier studies [29]. As we will show, for the dark matter mass below ∼ 1 GeV, inelastic up scattering dominates over elastic one if the former is kinematically open in the dark halo. This has important implications for constraining the parameter space.
We further consider the endothermic up-scattering process of dark matter particles and its influence on halo evolution and inner halo structure. If the heavy state decays back to the light state by releasing a massless species, the SIDM halo profile can become cuspy again, because the dissipative self-interactions can cool the inner halo and speed up the onset of core collapse [54,55]. Using dwarf galaxies that show density cores, Ref. [54] derives constraints on parameters that characterize the cooling rate of dissipative dark matter collisions. In this work, we will take the results in [54] to further narrow down parameter space of the inelastic SIDM model.
The paper is organized as the following. In Sec. 2, we present details of the model and outline the numerical method used in calculating elastic and inelastic dark matter self-scattering cross section. In Sec. 3, we present the astrophysical constraints on the parameter space. We conclude in Sec. 4. In the appendix, we provide details of substitutions and transformations exploited in this work for solving the Schrödinger equation with two states.

A. Scattering Cross Sections
We assume that the dark matter particle is a fermion (Ψ) and it interacts with a dark U(1) gauge boson (φ µ ). The model can be described by the following Lagrangian [32][33][34] where Ψ c is the charge conjugation of Ψ, φ µν is the field strength of φ ν , g χ is the gauge coupling constant, m is the Dirac mass of the dark matter state and ∆m is its Majorana mass. In this work, we assume m ∆m. Defining the Majorana mass eigenstates as χ 1 = i(Ψ − Ψ c )/ √ 2 and χ 2 = (Ψ + Ψ c )/ √ 2, we rewrite equation (1) as where m χ = m − ∆m/2 is the mass of the light state χ 1 .
Since the mass eigenstates are Majorana states, they carry no charge and only interact through an off-diagonal coupling. The relevant Feynman diagrams for both elastic and inelastic dark matter self-scattering are shown in Fig. 1. The tree-level elastic scattering process involves mixed initial and final states. All other elastic scatterings occur through high-order box or ladder diagrams. As the universe cools, the up-scattering process becomes kinematically unfavorable, driving the density of the heavy state down [29,56]. Thus, we assume dark matter is made of the light state in the halo and only consider elastic and inelastic scattering processes shown in the left and right panels of Fig. 1, respectively.
In the non-relativistic limit, we can apply the Schrödinger formalism. There are two wave functions coupled by a matrix potential of the form where µ = m χ /2 is the reduced mass, the vectorψ T = ψ 1 , ψ 2 detonate the wave functions for the two particle modes, and the matrix potential V is We have defined α χ ≡ g 2 χ /4π as the dark fine structure constant. The energy needed to create the heavy state as a pair is 2∆m. The numerical solution to this set of coupled differential equations gives the scattering cross sections through the method of partial waves.

B. Numerics
We assume that dark matter freezes out in the early universe with the relic abundance to be consistent with the observed density. In this paper, we set the dark fine structure constant to α χ = 0.01 (m χ /270 GeV) such that the annihilation cross section is 6 × 10 −26 cm 2 /g. The dark matter self-scattering cross sections, both elastic and inelastic, are in general velocitydependent. To capture the relevant physics on dwarf scales, we set the dark matter relative velocity to be 60 km/s in the halo throughout this paper unless otherwise stated. The model is left with three free parameters, the dark matter mass m χ , the mass splitting ∆m and the mediator mass m φ .
Performing separation of variables on equation (3), we have the radial equation where l is the angular momentum mode, k is the magnitude of the wave vector, R l,i (r) are the radial wave functions for i = 1, 2 and V i,j denotes components of the matrix (4). Defining the following dimensionless parameters and substitutions, we rewrite the radial equation (5) in the matrix form as The wave function can be expanded in terms of spherical waves, where F x/y,l are the scattering amplitudes for the two particle system and p in = a, c, depending on the initial state. If the initial state is χ 1 as we consider in this work, p in = a andψ T in = 1 , 0 . The differential cross section is given by, where p in , p out and F l again depend on the initial and final states. For elastic scattering To find the scattering amplitudes F x/y,l , we need to first find the wave function by numerically solving equation (7), and then map its form at large radii onto the spherical wave expansion in equation (8). However, a direct numerical solution to the wave equation (7) is unstable for a large part of the parameter space of interest. To tame these instabilities, we follow the procedure discussed in [29,57] and make a number of substitutions to transform the wave equation into a more manageable form; see the appendix for details. In this work, we calculate the viscosity cross section for dark matter self-interactions [53], which regulates both forward and backward scatterings. See the appendix for an explicit expression of the viscosity cross section in terms of phase shifts. In Fig. 2, we show the elastic (solid) and inelastic (dashed) self-scattering cross sections vs the relative velocity for a few representative cases, where we choose a wide range of dark matter masses from 10 MeV to 160 GeV. Overall the cross sections decrease as the velocity increases. For dark matter masses below 1 GeV (left panel), there is a clear indication of the threshold velocity below which inelastic up scattering is kinematically forbidden. In this mass range, when the coupling constant is set by the relic abundance constraint [58][59][60][61], i.e., α χ = 0.01(m χ /270 GeV), the inelastic scattering cross section is much larger  than the elastic one as long as the up-scattering channel is open. As we discussed, in this model, inelastic up scattering occurs at the tree-level, while elastic scattering at the high-order level. For small m χ below 1 GeV, the dark fine structure constant α χ is small as well, and the non-perturbative quantum effect is absent to enhance the elastic cross section. To demonstrate this point, we present another case, where we set α χ = 0.01 for m χ = 40 MeV. In this case, both elastic and inelastic cross sections are similar for the velocity larger than 18 km/s. For high dark matter masses (right panel), the elastic and inelastic cross sections become more compatible, aside from resonance peaks, and the up-scattering process is kinematically allowed in the plotted velocity range. As m χ increases, the gauge coupling α χ increase accordingly and the non-perturbative effect boosts the elastic scattering cross section significantly.

III. ASTROPHYSICAL CONSTRAINTS
Taking the benchmark m χ values shown in Fig. 2, we scan parameter space of the ∆m-m φ plane, such that the elastic cross section fall within the range of 1 cm 2 /g ≤ σ V /m χ ≤ 5 cm 2 /g for the relative velocity 60 km/s, a characteristic value for dwarf galaxies that prefer a dark matter density core. In Fig. 3, we show the resulting parameter space (shaded). For a given dark matter mass, there is a preferred range in the plane, where the elastic self-scattering cross section is large enough to thermalize the inner halo in accord with observations [3]. For all cases, if the mass splitting ∆m is small, the mediator mass m φ is almost a constant. While, as ∆m increases towards the high end, m φ must decrease to preserve the elastic cross section in the desired range, since the elastic process involves virtual up-scattering processes, as shown in Fig. 1 (left). The transition occurs when the mass splitting reaches the kinematic threshold, where up scattering is forbidden for larger values of ∆m, i.e., 2∆m = µv 2 rel /2 with v rel = 60 km/s. In Fig. 3, the orange shaded regions are where χ 1 χ 1 → χ 2 χ 2 is kinematically forbidden, while in the magenta and blue regions, the up scattering is allowed. Note in the case of m χ = 40 GeV there is more than one branch for the favored parameter space, because the scattering is in the strong resonance regime [53,62] and multiple ranges of the mediator mass is allowed; see also [29].
If the mass splitting is large enough and up scattering is forbidden, dark matter self-interactions are purely elastic and the condition of 1 cm 2 /g ≤ σ V /m χ ≤ 5 cm 2 /g is sufficient enough to specify astrophysical constraints. However, if χ 1 χ 1 → χ 2 χ 2 is allowed in the halo and the resulting χ 2 can further decay to χ 1 and some light species, this dissipative process may cool the inner halo and speed up the SIDM core collapse [54]. For the model we consider, m φ ∆m in the parameter regions of interest, e.g., the shaded regions in Fig. 3, hence the decay process χ 2 → χ 1 φ is kinematically forbidden. On the other hand, if we consider a more general setup, there are other interaction terms that may lead to dissipative decays of χ 2 . For example, Ref. [36] introduces a dimension-5 dipole operator (1/M )χ 2 σ µν χ 1 F µν , where M is the cut-off scale and F µν is the field strength of the standard model photon. With this operator, χ 2 can decay to χ 1 and γ. The rate is Γ χ2→χ1γ = 4∆m 3 /(πM 2 ), m χ = 10 MeV The shaded regions correspond to the accessible parameter space in the ∆m-m φ plane, where the elastic scattering cross section is in the range of 1 ≤ σV /mχ ≤ 5 cm 2 /g, favored by solving the small-scale issues. We have set the dark matter relative velocity to be 60 km/s, a characteristic value in dwarf galaxies that prefer a dark matter density core. In the orange regions, the up-scattering process (χ1χ1 → χ2χ2) is kinematically forbidden and dark matter self-interactions are purely elastic. In the magenta (blue) regions, the up-scattering process is allowed, and the dissipation process associated with the χ2 decay can lead to core collapse in dwarf galaxies with a timescale shorter (longer) than 10 Gyr. The starred and triangle points are references which show the mapping between the elastic scattering and the core-collapse constraints, as explicitly shown in Fig. 4. and χ 2 's lifetime is τ = 1/Γ χ2→χ1γ ∼ 0.5 sec (M /TeV) 2 (keV/∆m) 3 . For ∆m ∼ 10 −3 eV, it is comparable to the age of galaxies, ∼ 10 Gyr, for M up to ∼ 1 TeV. Thus, this dissipative decay is relevant to halo dynamics if the dipole operator is present. In addition, in atomic dark matter models, an excited atomic state can decay to a ground state by emitting a massless Mapping between elastic scattering and galaxy core-collapse constraints on the σ /σ-ν loss plane. In the gray regions within the black contour, the halo core collapse induced by dissipative dark matter self-interactions occurs within the age of galaxies; adapted from [54]. The dashed black lines are the linear extrapolation of the contour for σ /σ > 1. The magenta points that lie within the gray regions are disfavored by the core-collapse constraints; while the blue points outside are still allowed. The stars and triangles are reference points and their correspondences are shown on the ∆m-m φ plane in Fig. 3 dark photon.
In what follows, we assume that χ 2 can decay to χ 1 and a massless species that escapes the halo, and study additional astrophysical constraints on the parameter space. Using dwarf galaxies that show shallow density cores, Ref. [54] derives bounds on dissipative dark matter interactions by demanding the core-collapse timescale longer than the age of galaxies, ∼ 10 Gyr. In particular, it uses the energy loss per collision and the ratio of inelastic to elastic cross sections, i.e., ν loss = E loss /m χ and σ /σ, respectively, to characterize the cooling effect, and places constraints on their combinations. To apply the core-collapse constraints on our model, we set E loss = ∆m, calculate ν loss and σ /σ values for each favored model point shown in Fig. 3 (shaded), and then compare them with the limits on the σ /σ-ν loss plane from [54] as reproduced in Fig. 4 (gray shaded). In the magenta shaded regions of Fig. 3, the dissipative self-interactions are strong enough to cause core collapse in dwarf halos within 10 Gyr. While in the blue regions, inelastic up scattering can occur, but the overall cooling rate is small to trigger core collapse in the age of galaxies.
To better understand these constraints, we show the distribution of the model points in the σ /σ-ν loss plane for three benchmark cases in Fig. 4, along with the bounds from [54] (gray). All points (magenta) that lie within the gray regions are disfavored as they result in a core-collapse timescale too short to fit the observations; while the points (blue) outside are still allowed. We classify the model points shown in Fig. 3 using the same color scheme. Note we have extrapolated the disfavored parameter space following the trend beyond the upper limit of σ /σ in [54] (black dashed). This is reasonable, because the bounds should be stronger as σ further increases.
From Fig. 4, we see that as the dark matter decreases from 160 GeV, more of the parameter space is disfavored by the corecollapse constraints. When the mass approaches 10 GeV or smaller, the entire model points lie within the gray regions. Since the gauge coupling reduces as the dark matter mass decreases, the inelastic scattering gradually dominates over the elastic one if the former is open. For the dark matter below ∼ 10 GeV, only the portion of the parameter space, where inelastic scattering is kinematically forbidden, remains viable. While for the cases of m χ = 160 GeV and 40 GeV, in some parts of the parameter space, inelastic up scattering is allowed, but the cooling rate is not significant so they evade the collapse constraints, because either the inelastic cross section or the energy loss per collision is small. We demonstrate this by using the reference points in both Figs. 3 and 4 (stars and triangles), which are in one-to-one correspondence. In the case of m χ = 160 GeV, the two references have similar ∆m and σ/m, but the star point has much smaller σ /m than the triangle one as the former is closer to the threshold of the up scattering. While in the case of m χ = 40 GeV, the reference points mainly differ in ∆m, resulting different locations in the σ /σ-ν loss plane.

IV. CONCLUSION
We have studied an inelastic dark matter model with a light mediator, based on a U(1) gauge symmetry. The presence of a small Majorana mass splits a Dirac fermion into two Majorana states and the light one is the dark matter candidate. In this model, both elastic and inelastic dark matter self-interactions can be present in the halo. The former is mediated by a high-order process with multiple exchanges of the light mediator; while the latter by a tree-level process when it is kinematically allowed.
Using the technique of partial waves, we have developed a numerical procedure to calculate the elastic and inelastic scattering cross sections.
We have explored astrophysical constraints on the model parameters, i.e., the dark gauge coupling constant, dark matter mass, mediator mass and mass splitting between the two Majorana states. We first imposed the relic density constraint on the coupling constant by assuming the standard freeze-out scenario, then chose six benchmark cases that cover a wide range of the dark matter mass, 10 MeV-160 GeV. For each case, we have found the parameter regions, where the elastic scattering cross section falls within the range of 1 cm 2 /g-5 cm 2 /g in dwarf galaxies in order to solve the small-scale issues. Our analysis shows that if the mass splitting gets too large, the kinematic suppression of the intermediate virtual processes demand that the mediator become lighter to preserve the desired elastic cross section. We also found that when the dark matter mass decreases the inelastic scattering cross section dominates over the elastic one. This is because the coupling constant becomes smaller as the mass decreases and the non-perturbative quantum enhancement for the elastic cross section diminishes accordingly.
If the heavy state can decay to the light state and a massless degree of freedom, inelastic dark matter self-interactions may induce a dissipative process that cools the inner halo and leads to SIDM core collapse. Observations of dark matter density cores in many low surface brightness galaxies put a constraint on the rate of energy loss. We studied its implications for the dark matter model we consider, and found it eliminates the majority of the parameter space for dark matter masses below ∼ 10 GeV, unless the mass splitting is large enough so that the up-scattering process is forbidden. For a higher mass, there are parameter regions where the model evades the core-collapse constraints while the inelastic scattering is kinematically allowed. Our work demonstrates that astrophysical observations can provide powerful tests for inelastic dark matter models with a light mediator. The analysis can be used to constrain models where dark matter is made of a composite state, such as dark atoms and nuclei. It is also interesting to test those models using observations of galaxy clusters that show evidence of a density core in their inner halos [2,20].
Appendix A: Cross Section Calculations

Reformulation
In this appendix, we outline the variable phase space approach to reformulate equation (7); see [57] for a more detailed discussion on this method. The basic idea is to build a solution to equation (7) using solutions to the free-particle case (α χ → 0). In the non-interacting limit, the solutions can be written as superpositions of spherical Bessel and Neumann functions with constant coefficients. To build solutions to (7), we use superpositions of free-particle solutions and upgrade the coefficients to functions, i.e., where χ i (x) are numerical functions, p i = a, c, depending on the particle state, and f (l) (p i x) and g (l) (p i x) are the free-particle solutions. They obey the differential equation, where z takes the place of f or g. The function f (l) (p i x) is defined to be regular at the origin and g (l) (p i x) is irregular as x → 0.
To form a general solution to equation (7), we must solve two coupled second-order differential equations. We therefore require four linearly independent solutions. The expression (A1) represents only one of the four solutions but has four degrees of freedom; two of them come from α i (x) and the other two from the normalization of f (l) (p i x) and g (l) (p i x). We must impose constraints to reduce the extra degrees of freedom. Suppose that dχ i (x), which is trivially true for constant coefficients. This requires that We set the normalization of f (l) (p i x) and g (l) (p i x) by defining the Wronskian of the system to be After imposing the constraints, we have only one degree of freedom and an overall constant per linearly independent solution.
A consistent choice for f (l) (p i x) and g (l) (p i x) is where j l (p i x) is the spherical Bessel function and h (l) l (p i x) is the spherical Hankel function of the first kind. To keep track of the linearly independent solutions, we introduce a new subscript, where we have dropped the angular momentum label l for brevity, n = 1, 2 for the two independent solutions for a given i = 1, 2, which labels the particle state. Defining we can rewrite equation (A1) in a compact form, Further defining we have Taking the x → ∞ limit of the choices for f (p i x) and g(p i x), we can see the virtue of the conventions and definitions employed so far, Inserting (A11) into (A10) and comparing with equation (8), one can find that the components of M are related to the scattering amplitudes as where the incoming scatterers are of type 1. Similarly, for incoming particles of type 2. The i th column of ξ is interpreted as the scattered wave functions for the two particle states where the incoming states are of type i. Next we make the definition, Using the formalism developed in this subsection, we can derive the following first-order different equation for U ij (x), and g (p i x) ≡ dg(p i x)/d(p i x). As x → 0, β ij (x) → 0 since the solution χ ij (x) must be regular at the origin and we take α ij (x) → δ ij . Therefore, M ij (x → 0) = 0 and the initial condition for U ij (x) becomes, The advantage of this differential equation is that only logarithmic derivatives of the free solutions enter into the equation greatly increasing its numerical stability. We can now solve equation (A15) using the initial condition (A16) for U ij (x). Once U ij (x) is known then the scattering amplitudes ∼ M ij (x) can be obtained using the definition (A14). Finally, the scattering cross section can be calculated using equation (9). It is useful to note that the transformation ∆m → −∆m changes the incoming particles from one type to the other (a ↔ c). Therefore we only need to solve for M 11 (x) and M 21 (x), once ∆m is changed to −∆m, in order to obtain all scattering cross sections. Also, there is an equation for α ij (x) (β ij (x)) which carries information needed to solve for the Sommerfeld enhancements but are not needed in calculating the self-scattering cross sections [29].