Self-Scattering for Dark Matter with an Excited State

Self-interacting dark matter scenarios have recently attracted much attention, as a possible means to alleviate the tension between N-body simulations and observations of the dark matter distribution on galactic and sub-galactic scales. The presence of internal structure for the dark matter --- for example, a nearly-degenerate state in the spectrum that could decay, or be collisionally excited or de-excited --- has also been proposed as a possible means to address these discrepancies. Such internal structure can be a source of interesting signatures in direct and indirect dark matter searches, for example providing a novel explanation for the 3.5 keV line recently observed in galaxies and galaxy clusters. We analyze a simple model of dark matter self-scattering including a nearly-degenerate excited state, and develop an accurate analytic approximation for the elastic and inelastic s-wave cross sections, which is valid outside the perturbative regime provided the particle velocity is sufficiently low (this condition is also required for the s-wave to dominate over higher partial waves). We anticipate our results will be useful in incorporating inelastic self-scattering into N-body simulations, in order to study the quantitative impact of nearly-degenerate states in the dark matter spectrum on galactic structure and dynamics, and in computing the indirect signatures of multi-state dark matter.


Introduction
The verification of the existence of dark matter (DM) on wildly disparate scales is one of the greatest triumphs of modern astrophysics. It also poses one of the greatest outstanding questions of modern particle physics, due to the lack of a viable DM candidate within the Standard Model (SM). As one of the most promising potential windows onto new fundamental physics, there are many efforts underway to probe the particle nature of DM and its interactions with the SM, but to date it has only been detected through its gravitational properties. The formation and dynamics of DM structures can be highly sensitive to the microphysics of DM. There are stringent constraints on DM interactions with SM particles, but the interactions of DM within its own "dark sector" (at most weakly coupled to the SM) are far less constrained. Increasingly detailed probes of the distribution of dark matter in halos may allow us to explore and constrain such interactions, independent of the coupling of the dark sector to known particles. It is therefore important to understand the potential observational signatures of dark sector physics.
It is possible these discrepancies reflect the gravitational influence of baryonic matter, which is not included in CCDM-only simulations [10][11][12][13][14][15]. However, whether baryonic effects can explain all the observations remains an open question, and the discrepancies might also be clues to DM microphysics. DM self-interaction (via some novel dark sector physics) can alleviate the discrepancies by increasing the scattering of DM particles out of the dense central regions of halos [16]; recent simulations have confirmed the efficacy of this process [17][18][19][20]. Consistency with existing limits on dark matter self-interaction from large scale structures (which have much greater virial velocities than dwarf galaxies) is most easily achieved if the scattering cross-section has a velocity dependence (growing larger at low velocities); however, constant cross sections in the range σ/m χ ∼ 0.1 − 1 cm 2 /g remain viable [18].
Nearly-degenerate excited states are another simple modification to DM physics with potentially striking observational signatures, especially in the presence of non-negligible selfinteraction. Exothermic scatterings from an excited state would deliver velocity "kicks" to DM particles, which could be comparable to or larger than the escape velocity in bound DM structures, especially in the slow-moving environs of dwarf galaxies. In this way inelastic DM scattering can dilute dense cusps or dissipate dwarf halos entirely [21]. (Velocity kicks from late-time decays of a metastable excited state have been considered in e.g. [22][23][24][25][26][27].) Inelastic DM scattering from either DM or baryonic matter, exothermic or endothermic, can also have interesting signatures in direct and indirect searches for DM (e.g. [28][29][30] for direct detection, [31][32][33] for indirect detection).
Such scenarios are easily constructed from a theoretical perspective. If DM is charged under some new dark gauge symmetry which is broken, then the states in the dark matter multiplet can naturally acquire a small mass splitting, regardless of whether the gauge group is Abelian or non-Abelian (e.g. [34][35][36][37]). In any case where the mediator of the DM selfinteraction is light, the scattering cross section is automatically enhanced at low velocities, ensuring a larger impact on dwarf galaxies relative to Milky-Way-size galaxies and clusters.
Previous studies of self-interacting DM (both analytic and simulation-based) have fo-cused on the case where only degenerate dark matter states participate in the interaction (i.e. either DM is a single Majorana fermion, or a Dirac fermion with particle-particle and particle-antiparticle scatterings). Even in situations where the scattering is purely elastic (i.e. there is no transition to a state of different mass), as we will show in this article, the presence of a nearly-degenerate state in the spectrum can significantly modify the large resonances present in the low-velocity scattering rate. However, the addition of a second state adds at least one additional parameter to the problem (the mass splitting), making numerical analysis computationally expensive and, in parts of parameter space, unstable.
The main goal of this article is to work out an analytic approximation for the DM-DM scattering cross section induced by an off-diagonal Yukawa interaction, in the presence of a single excited state. This corresponds to the case where the dark gauge group is U (1) and provides a simple and illustrative toy model for inelastic dark matter self-scattering more generally. As we will show, our expressions give good agreement with numerically solving the appropriate multi-state Schrödinger equation, but are very quick to compute and provide intuitive insight into the scattering behavior. Our results are applicable to the low-velocity regime where s-wave scattering dominates the total cross section; we defer study of the "classical" regime, where many partial waves contribute, to later work. We identify regions of parameter space with particular relevance to dwarf-galaxy-sized halos, and consequently to the discrepancies described above. We also determine regions of parameter space which could potentially provide a DM explanation for the observed 3.5 keV spectral line from clusters [38,39] via collisional excitation followed by decay, as in the "XrayDM" scenario [33].
We begin, in Section 2, by describing our toy model and approach to computing the scattering cross sections (the method follows that of [40]). We discuss our analytic approximate scattering cross sections, and their behavior in several limiting regimes of interest, in Section 3. We simultaneously verify the validity of our approximate solutions by direct comparison to the numerical results. In Section 4 we provide first simple estimates of the parameter space where such scatterings could affect the internal structure and dynamics of dwarf galaxies, or yield the appropriate upscattering cross section to account for the 3.5 keV line in the XrayDM scenario. Concluding remarks follow in Section 5.

A Simple Model
We consider the case of a Yukawa-like interaction coupling two states with some small mass splitting, δ. We take the dark matter to be a pseudo-Dirac fermion charged under a dark U (1). At high energies, where the U (1) symmetry is unbroken, the dark matter is a charged Dirac fermion; at low energies, a small Majorana mass splits the Dirac fermion into two nearlydegenerate Majorana states χ 1 and χ 2 . Since these states cannot carry conserved charge, their couplings with the U (1) vector boson A D are purely off-diagonal: that is, there is no vertex of the schematic form We will not specify the high-energy physics that gives rise to the Majorana masses, since we restrict ourselves to the extremely non-relativistic environment relevant for present-day DM scattering. The phenomenology of such models in the context of dark matter annihilation and inelastic scattering has been discussed in e.g. [28,31,34,40,41]. Inelastic scattering between states with small mass splittings can also naturally be generated in other contexts, for example composite dark matter [42] and sneutrino dark matter [43].
Since we work always in the non-relativistic limit, the interaction between the DM-like states can be characterized entirely in terms of a matrix potential that couples two-body states. The scattering cross sections can then be calculated using the methods of quantum mechanics (we refer the reader to e.g. [44] for a field-theoretic discussion using the Bethe-Salpeter formalism, and [45] for a derivation using effective field theory). The potential matrix coupling the |11 and |22 two-body states (corresponding to both particles being in the ground state or both particles being in the excited state) is [40]: where α is the coupling between the dark matter and the mediator, m φ is the mass of the mediator, δ is the mass splitting between the ground and excited states, r is the separation between the particles, and the first row corresponds to the ground state |11 . The two-body Schrödinger equation for the relative motion (factoring out the overall free motion of the system as usual) then takes the form: Here m χ is the mass of the dark matter, v is the individual velocity of either of the dark matter particles in the center-of-mass frame (half the relative velocity), and Ψ( r) is the wavefunction. As mentioned above, the interaction between the ground state, |1 and the excited state, |2 , is purely off-diagonal. This is an automatic consequence of taking the force carrier to be a vector, as the mass eigenstates are 45 • rotations of the high-energy gauge eigenstates, and do not carry a conserved charge. As a result, two particles initially in the same state (ground or excited) can only scatter into the two-body states where they are both in the ground state, or both in the excited state. Such scatterings are described by the off-diagonal terms in the V(r) matrix. If the initial state is |12 (i.e. one particle is in the ground state and the other in the excited state) then their scattering decouples from the other two-body states and is elastic, with the final state being |12 or |21 . This case can be treated by the existing methods in the literature (e.g. [46]).
In this article we will treat only s-wave scattering (the methods developed in [40] do not generalize to higher partial waves). Accordingly we will make the ansatz of spherical symmetry in solving the Schrödinger equation (more formally, as discussed in [40], we expand in partial waves and only keep the = 0 term), and write rΨ( r) = ψ(r).
We define the dimensionless parameters: so that rescaling r by αm χ c/ gives the s-wave Schrödinger equation: In the last line we have defined the matrix potentialV (r) and the scalar Yukawa potential V (r) = e − φ r /r.

Calculating Approximate Wavefunctions
The eigenvalues λ ± and eigenvectors ψ ± of the matrixV (r) are There is a transition in the behavior of the eigenvalues and eigenvectors when e − φ r r ∼ 2 δ 2 . In the regime where e − φ r r 2 δ 2 , the eigenvalues and eigenvectors can be approximated as When e − φ r r 2 δ 2 , the eigenvalues and eigenvectors can be approximated as In other words, at small radii the symmetry is restored and the energy eigenstates become the gauge eigenstates, as at high energies; at large radii the energy eigenstates are the mass eigenstates. At small and large r, the diagonalization of the potential matrix is roughly independent of r. If φ ± (r) = λ ± φ ± (r), then φ ± (r) ψ ± is an approximate solution to the matrix Schrödinger equation in those regimes, and we can use standard techniques -in particular, the WKB approximation -to solve for the scalar wavefunctions φ ± (r). However, in the transition region where e − φ r r ∼ 2 δ 2 , the eigenvectors will vary as a function of r. Provided that 2 δ /2 φ , this transition occurs at r 1/ φ , where the exponential behavior of the potential dominates: in this case we can use an exact solution for the twostate exponential potential (taken from [47]) to cross the transition region and match the WKB solution at small r to the large-r asymptotic solution. Where 2 δ /2 > φ , we can still use this approach, but the exponential potential is a poor approximation to the true V (r) at the radius corresponding to the transition region, and so the results are less reliable.
The full procedure for determining the approximate wavefunction ψ(r) can be summarized as: • For r 1, we approximate the Yukawa potential as V (r) ∼ 1/r, assume the potential term dominates (since 2 v and 2 δ are assumed to be small), and solve the resulting Schrödinger equation exactly.
• For r 1, we propagate the small-r solution outward using a WKB approximation. The validity of the WKB approximation in the potential-dominated regime requires that V (r)/V (r) = 1 2 √ r e φ r/2 ( φ + 1/r) 1, which is just equivalent to requiring that the spatial variation of the local de Broglie wavelength is sufficiently gradual. For r > 1, the condition for validity of the approximation is that φ / V (r) 1, so the approximation breaks down when V (r) ∼ 2 φ .

Regimes for the wavefunction
Conditions on v , φ , and δ 2nd and 3rd regimes above overlap Scattering cross section enhanced by dark force v , φ , δ 1 Table 1. A summary of the various approximations and assumptions used for deriving the wavefunction in different regimes. The top four rows characterize the different regimes in r where various techniques can be applied to approximate the wavefunction in the vacinity of a Yukawa potential; the lower three rows describe constraints on the parameters which allow the approximation to hold and cause the scattering cross section to be enhanced. Despite these restrictions, the various regimes of validity overlap substantially, making these approximations useful in large swatches of parameter space.
• Where the WKB approximation breaks down (at V (r) 2 φ ) or where the diagonalization approximation fails (at V (r) 2 δ /2), we must match the WKB solution onto a large-r solution. Therefore, in this work, we will choose the matching radius r M such that e − φ r M r M = max 2 δ 2 , 2 φ . For large r (r 1/ φ , so the exponential behavior dominates) we can approximate the Yukawa potential as an exponential potential V (r) ∼ V 0 e −µr , which leads to a Schrödinger equation that can be solved exactly even in the two-state case [47].
We summarize the different regimes in Figure 1. The large-r solution involves hypergeometric functions which can only be analytically matched at r M by using an asymptotic expansion. As alluded to in Table 1 (noting that µ ∼ φ ), this matching is simplest if v µ; if v µ, there is a term in the expansion that appears exponentially suppressed but can ac-quire an exponentially large prefactor. The mathematical condition that v µ is physically equivalent to ignoring the contribution from the small-r eigenstates of the 1/r potential that experience a repulsive interaction. (For further discussion, see Appendix A. The breakdown of the approximation is not inevitable in the case where the mass splitting is substantial, as discussed in Appendix E). Additionally, for v µ, we expect the s-wave term we com- Figure 1. Example of the different r-regimes and matching points for a sample parameter set ( v = 0.1, δ = 0.02, φ = 0.05), following [40]. The plot shows the exact Yukawa potential (solid black line) and the approximate potentials we employ, in their regimes of validity. In the r r M region where the eigenstates are decoupled, for r 1 the potential dominates the kinetic energy and mass splitting and is well approximated by V (r) ≈ 1/r (red dotted line), whereas for 1 r r M the WKB approximation is employed to obtain an approximate wavefunction. At r r M the WKB approximation may break down, but there V (r) ≈ V 0 e −µr (dashed blue line).
pute here to be subdominant [46] (see also Appendix E) 1 . The method developed here does not generalize straightforwardly to higher partial waves [40] -it is useful for the "resonant" regime described by [46], not the "classical" regime, which even in the elastic-scattering case demands a different approach. Consequently we focus on the part of parameter space where v µ ∼ φ . We also require that v , δ , and φ are all less than 1 in order to see substantial enhancement to the s-wave cross section, beyond the geometric cross section associated with the mass of the DM particle (not the force carrier). If v 1, the kinetic energy is large compared to the potential energy; if φ 1, the range of the interaction is short; in both cases the presence of the potential does not significantly deform the wave-function. If δ 1, the mass splitting is large compared to the Bohr potential energy (∼ α 2 m χ , leading to a suppression of virtual excitations.) Since the potential is purely off-diagonal, this suppresses the elastic scattering cross section as well.
For completeness, we include a full description of the approximate wavefunctions in Appendix A, including the WKB matching between the small-r and large-r wavefunctions. A more in-depth derivation of these wavefunctions (including more extensive discussion of the regimes of validity) can be found in [40], albeit with different boundary conditions.
In this case, we use the regular boundary conditions, φ + (0) = φ − (0) = 0. This corresponds to setting ψ(0) = 0, which is equivalent to requiring that the physical wavefunction Ψ(0) is finite at the origin. (In [40], the Sommerfeld enhancement to annihilation was extracted from irregular solutions with φ ± (0) = 0.) Additionally we will impose one of two sets of boundary conditions: the radially ingoing particles will either be purely in the ground state or purely in the excited state (since the case of a ground state particle interacting directly with an excited state particle can be treated using existing methods in the literature [46].) We include a derivation of the dimensionless transfer cross sections in Appendix B -these must be multiplied by 2 /(c 2 α 2 m 2 χ ) to obtain the physical cross sections, since we initially rescaled r by αm χ c/ .

Analytic Results
While the interaction is purely off-diagonal, it is still possible for particles to scatter "elastically", i.e. two particles incoming in the ground (excited) state may scatter off each other while remaining in the ground (excited) state. This process does not occur at tree level, in the perturbative regime (see Appendix C), but beyond the perturbative regime it is not suppressed. Figure 2 shows schematic perturbative Feynman diagrams for elastic and inelastic scattering; the non-perturbative regime requires resumming such ladder diagrams with arbitrary numbers of vector boson exchanges. Thus, below we present results for ground-to-ground and excited-to-excited elastic scattering, as well as the upscattering and downscattering rates. Schematic Feynman diagrams at lowest order for (left panel) "ground-to-ground" scattering (i.e. χ 1 χ 1 → χ 1 χ 1 ) and (right panel) "ground-to-excited" scattering (i.e. χ 1 χ 1 → χ 2 χ 2 ): swapping the fermion labels 1 ↔ 2 gives the "excited-to-excited" and "excited-to-ground" diagrams. The u-channel diagrams also contribute.
Denoting the initial and final states as "gr" for "ground" and "ex" for "excited", we obtain the following dimensionless scattering cross sections: where we have defined ∆ ≡ 2 v − 2 δ , and µ and V 0 are the defining parameters for the exponential potential V 0 e −µr , given by: The terms Γ v and Γ ∆ come from matching the WKB wavefunction onto the wavefunction for the exponential potential, and are defined by, where Γ denotes the gamma function. Finally, ϕ is a phase that comes from extending the WKB solution to the matching region. Its full form is rather complicated, but where δ φ and our other approximations hold (i.e. φ 1, v µ), the phase ϕ can be accurately approximated by ϕ = 2π/ φ (as noted in [40]). More generally, it is given by: Here r s must be chosen such that V 0 e −µrs µ 2 , 2 v , 2 δ (see Appendix A.3), but formally can (and should) be taken to −∞. To see this, note that the phase can be rewritten in the form: We can now simply take z s → ∞. Furthermore, the second integral is analytically tractable, yielding (noting that λ − is always negative): where we have chosen the convention for the branch cut such that arctanh(x) → −iπ/2 as x → ∞, for x on the positive real line. There is still one numerical integral to perform, but it is fast and stable. We emphasize again that these cross sections only hold in the regimes described in the lower half of Table I. Since we are only computing the s-wave piece of the scattering amplitude, which is angle-independent, the viscosity and transfer cross sections are trivially related to σ (see [46] for a discussion of the different cross sections and their relevance for the problem at hand).
The astute reader will notice the following salient feature of the scattering cross sections: the elastic and inelastic cross sections are the same whether the system starts in the ground state or the excited state, modulo a swap of v with ∆ (assuming that ∆ is the same in both cases, which requires the system to be above threshold.) The scattering amplitudes are identical. This reflects the identical interactions of the ground and excited states: swapping v ↔ ∆ simply corresponds to relabeling the states. The result also agrees with our intuition from quantum mechanical scattering off a 1D step potential: the transmission and reflection probabilities are the same for "downhill" and "uphill" scattering when the particle's energy is greater than the potential barrier, and the same is true in this system above the mass-splitting threshold.

Numerical Cross-Checks
Here we compare the results of our analytic approximation to the exact s-wave cross sections, obtained by numerically solving the matrix Schrödinger equation using Mathematica. Numerically solving for the scattering amplitudes proved computationally expensive and/or unstable in certain regions of parameter space, and was in all cases several orders of magnitude slower than computing the semi-analytic results, further motivating the use of these approximations. The results are shown in Figures 3-4, for two sample choices of the mass splitting parameter δ = 0.01, 0.05, and for elastic scattering in the ground and excited states (i.e. "gr→gr" and "ex→ex" respectively), upscattering ("gr→ex") and downscattering ("ex→gr"). We show only results for v > δ in this figure, since if this is not the case only ground-to-ground scattering is possible, and also restrict v < φ since otherwise we expect both that our approximation may break down and that higher partial waves will become important. We note that in all cases, notable resonances and antiresonances develop at particular values of φ .
We find that our approximations agree with the numerics to within 10% away from resonances (at resonances, minor shifts can cause huge fractional disagreement in spite of the fact that the approximation actually does capture the behavior quite well), and correctly describe the resonance positions. (Note that the minor numerical artifacts in these plots reflect the instability of the numerical calculation, and should be ignored for purposes of comparison.) Figure 5 shows the analytic calculation of the effect of a non-zero mass splitting on elastic scattering in the ground state, including for v < δ , for the same two choices of δ . For ground-state elastic scattering, the main effect of the mass splitting is to cause the positions of the resonances to shift. The effect is much more pronounced for small φ . We will understand this behavior in the following subsection, by studying analytically tractable limits of eqs. 3.1-3.4. The numerical results again agree with the analytic results in this regime (we will show an explicit comparison in the case of small v in Figure 7).    Figure 1 of [46]. The black areas to the lower right of the diagonal on the plots indicates where our approximation is no longer valid because of the repulsed eigenstates and higher partial waves (i.e. when v > φ .) The black cross-hatching on the bottom of the δ = 0.05 plot indicates where φ < 2 δ /2 which is a region where our approximation is less accurate (see Table 1.) Note that the use of "log" indicates the base 10 logarithm.

Features and Limits of the Scattering Cross Sections
We will now study various limits of the scattering cross sections using our analytic approximation. First, however, let us emphasize a point regarding our definition of v . This quantity relates to the energy of the particle relative to the ground state. A particle freely propagating in the excited state will always have v > δ . While, as noted above, the "uphill" and "downhill" scattering amplitudes are identical above the kinematic threshold given by the mass splitting, the cross section for downscattering diverges as 1/ ∆ close to threshold (except at anti-resonances) while the cross section for upscattering goes identically to zero. These different behaviors originate solely from the very different phase-space factors near threshold, and are demonstrated in Figures 3-4. For larger v , upscattering and downscattering are equally likely, since there is relatively little energetic overhead to upscattering. As well as Figures 3-4, the threshold behavior is shown for a constantφ slice in Figure 13.
These cross sections, with the same v , correspond to different physical scenarios in the context of a virialized DM halo -for instance, an excited state particle in a virialized halo will give a larger ground state velocity and thus correspond to a larger resulting v than for virialized ground state particles in the halo. We further discuss the astrophysical relevance of these scenarios in Section 4.

The Degenerate Limit
In the limit where δ → 0, analytic expressions for the scattering cross section in a repulsive or attractive potential have been previously derived [46]. In this limit, our "elastic" and "inelastic" cross sections refer to scatterings between particular linear combinations of the attracted and repulsed two-body eigenstates, and it is natural to switch to the basis of (dark) charge eigenstates accordingly (which in this limit are also mass eigenstates). In Appendix D we show in detail how this conversion is done, and find agreement with previous results [46]. In particular, we find that the resonance positions occur when ϕ = nπ and n is an integer; using the approximation ϕ ≈ 2π/ φ (valid for δ φ ), meaning the resonances occur at φ ≈ 2 π 1 n 2 . In the analysis of [46] the resonances occur at φ = 1 κ 1 n 2 where the parameter κ is chosen to be 1.6, in close agreement.
For completeness, we include the δ → 0 expansion here, as well as in Appendix D. The "elastic" scattering cross section to first order in δ (for both the ground and excited state, since they are now degenerate) becomes: (3.10) Off-resonance, in the limit as v → 0, sinh π v µ − iϕ → −i sin ϕ, and the cross section scales as 1/µ 2 . More precisely, a Taylor expansion yields: where γ is the Euler-Mascheroni constant. On-resonance, where sin ϕ = 0, the "elastic" scattering amplitude is simply 1, and the dimensionless cross section is accordingly σ "elastic" = π/ 2 v . Meanwhile, for the "inelastic" case, setting v = ∆ yields . (3.12) Off-resonance, as v → 0, this cross section approaches On-resonance, where cos(2ϕ) = 1, the "inelastic" scattering amplitude approaches 1, and again σ "inelastic" → π/ 2 v . We plot our results for the cross section in this limit in Figure 6, again finding good agreement between our full semi-analytic approximation and the numerical results for φ < 1. These results can be used more generally for elastic and inelastic scattering where δ is small relative to all the other parameters in the problem.

The Low-Velocity Limit
Now let us consider the case where v → 0, without first setting δ → 0. In this limit, scattering into the excited state is forbidden, so we will only examine elastic scattering in the ground state. Expanding σ gr→gr to first order in v yields where ψ 0 is the digamma function and γ is the Euler-Mascheroni constant. We see that the cross section does not vary with v in this low-velocity limit, with the cross section approaching the expected geometric size of ∼ π/m 2 φ , up to constant prefactors, once we convert to dimensionful parameters (making the approximation µ ∼ φ = m φ /(αm χ ), and then multiplying the dimensionless cross section by 1/(αm χ ) 2 as usual). Resonances occur when ϕ = (n + δ /2µ)π, and in the case where δ = 0 (as mentioned previously) the resonance positions are ϕ = nπ. Thus the presence of a mass splitting induces a shift to the resonance positions at velocities below the threshold, to smaller φ : Figure 7 shows the impact of the mass splitting on the resonances at low velocity for several choices of δ. The effect on the resonances is the same as found for the case of annihilation [40]. We also demonstrate in Figure 7 that our analytic approximation accurately reproduces the numerical results for lowv scattering below threshold.
Except for the shift in resonance positions, this cross section is very similar in form to (3.11); in the limit as δ → 0 (but v < δ ) they are identical, except that (3.14) has a 4γ term rather than 2γ (two extra γ's come from ψ 0 ( 1 2 )). This is a subdominant correction; generally larger contributions will arise from the cot and log terms. So we see that for elastic scattering in the ground state, the effect of the mass splitting is primarily just to shift the resonance positions; this contrasts with the case of annihilation where switching on the mass splitting can lead to a generic enhancement of the cross section by a factor of 2-4 at low velocities [40].

The Threshold ( v = δ ) Limit
Scattering amplitudes involving the excited state will be suppressed by ∆ as ∆ approaches zero from above, but the corresponding cross sections need not vanish. The case where v ≈ δ corresponds, for particles initially in the excited state, to very low physical velocities. We perform a Taylor expansion in small (but real and positive) ∆ , finding for the cross sections: where for convenience, we have defined the subdominant O(1) term which is a real quantity because the term with the hyperbolic tangent cancels out the imaginary part of the digamma function. We see that in all cases there is a potentially large enhancement corresponding to the zero-δ resonances, ϕ = nπ so cos(2φ) = 1. The cross section does not actually diverge at these pseudo-resonances, but scales as 1/(cosh(π δ /µ) − 1), and so is large when δ µ.
The upscattering cross section vanishes as ∆ → 0, as expected, as the phase space for newly-excited particles shrinks to zero. The elastic scattering cross section for the particles in the excited state scales parametrically as 1/µ 2 , except close to the resonances, where it instead scales as 1/ 2 δ if δ µ. Both these behaviors correspond to geometric cross sections, one governed by the range of the force and one by the momentum transfer associated with virtual de-excitation to the ground state. The downscattering cross section -which is perhaps most interesting for scenarios where an abundant relic population of dark matter in the halo exists in the excited state -diverges as 1/ ∆ , meaning that if v ex = α ∆ is the physical velocity of the incoming particles in the excited state, σv ex will approach a constant value at low velocities. For δ µ, the cross section scales as δ /( ∆ µ 2 ) away from the resonances, and 1/( ∆ δ ) close to the resonances. Inserting the dimensionful prefactors, the physical cross sections for downscattering and elastic scattering in the excited state have the following scaling behavior: We note that for slow-moving particles initially in the excited state, inelastic downscattering will generally dominate over elastic scattering (due to the 1/v scaling). The constant σv for downscattering implies that the argument given in [21] (which predicts a constant density core in dwarf galaxies as a direct result of a constant σv for exothermic interactions) holds even at low velocities where the perturbative approach used in that work is not valid. However, the scaling of the constant σv with the parameters of the model is quite different to the perturbative case. Note in particular that in regions of parameter space close to a resonance, large scattering cross sections can be achieved even for large m φ (provided m φ αm χ so our approximation holds), depending only on the mass splitting and the dark matter mass rather than the mediator mass.
The cross section for elastic scattering in the ground state does not have a simple behavior close to threshold, since there is nothing special about v ≈ δ from the perspective of the ground state. Setting ∆ = 0 we obtain: The term involving sinh's approaches 1 when ϕ → nπ and -1 when ϕ → (n+1)π/2, which gives rise to the characteristic resonances and anti-resonances in the low-v limit. More explicitly, we can perform a Taylor expansion in the low-velocity limit v µ (here having already set v = δ ), obtaining: σ gr→gr ≈ π µ 2 π cot(ϕ) + 2 ln (3.23) We see that this cross section has the same form as the other low-velocity and low-masssplitting limits we have studied; it is identical to the expression obtained by first taking v → 0 and then δ → 0. On dwarf galaxy scales, an elastic scattering cross section of roughly σ/m χ 0.1 cm 2 /g is required in order for dark matter self-scattering to have a significant impact on the internal structure [18]. This corresponds to particles in the core interacting once on average over the age of the universe [21], and so is likely also a necessary condition for exothermic downscattering to be relevant. We will thus use this cross section as a benchmark.
As discussed in [21], requiring a significant relic population of particles in the excited state at late times (that was not depleted by scatterings in the early universe) requires m χ at the MeV scale or lighter. However, the excited state might be populated non-thermally, in which case much heavier DM masses might also be viable.
For the non-degeneracy of the excited state to have a significant impact on scattering in dwarf galaxies, the mass splitting should be significant compared to the typical kinetic energy of the dark matter particles. Taking the typical velocity in dwarf galaxies to be 10 km/s∼ 3 × 10 −5 c [49], this implies δ 10 −9 m χ in order to see differences from purely elastic scattering. Our requirements that δ 1, φ 1 impose that δ α 2 m χ and αm φ α 2 m χ 2 . So in order for our approximation to be valid and the mass splitting to be interestingly large, we will focus on the range α 10 −4 (or higher for larger δ: α δ/m χ ), which will also guarantee v = v/α 1 as required. For a vector mediator, this is in agreement with broad expectations from the Standard Model if the coupling is not fine-tuned to be small.
In general, we will treat α as a free parameter within the range 10 −4 α 1, since the constraints on it are rather model-dependent. In particular, we do not impose the constraint that annihilations of the DM to the force carrier should generate the correct relic density. The self-interacting DM could be non-thermal or a sub-dominant component of the total DM if α is higher than the value expected for a thermal relic, or annihilation channels not involved in the self-scattering could prevent the overclosure of the universe if α is too small. However, for calibration, [41] found typical values for α (yielding the correct relic density) of a few times 10 −2 for TeV-scale DM in a model with the same potential as the one employed in this work, and the annihilation of the dark matter to the force carriers has a cross section scaling as α 2 /m 2 χ at freezeout, so lighter DM will generally imply a smaller value of α if the DM is indeed a thermal relic.
There are few model-independent constraints on the mediator mass m φ ; the coupling of the force carrier to Standard Model particles is independent of its role here of mediating dark matter scattering. For significant scattering we require that m φ αm χ , and for swave scattering to dominate and our approximation to be valid we will generally require that v φ in dwarf galaxies, i.e. m φ m χ v ∼ 3 × 10 −5 m χ . When we consider the exothermic scenario with a significant population of particles initially in the excited state, their scatterings have v ∼ δ in our notation, assuming the kinetic energy of the excited-state particles (limited by the escape velocity of the dwarf) is small compared to the mass splitting. Thus for this scenario we will also require δ φ , i.e. 2δ/α 2 m χ m 2 φ /α 2 m 2 χ ⇒ m φ δm χ . The requirement that m φ αm χ means that the higher α is above its lower bound of δ/m χ , the more valid parameter space there will be for m φ (although raising m φ above the upper bound simply means we should use the Born approximation, as in Appendix C, rather than the approximations derived in this article.)

Results
With the above reasoning as our guide, we performed several cuts through parameter space to estimate phenomenologically interesting regions in DM mass, mediator mass, mass-splitting, and coupling. We selected a few fiducial velocity scales which correspond to the virial velocities of dwarf galaxies, galaxies the size of the Milky Way (MW), and clusters. Practically speaking, the DM scattering rate within a halo is not determined by one fixed velocity, but rather some distribution of velocities as DM particles pass through different parts of the halo. In order to understand these dynamics in detail, one must perform numerical simulations which lie outside the scope of this work. However, the simple estimates presented here can serve as benchmarks for simulations of the dynamical behavior of inelastically scattering DM in an actual halo.
The results of our parameter sweeps for four different scattering scenarios (ground to ground, ground to excited, excited to excited, and excited to ground) can be found in Figures  8-11, respectively. All plots depict σ T /m χ in nine slices of the m χ vs. m φ plane with α and δ held fixed (note we label the cross section as the transfer cross section σ T for comparison to the literature, but in our s-wave approximation it is identical to the cross section we have discussed so far; see Appendix B.3.3 for a discussion). Each slice has different values of α and δ, with α decreasing from left to right and δ increasing from top to bottom (with each subplot labelled accordingly.) We also include three velocity scales: 10 km/s for dwarf galaxies, 200 km/s for MW-sized galaxies, and 1000 km/s for clusters. In terms of our model, these velocities correspond to v for situations where the particles are initially in the ground state and ∆ for situations where the particles are initially in the excited state. Additional features of the plots are listed as follows: • Black regions show where our approximation breaks down. Since we know that our approximations will give wrong or misleading cross sections in these regions, we cover them entirely. Black, triangular regions in the lower right corner show when φ > 1. Black, horizontal strips across the bottom of the plots (when present) show when δ > 1. Black, vertical strips across the left side of the plots (when present) show when 2 δ /2 > φ . Black regions in the upper left show when higher partial waves and repulsed eigenstates can no longer be ignored, which happens when v > φ . Note that we chose to black out the region with v corresponding to dwarf scales, e.g. the blacked out regions have φ < 3 × 10 −5 /α for scenarios where the particles are incoming in the ground state. Since this criterion is velocity-dependent, we only plot cross sections for DM in the MW and clusters when our approximations are valid at those velocity scales (rather than blacking out regions where the approximation is not valid for those scales.) • Yellow regions show 0.1 cm 2 /g < σ T /m χ < 1 cm 2 /g (light), 1 cm 2 /g < σ T /m χ < 10 cm 2 /g (medium), and 10 cm 2 /g < σ T /m χ < 100 cm 2 /g (dark) on dwarf velocity scales with v ∼ 10 km/s.
• Blue contours represent σ T /m χ = 0.1 cm 2 /g (light) and σ T /m χ = 1 cm 2 /g (dark) on MW velocity scales with v ∼ 200 km/s. We emphasize that the contours are only plotted when v < φ on these scales, which ensures that the repulsed eigenstates and higher partial waves are safely neglected.
• Red contours represent σ T /m χ = 0.1 cm 2 /g (light) and σ T /m χ = 1 cm 2 /g (dark) on cluster velocity scales with v ∼ 1000 km/s. We emphasize that the contours are only plotted when v < φ on these scales, which ensures that the repulsed eigenstates and higher partial waves are safely neglected. The strongest robust constraints on DM selfinteraction come from observations of merging clusters, with the upper bound being around σ T /m χ ∼ 1 cm 2 /g (e.g. [50]). Thus, parameter space to the lower left of the dark red contours is effectively ruled out for our model at the selected masses and couplings.
• Teal contours represent the Born approximation with σ T /m χ = 0.1 cm 2 /g (light) and σ T /m χ = 1 cm 2 /g (dark). We present a full derivation of the Born cross sections in C but report them here. For elastic scattering in the limit of slow-moving initial particles, the cross sections for ground-state elastic scattering and excited state elastic scattering are respectively (4.1) The Born approximation for the elastic cases is shown in dot-dashed lines (we wish to emphasize that we have taken a low-velocity limit in computing the Born approximation, so the result does not actually depend on the velocity of the incoming particles; it is most likely to be valid for dwarfs where the virial velocities are very low).
Meanwhile, for the inelastic cases, the upscattering cross section is and the downscattering cross section is For the inelastic cases, the Born approximation is shown in dotted lines for MW velocities and dashed lines for cluster velocities.   Figure 8, except for inelastic scattering from the ground state to the excited state (upscattering). Upscattering in dwarf-sized halos, even where kinematically allowed, was never significant for the parameters we sampled (i.e. σ/m χ 0.01 cm 2 /g for all of the parameter space.) Note that the horizontal cutoffs of the contours come from the mass splitting threshold, as upscattering does not occur below threshold.

Discussion
The results in Figures 8-11 show that in the parameter region where our approximation holds, parameters below the dark red line are definitively ruled out by cluster mergers while yellow regions are favored because they potentially solve the small-scale structure anomalies. Interestingly, tentative evidence from cluster mergers suggests that the dark matter selfinteraction cross section may be nonzero, and the favored value of σ/m χ is 0.8 cm 2 /g at around 1σ [51]. Thus, the discovery of future mergers may further constrain self-interacting dark matter or even possibly pick out a strongly-preferred interaction cross section at that velocity scale. The resonances occur in thin bands of roughly constant φ from the lower left corner of the m χ vs. m φ plane up to the upper right corner. Decreasing α spaces the resonances further apart and broadens them due to the weak dependence of the ratio δ /µ (which appears in our expressions for predicted resonance positions) on α. Increasing δ shifts the resonance positions slightly. For elastic scattering in the ground state, increasing δ elongates the resonance bands and makes them slightly more curved; the most pronounced curvature develops at small m φ , as expected from the discussion in Section 3.3.2. Increasing δ has the opposite effect on elastic scattering in the excited state, in that it dampens the resonances (as discussed in Section 3.3.3) and so reduces their effect on the favored regions. Our approximations match smoothly onto the perturbative expressions as φ → 1, as expected.
We see that for these mass splittings, in the region of greatest astrophysical interest the ground-ground elastic scattering cross section is generally similar to the case where the DM states are degenerate [46]. Additionally, the favored regions are quite comparable for elastic scattering and downscattering, so in the context of this model it is possible for both effects to simultaneously contribute to the dynamics of DM within halos and help alleviate small-scale structure issues. This is a consequence of our choice of mass splittings fairly similar to the kinetic energy of virialized DM particles, since the (non-resonant) ratio of the scattering cross sections at low velocity scales as δ/(m χ v 2 ) (see Eq. 3.19). However, for much higher mass splittings the scattering will most likely be in the classical high-velocity regime for most of the parameter space of interest, unless the mediator mass is also raised. However, as expected from our earlier discussion, at small m φ the resonances occur for different parameters for downscattering, compared to elastic scattering in the ground state. Consequently, there are regions of parameter space where the ground-state elastic scattering cross section is large and the downscattering rate small, and vice versa.
Interestingly, for the case of upscattering there are no regions in the sampled parameter space where upscattering was significant in dwarf halos (i.e. there were no regions where upscattering exceeded 0.1 cm 2 /g.) From the perspective of avoiding dark matter "cooling" in dwarf halos (which could potentially worsen the core-cusp problem, etc.) the ease of suppressing upscattering is an appealing feature. However, there are substantial regions of parameter space where upscattering is significant for the MW and for clusters. It is possible that upscattering could contribute to steepening of the density profile in the central parts of such large halos.

Regimes of Interest for XrayDM
There has been a great deal of recent interest in the detection of an apparent ∼ 3.5 spectral keV line in radio observations of galaxies and galaxy clusters, as a potential signal from dark matter [38,39]. However, there appears to be some tension between the interpretation of this line as originating from the decays of keV-scale dark matter, and its non-detection in the Virgo Figure 12. Contours satisfying (4.4) for the characteristic velocity scale of clusters, v = 1000 km/s. These lines serve as illustrative benchmarks -in reality, there will of course be some distribution of velocities in any cluster. The region where our approximation breaks down because v > φ is masked out, but it is likely that this part of parameter space can also furnish appropriate cross sections. The horizontal dashed line indicates the threshold for upscattering to be kinematically allowed, v > δ .
cluster. It has been proposed that the signal could instead originate from the upscattering of (weak-scale) dark matter to an excited state ∼ 3.5 keV heavier than the ground state, followed by decay back to the ground state with emission of a photon [33]. This scenario was termed "XrayDM". Since the upscattering process involves two particles, the rate of excitations (and hence decays) scales with the density squared rather than the density (as would also be the case for annihilation, e.g. [52]), and also depends on the typical velocity of the DM particles: this can modify the relative strength of the signal in different regions 3 .
The example model employed in the XrayDM scenario of [33] is exactly the simple model studied here: accordingly, we can now use our approximation to calculate which regions of parameter space can give rise to a sufficiently large cross section to explain the 3.5 keV line. Due to the large virial velocities of clusters, for small m φ and/or large m χ our approximation becomes invalid (as higher partial waves become important), but there is a significant region of interesting parameter space where the s-wave contribution dominates, as shown in Figure  12. Here we have imposed the criterion that: (4.4) following [33], with v = 1000 km/s and δ = 3.5 keV. As previously, we only show the result for a single velocity, rather than integrating over a distribution, since our purpose here is to provide estimates rather than a detailed analysis of the allowed parameter space. Within the regime of validity of our approximation, we find that DM masses in the range of a few GeV to a few tens of GeV, moderately large values of α, and mediator masses in the range of 10 MeV -1 GeV can naturally produce the required flux. We expect that similar DM masses and lower mediator masses will also provide viable explanations, but higher partial waves become important in these cases. We do not show results for smaller α as in that case, for interesting cross sections, the perturbative regime transitions directly to the classical regime (where high partial waves are important) without an intermediate resonant regime -this does not mean no parameter space is open for smaller α, just that it does not require the use of our results.
A similar analysis could be performed for the original "XDM" scenario [31]; we do not carry it out here because in that case the mass splitting is quite large (∼ 1 MeV), meaning that only particles in the high-velocity tail of the velocity distribution (for DM in the Milky Way) are typically able to upscatter. The viability of this scenario thus depends critically on the details of the velocity distribution, and also typically requires a more complex model where both particles in a collision do not need to excite simultaneously (to reduce the energy requirement for upscattering).

Conclusions
We have derived a semi-analytic approximation for the scattering cross sections for dark matter interacting via an off-diagonal dark Yukawa potential. This physical model demonstrates several new features relative to a system with purely elastic scattering, and we have presented simple analytic forms for the upscattering and downscattering rates in the nonperturbative resonant regime. Within the regime of validity of our approximations, we readily obtain cross sections that appear large enough to modify the inner regions of dwarf halos or give rise to interesting signatures in indirect DM searches.
Our approximations are valid when v m φ /m χ α and (for particles initially in the excited state) δ m 2 φ /m χ α 2 m χ , where m χ is the DM mass and m φ the mass of the mediator. Within this regime, the scattering cross section exhibits resonant enhancement at particular values of the parameter m φ /(αm χ ). For particles initially in the excited state, or scatterings of ground-state particles with energies above the excitation threshold, the resonance positions are very similar to those in the case with no mass splitting: however, for scatterings in the ground state below the excitation threshold, the resonances are shifted by the presence of the mass splitting when δ becomes comparable to m 2 φ /m χ . Consequently, the relationship between the elastic and inelastic scattering cross sections can be quite complex. Away from the resonances, we recover the expected geometric cross sections at low velocities: σ ∝ 1/m 2 φ for elastic scattering, σv ∝ δ/m χ (1/m 2 φ ) for downscattering. Our approximations match smoothly onto the corresponding perturbative expressions when m φ ≈ αm χ .
We hope that in future work, the incorporation of these semi-analytic scattering cross sections into numerical simulations will allow the first detailed studies of halos containing inelastically scattering dark matter.

A.1 Approximate Small-r Wavefunctions
For small r, we can approximate the Yukawa potential as 1/r, which therefore dominates the eigenvalues at small r so λ ± ≈ ±1/r. The s-wave solutions to the rescaled Schrödinger equation can be expressed in terms of Bessel functions: where A + and A − are the coefficients of the repulsed and attracted eigenstates, respectively. Moving radially outward (but still within the regime of validity for the V ∼ 1/r approximation), the large-r asymptotics of the Bessel functions give φ + (r) = 1 (A.2) If the particle velocity is high enough (above threshold) such that there exists a radius r * in this regime where V (r * ) = v 2 v − 2 δ then λ + (r * ) = 0 and we must perform a WKB approximation about the turning point. Linearizing the potential and matching the wavefunctions on either side using the connection formulae yields: In either case (above or below threshold), the small-r wavefunctions are in a form that will match smoothly onto the WKB wavefunctions.

A.2 Approximate Large-r Wavefunctions
At large r, we can approximate the Yukawa potential as a purely exponential potential of the form V 0 e µr . We impose conditions on V 0 and µ by requiring that the exponential potential mimic the Yukawa for r > r M , where r M is the matching radius which we have chosen such that V (r M ) = max from both potentials match to first order in the coupling constant α (following [44]). The parameters µ and V 0 are therefore given by The wavefunctions for an exponential potential can be solved for exactly in terms of 0 F 3 hypergeometric functions [47] as follows: The wavefunctions are expressed in terms of four linearly-independent solutions, corresponding to ingoing or outgoing particles in the ground or excited states. In particular, • the C 1 term represents an ingoing wave in the ground state, • the C 2 term represents an outgoing wave in the ground state, • the C 3 term represents an ingoing wave in the excited state, • the C 4 term represents an outgoing wave in the excited state.

A.3 WKB Approximation for the Intermediate-r Wavefunctions
To match the large-r wavefunctions to the small-r wavefunctions, one can use the WKB approximation to propagate the known wavefunctions of the exponential potential into the transition region. We write the large-r WKB solutions as whereλ ± are the eigenvalues of the matrix Schrödinger equation with an exponential potential rather than a Yukawa. In order to match the WKB solution with the exact solution for the large-r exponential potential, we define the following convenient quantities: Then, deriving expressions for the WKB coefficients E ± and F ± to match onto the exponential wavefunctions is a matter of using the asymptotic behavior of the 0 F 3 hypergeometric functions in the r → −∞ limit (see [40] for details), and using the WKB approximation again to propagate these asymptotic solutions into the matching region. We find that where r s is some (possibly negative) radius chosen such that (a) V 0 e −µr µ 2 (required for the asymptotic expansion to be valid) and (b) the potential dominates the kinetic and masssplitting terms, i.e. V 0 e −µrs 2 v , 2 δ . There is no reason not to choose r s arbitrarily large and negative, however, and as discussed in the main text, we generally take r s → −∞. (Of course, this does not correspond to a physical region of real space, but this matching is simply a mathematical trick to translate the hypergeometric functions into a form that facilitates matching to the previously derived WKB solutions.) Note that this procedure sets E + = 0. This originates from the neglect of a particular (exponentially suppressed) term in the asymptotic expansion for the hypergeometric functions. However, this term can gain an exponentially large prefactor when v µ, as discussed in [40]; we shall explore the consequences of neglecting this term further in Appendix B.
In order to match the WKB solution with the small-r solution below threshold, we equate (A.7) with (A.2), which gives and similarly, above threshold, equating (A.7) with (A.3) gives (A.11) where r † is the radius above threshold at which the eigenvalue of the repulsed eigenstate of the exponential potential passes through zero, defined by V 0 e −µr † = v 2 v − 2 δ . We can then equate the coefficients from matching the WKB solution onto the small-r solution with the coefficients from matching the WKB solution with the large-r solutions. After imposing appropriate boundary conditions, we can then determine the coefficients for the various physical solutions to the Schrödinger equation. For instance, we can calculate the coefficients for the repulsed and attracted eigenstates for the small-r solution. We can also determine the coefficients for the ingoing and outgoing spherical waves in the ground or excited states for the large-r solution.

B.1 Matching the Wavefunctions Using the Boundary Conditions
As mentioned in Section 3, we require that φ + (0) = φ − (0) = 0, since φ(r) is the radial wavefunction rescaled by r and the wavefunctions should be regular at the origin. For utility, we then define the following useful phases: where r s is some radius chosen such that V 0 e −µrs 2 v , 2 δ , as in (A.9). Then, for the below-threshold case, equating (A.9) with (A.10) at the matching radius, r M gives where we have defined ∆ ≡ 2 v − 2 δ and Γ v and Γ ∆ are defined in (A.8). Similarly, above threshold, equating (A.9) with (A.11) gives Both above and below threshold, the E + equation gives us A + = 0, which is akin to neglecting the contribution from the repulsed eigenstate at small radii. This indicates that repulsive scattering occurs most significantly at large r, where the exponential part of the Yukawa dominates the behavior of the wavefunction. This makes intuitive sense since we are interested in the low-velocity limit, which means that incoming particles must climb up to the classically disallowed region of a repulsive potential in order to even reach the small-r region. Recall that the result E + = 0 was a consequence of assuming v µ: when this assumption is not satisfied, we cannot expect A + = 0 on physical grounds.
Since A + = 0, the F + equation gives both above and below threshold.

B.2 General Scattering Amplitudes
Once we have imposed the relevant boundary conditions on the large-r wavefunctions, we can extract the scattering amplitudes by reading off the coefficients of the outgoing solutions.
Since the hypergeometric functions asymptote to 1 as r → ∞, the large-r ground and excited wavefunctions approach ingoing and outgoing spherical waves: More generally, consider wavefunctions for nondegenerate states X and Y given by where the A terms represent the unscattered wavefunction in state X, the B term represents the elastically scattered wavefunction in state X, and the C term represents the inelastically scattered wavefunction in state Y (hence the wavenumber k as distinct from k.) In this case, conservation of probability current dictates that k |A| 2 = k |A + B| 2 + k |C| 2 . We can reformulate these wavefunctions (recall that all wavefunctions used in this paper are rescaled by r) in the context of a 3-dimensional scattering problem as an ingoing cylindrical wave and a scattered outgoing spherical wave: where the second equality follows because for the s-wave, e ikz → j 0 (kr)P 0 (cos θ) = sin kr kr . (B.8) In general to get a differential cross section, we relate the incident probability flux through an area to the scattered outgoing probability flux through a solid angle: Equating the top row of (B.7) with ψ X from (B.6) gives (B.10) Similarly, equating the bottom row of (B.7) with ψ Y from (B.6) gives If we apply the analogy to our wavefunctions for the case where we begin purely in the ground state (which corresponds to setting C 3 to zero), then the elastic scattering cross section is and the inelastic scattering cross section is Similarly, for the case where we begin purely in the excited state (which corresponds to setting C 1 to zero), then the elastic scattering cross section is and the inelastic scattering cross section is We will impose boundary conditions such that the ingoing wave is purely in the ground state, which implies that C 3 = 0. We are free to set C 1 = 1 up to some overall normalization.
Dividing the E − equation by the F − equation of (B.2) yields and combining this with (B.4) gives So by (B.12) the elastic scattering cross section is and by (B.13), the inelastic scattering cross section is .

B.3.2 Incoming in the Excited State
We will impose boundary conditions such that the ingoing wave is purely in the excited state, which implies that C 1 = 0. We are free to set C 3 = 1 up to some overall normalization.
Dividing the E − equation by the F − equation of (B.2) yields 20) and combining this with (B.4) gives So by (B.14), the elastic scattering cross section is and by (B.15), the inelastic scattering cross section is . (B.23)

B.3.3 Relation to the Transfer Cross Section
When considering the effects of DM scattering on structure formation in general, the total cross section may acquire an unphysical forward divergence. Generally, the literature has instead employed the transfer cross section σ T , which determines the longitudinal momentum transfer: (B.24) Since our differential cross sections are angle-independent, we can pull those out of the integral. Since the cos θ term is orthogonal to the sin θ term in the dΩ Jacobian, the remaining integral just gives which is the same as the cross section that we computed. (Note that [46] argues for the use of the viscosity cross section instead, but computes σ T in order to make contact with the literature -we follow their approach, albeit we reiterate that for the s-wave component of the amplitude, the distinction is a trivial one.)

C The Perturbative Regime
The Born approximation can be straightforwardly applied to this multi-state system, by writing the three-dimensional Schrödinger equation in the form: is the full, three-dimensional physical wavefunction for the twobody state; for the spherically symmetric potential and s-wave scattering we consider, it is related to the wavefunction ψ(r) by ψ(r) = rΨ( r). Ψ 0 ( r) is the equivalent wavefunction for the unperturbed system with no potential, satisfying ∇ 2 Ψ 0 ( r) = − k 2 0 0 k 2 · Ψ 0 ( r). (Note, since Ψ( r) describes a two-body state, we have replaced the mass in the Schrödinger equation with the reduced mass m χ /2.)V ( r) is the matrix potential given in Eq. 2.4, but now with the diagonal 2δc 2 term omitted (we instead incorporate the effect of this term in the difference between k and k). Here we use units where c = = 1.
The first-order Born contribution to the scattering amplitude is solely inelastic, since the matrixV is purely off-diagonal. For an initial state consisting of a plane wave in a single state (purely for notational purposes, we here choose the upper row of the matrix to correspond to the initial state), we obtain: Writing k = kr, k = k r, p = kẑ, p = k ẑ, and making the long-distance approximation e ik| r− r 0 | /| r − r 0 | ≈ (e ikr /r)e −i k· r 0 , we obtain: We can read off the scattering amplitudes as (in the notation of Appendix B): where θ describes the angle betweenr andẑ, and as in Appendix B we have: Performing the angular integral yields: This result applies to both particles initially in the ground or excited states, with the proper choices of k and k . For upscattering, k = m χ v and k = m 2 where v is the speed of a single particle in the center-of-mass frame in the initial state. If we use the convention in the body of the text where v always refers to the speed of the particle in the ground state, then for downscattering these choices of k and k are simply reversed. (If instead we take v to be the speed of the particle in the initial state, for downscattering we have k = m χ v and k = m 2 χ v 2 + 2δm χ .) Thus we obtain the upscattering cross section, and the downscattering cross section (taking v to be the speed of the ground-state particle): The corresponding dimensionless cross sections are respectively: The lowest-order contribution to elastic scattering comes via the second term in the Born series. For the initial condition above, the contribution to the wavefunction is: For scattering where the initial particles are in the excited state, k = 2δm χ for k = 0, and the amplitude is given by: (C.14) For scattering where the initial particles are in the ground state and k = 0, k = i 2δm χ , and we obtain: Thus the cross sections for ground-state elastic scattering and excited-state elastic scattering, in the limit of slow-moving initial particles, are respectively, (C. 16) The corresponding dimensionless cross sections are: Note that these second-order cross sections become comparable to (or larger than) the first-order inelastic scattering cross sections for φ 1, signaling the breakdown of the Born approximation and the need to transition to the resonant regime described in the main text.

D Comparing our δ → 0 Limit with Previous Results
An analytic approximate form for the phase shift due to elastic scattering, for both attractive and repulsive potentials, has previously been presented in the literature [46]. We show here how to recover the analogous result in our approximation.
For the − state, let the scattering solution φ − (r) (for the coefficient of the ψ − eigenvector) have the asymptotic form φ − (r) = e i(kr+2δ − ) + e −ikr . The phase shift δ − characterizes the scattering amplitude, which is given by 1 − e 2iδ − 2 . Likewise, for the + state, let the phase shift be δ + .
Since the differential equation is linear, any linear combination of these solutions (Aφ − (r)ψ − + Bφ + (r)ψ + ) is also a solution. In particular, if we set A = B = 1/ These correspond to the cases we studied above, where the particles are initially purely in the ground or excited states. So we see that by calculating the phase shifts for these initial conditions (given by the A, B and C coefficients in Eq. B.10 and Eq. B.11), we can recover the values for δ − and δ + , and vice versa. Our cross sections are given by: (D.2) and in this case, since δ = 0, swapping the identifications of "ground" and "excited" states has no effect. Note that the sum of these cross sections gives σ tot = π k 2 1 − e 2iδ − 2 + 1 − e 2iδ + 2 = σ − + σ + , as it must -the total scattering rate cannot depend on the choice of basis.
Off-resonance, as v → 0, this probability approaches σ inelastic → π π cot ϕ µ 2 . (D.11) On-resonance, where cos(2ϕ) = 1, the inelastic scattering amplitude simply approaches 1. We see that these would agree precisely with the approximate forms of the cross sections derived from the results of [46] if we made the replacements: These replacements are parametrically correctµ ∼ φ ∼ 1/c up to O(1) factors, likewise V 0 ∼ µ up to O(1) corrections. The results are most sensitive to the identification ϕ → π √ c, since this sets the resonance positions: taking c = 1/(κ φ ), and our approximate expression ϕ ∼ 2π/ φ , we see that they agree exactly if κ = π/2 ≈ 1.57. The value of κ = 1.6 chosen by [46] therefore leads to percent-level agreement in the resonance positions.
Perfect agreement between the two analyses should not be expected, since they use different potentials (albeit with similar properties), but our approach agrees both qualitatively and quantitatively with the results of [46] in the region of parameter space where they can both be used.   Figure 13 shows a particular slice through the parameter space of Figure 3, for inelastic scattering, but extending to v > φ . As expected, our analytic approximation breaks down in this regime (note that µ and φ are generally equal up to a O(1) factor). Conveniently, v φ is also precisely the condition for s-wave scattering to dominate over the higher partial waves. Consequently, while a more careful treatment of the matching between the WKB and large-r regimes (see Figure 1) might allow extension of our approximation for the s-wave to the region with v φ , at that point it would be necessary to include the higher partial waves as well.

E The High-Velocity Regime
This can be easily seen by comparing the relevant length scales: for the th partial wave, the vacuum solution is proportional to the Bessel function j ( v r), which peaks when ∼ v r, i.e. r ∼ / v . In order for scattering of the th partial wave to be significant, this peak must lie within the range of the potential, i.e. r 1/ φ , and so we must have v / φ . If v / φ 1, then only the s-wave term can penetrate the potential far enough to experience significant scattering.
In the case where δ is non-zero, this argument still holds -for particles in the excited state, the asymptotic wave function is now j ( 2 v − 2 δ ), but since 2 v − 2 δ < v , requiring v φ is certainly sufficient to ensure that the potential cuts off at smaller r than the peak of the higher-wavefunctions. Since whenever particles in the excited state are present and the scattering rate is significant, their downscatterings will populate the ground state, we will generally consider v φ to be both a necessary and sufficient condition for our approximate solution to be useful.

E.2 The Adiabatic Regime
However, there is a regime where v φ but our approximate solution remains valid, although the s-wave does not generally dominate scattering in this part of parameter space, and so we caution that our s-wave result should not be used as a proxy for the total scattering cross section. However, it can be used as a lower bound.
As discussed in Appendix B, our method essentially neglects repulsive scattering at small distances, which is valid when the range of the potential is relatively short and so the scattering wavefunction for the repulsed eigenstate is peaked outside its range (note this is the same reason we can ignore the higher partial waves in this regime).
There is another regime where this approximation is valid, for a different reason. Suppose the system starts with both dark matter particles in the ground state (i.e. the state of lowest energy). At short distances, the lowest-energy eigenstate is the one that experiences an attractive potential (corresponding to the +− two-body state at high energies). If the transition from long distances to short distances is adiabatic -i.e. this transition occurs slowly relative to the scale associated with the splitting between the eigenstates -then particles in the lowest-energy eigenstate at long distances will find themselves entirely in the attracted eigenstate at short distances, in analogy to the adiabatic theorem, and ignoring the repulsed eigenstate will be valid because it will simply never be populated.
The splitting between the eigenstates corresponds to an energy scale of 2 δ in our dimensionless coordinates, and hence to a time scale of 1/ 2 δ ; the corresponding distance scale, for an inward-moving wavepacket, would be ∼ v / 2 δ . The rotation of the eigenstates (from mass eigenstates to gauge eigenstates of the high-energy potential), as described in Section 2.1, occurs when V (r) = e − φ r /r becomes comparable to 2 δ /2. If the cause of this transition is the exponential cutoff, i.e. r ∼ 1/ φ , then the transition occurs over a radius ∆r ∼ 1/ φ ; if φ 2 δ then it occurs when V (r) ∼ 1/r and over a range ∆r ∼ 1/ 2 δ . So the criterion for adiabaticity is v / 2 δ 1/ φ if φ 2 δ , or v / 2 δ 1/ 2 δ otherwise. In the first case, the transition is adiabatic for v φ 2 δ ; in the second case adiabaticity always holds for v 1 (however, note that φ 2 δ is a regime where the approximations we use are known to be less accurate [40]), which by the condition φ 2 δ implies φ v 2 δ . We can summarize this by saying adiabaticity holds if and only if φ v 2 δ , provided our other assumptions hold (that is, v , δ , φ 1).
This mechanism is also responsible for the enhancement in the annihilation rate noted in [40] for the case with a mass splitting, compared to the case where the mass splitting is  Figure 14. A scan through v with φ = 0.03 and δ = 0.04. This demonstrates the shift from the transition to small-r being adiabatic vs. nonadiabatic. We can see the breakdown near 2 δ ∼ v φ , which happens near v ∼ 0.06. negligible relative to the kinetic energy and can be ignored; under adiabatic conditions, the presence of the mass splitting causes particles initially in the ground state to transition into a purely attracted state, rather than a equal linear combination of attracted and repulsed states.
This argument cannot be applied to scattering from the excited state or into the excited state, as if the excited state is populated then this implies the repulsed eigenstate will also be populated and cannot be ignored. But provided we are only interested in elastic scattering from the ground state (i.e. for the below-threshold case v δ ), we expect our results to be accurate (for the s-wave) even when v φ , in the event that φ v