WIMPonium and Boost Factors for Indirect Dark Matter Detection

We argue that WIMP dark matter can annihilate via long-lived"WIMPonium"bound states in reasonable particle physics models of dark matter (DM). WIMPonium bound states can occur at or near threshold leading to substantial enhancements in the DM annihilation rate, closely related to the Sommerfeld effect. Large"boost factor"amplifications in the annihilation rate can thus occur without large density enhancements, possibly preferring colder less dense objects such as dwarf galaxies as locations for indirect DM searches. The radiative capture to and transitions among the WIMPonium states generically lead to a rich energy spectrum of annihilation products, with many distinct lines possible in the case of 2-body decays to $\gamma\gamma$ or $\gamma Z$ final states. The existence of multiple radiative capture modes further enhances the total annihilation rate, and the detection of the lines would give direct over-determined information on the nature and self-interactions of the DM particles.


I. INTRODUCTION
One of the most promising ways in which to probe the nature of weakly-interacting-massive-particle (WIMP) dark matter is provided by indirect detection experiments in which the annihilation products of dark matter in astrophysical contexts are observed. Indirect detection experiments such as AMS [1], ATIC [2], EGRET [3], GLAST [4], HEAT [5], HESS [6], INTEGRAL [7], PAMELA [8], and VERITAS [9] look for signals ranging across energetic gamma rays, positron excesses, and antiproton fluxes, and hints of deviations from background expectations now abound. A common feature in the interpretation of these experiments is that the WIMP annihilation proceeds via simple, almost free-particle annihilation leading to a rate that depends on the WIMP relative velocity in only a very simple, essentially structureless way, and moreover, is directly related to the crosssection that led to the DM density at freeze-out. This assumption then implies that the primary astrophysical quantity determining the fluxes from DM annihilation is the local DM squared-density, ρ(x) 2 , with the velocity distribution of the DM being essentially irrelevant. In addition the assumption of almost free-particle annihilation gives a relatively simple resulting energy spectrum for the annihilation products (though of course features such as steep falls at kinematic thresholds are generic). The dependence of the annihilation fluxes on ρ(x) 2 has been widely employed in the suggested interpretations of the various experimental anomalies as, very frequently, a so-called "boost-factor" in the annihilation rate is required to match the observed flux, and this is assumed to come from local over-densities of DM.
In this paper we point out that annihilation of WIMP dark matter via intermediate long-lived "WIMPonium" bound states, Ω (n, ) , is possible in many particle physics models of DM [10] (see [11] for another recent discussion of WIMPonium). As we argue below, the WIMPonium bound states can occur at or near threshold in which case they lead to potentially very substantial (factors of 10 3 to 10 5 ) enhancements in the DM annihilation rate, closely related to the well-known Sommerfeld non-perturbative enhancement [12] that has been applied to freeze-out calculations and indirect dark matter signals [11,[13][14][15][16]. The existence of this ρ 2 -independent dynamicallyinduced "boost factor" has important implications: First large amplifications can occur in low-velocity-dispersion systems without a large (eg, cusp-like) density enhancement 1 , possibly preferring colder less dense objects such as dwarf galaxies as more promising places to search for DM signals compared to the higher density but higher velocity galactic center 2 . Second, a knowledge of the full phase space density of the DM is necessary to reliably compute the DM annihilation rate, which strengthens further the case for detailed realistic simulations of the DM distribution in our galaxy.
In addition, the deeply-bound WIMPonium spectrum is often quite involved, and the radiative capture to and transitions among the various states generically lead to a rich energy spectrum of annihilation products, with many distinct lines (in the case of decays to γγ or γZ 1 Furthermore cusp-like dark matter distributions seem not to be favoured at present, see for example [17] final states, similar to that found in WIMP annihilations [19]) possible. Although the observability, or otherwise of these distinct lines depends on the DM model being studied, and especially the resolution of the detector, the existence of all the various radiative capture modes further enhances the total annihilation rate, and the detection of even just a few of the lines would give direct information on the interactions of the DM particles.

II. THRESHOLD BOUND STATE BASICS
In order to demonstrate the important features of WIMPonium bound states it is useful to examine a simple, almost model independent, example. We introduce two complex scalar fields s and n with the following interactions and masses where λ is a dimensionless coupling and µ, m n and m s are mass dimension one parameters. We assume that ≡ m n /m s 1 (typically we will take m n ∼ m Z , while 500 GeV < ∼ m s < ∼ 30 TeV). A discrete symmetry is imposed forcing s particles to appear in the Lagrangian in pairs so the s scalar will be the stable dark matter WIMP.
As we have indicated in the introduction we are interested in two related phenomena: 1) The scattering of slow moving particles near bound state thresholds leading to amplification of the direct annihilation rate; 2) The radiative capture to and decay of deeply bound states and transitions between different bound states.
Since we are interested in low-velocity processes we can proceed by solving the Schrödinger equation for the two-s-particle system with the Yukawa potential V = −λ 2 e −mnr /8πr that follows from the exchange of n-scalars. Here λ ≡ λ µ/m s . Semi-classical considerations show that the number of bound states N of given orbital angular momentum satisfies [20] (2 + 1)N < 2M r r|V (r)|dr, where M r is the reduced mass of the two particle system, so, for example, the number of Swave bound states in our case satisfies N 0 < ms mn α, where α ≡ λ 2 /8π. A more precise condition for there to be at least one bound state follows from numerical methods giving [21] α ≥ 0.84m n /m s . Thus if we want to have a rich structure of energy levels our dark matter particles will either need to have large couplings with the n scalar or need to have m s /m n 1, or both. We emphasize, however, that the most important phenomenology -the large amplification of the DM annihilation rate -requires only a single at-threshold bound state, and so imposes only a mild condition on the coupling. For instance for m n = m Z , m s = 500 GeV, we require α ≥ 0.15 which is well within the perturbative regime α < ∼ 2π. Note that, upon writing the complex scalar s in terms of its CPeven and odd parts s = φ s + ia s , we have, due to Bose symmetry, that the φ s φ s and a s a s bound states can only possess even orbital angular momentum = 0, 2, ..., while the angular momentum of the φ s a s bound states is unrestricted.
Turning to the scattering of two slow moving s particles by the Yukawa potential arising from the n scalar exchange interaction, both elastic and inelastic scattering cross sections (such as radiative capture) are amplified by a non-perturbative Sommerfeld-like enhancement. This enhancement can be formulated in terms of a nonrelativistic quantum two-body problem with a potential acting between the incoming particles. To a good approximation this leads to a dressing of the dominant S-wave part of the tree-level cross sections by a multiplicative factor, An exact analytic calculation of R for a Yukawa potential is not possible (although we give a close analytic approximation later) and so we must proceed numerically. The Schrödinger equation for the radial part of the two dark matter particle state, ψ(r), with l = 0, reads −ψ (r)/m s + V (r)ψ(r) = Eψ(r), where E = m s β 2 is the kinetic energy of the two dark matter particles in the center-of-mass frame, where each dark matter particle has velocity β. Using the outgoing boundary conditions, Rewriting r = y/m n and letting = m n /m s we can rewrite the Schrödinger equation as a function of the two ratios α/ and α/β, viz The resulting numerical solutions for R are functions of t ≡ α/ and u ≡ α/β with the 2d contour plot shown in ∼ ) there is a relatively flat region -the Coulomb region. This is the part of parameter space relevant for freeze-out, the magnitude of the enhancements being at most a factor of 3 to 5 [13,14]. More interesting is the low-velocity region in which we see the effect of resonance peaks. These peaks correspond to the formation of l = 0 bound states at threshold (E = 0), and are the focus of this analysis.
Provided we are sufficiently close to a resonance peak, the dependence of R on t = α/ and u = α/β is described by a modified Breit-Wigner resonance formula applicable for threshold resonances due to Bethe and Placzek (BP) [22,23]. Considering only elastic scattering, the Breit-Wigner resonance cross section is where Γ e and Γ are the elastic dissociation and total width for the resonant bound state and ε 0 is the distance of the resonance from exact zero-energy and is independent of β. Following BP, for small β near a threshold resonance the BW expression is modified by replacing Γ e = √ Eγ e with γ e not depending on E . Using (as R has been found by including only elastic scattering we have set Γ = Γ e ).
Eq. (5) shows that for ε 0 = 0 the cross-section first increases as 1/β 2 and then becomes independent of β for β 1 as shown in Fig 3. The plateau begins at β ∼ |ε 0 |/m s , and corresponds to σ BP e 4π/(m s |ε 0 |). Comparing the numerical calculation of the elastic cross section to σ BP e in the low-velocity limit we can extract the numerical values of γ e and ε 0 as a function of   the parameters α and for each of the possible nearthreshold bound state resonances. The plateau arises as β → 0 (for ε 0 = 0) with asymptotic value of the cross section Since The function R(t) has the form depend upon the principal and orbital angular momentum quantum numbers (nl) of the resonance that is close to threshold as t is varied. Table  I gives numerical fits for the S-wave resonances.
We remark in passing that an exact analytic treatment is possible for the Hulthén potential V H (r) = −αm n e −rmn /(1 − e −rmn ) which has similar r → 0 and r → ∞ behaviour to the Yukawa potential. The S-wave phase shifts are (see also [24]) where k = m s β is the momentum of the scattered state and A can be thought of as t = α/ up to a factor of two. From this the S-wave elastic cross section σ e = 4π |f | 2 follows using f = e 2iδ0 − 1 /2ik. Taking the β → 0 limit of σ e , the dominant behaviour is . Exactly on one of the resonances, A = A 0 , the dependence on β is σ e ∝ 1/β 2 which agrees with the BP form.
So far we have only discussed the case of elastic scattering. However, as we are interested in the indirect signals coming from dark matter annihilations we need to examine the case of inelastic scattering. In order to write down the inelastic cross section in the presence of the long range enhancements we again follow BP by writing where Γ i is the inelastic width associated with the direct annihilation or radiative capture of the incoming s-pair. This form is again only applicable when we are suitably close to one of the S-wave resonances. Following BP the inelastic width is a constant as opposed to the energy dependent elastic width. Substituting the form for the elastic width Γ e = √ Eγ e into Eq. (9) we have ] .
The first point to note here is the 1/β enhancement of the inelastic cross section compared to the elasticthe usual Bethe 1/v law [25]. This follows if we take the β → 0 limit of Eq.(10) (assuming that we are not exactly on resonance) compared to the β → 0 elastic cross section, σ BP e = πγ 2 e /m s (ε 2 0 + Γ 2 i /4). Second, for m s β 2 > ε 0 , Γ i , the inelastic cross section Eq.(10) shows a 1/β 3 dependence which plateaus to the standard 1/β dependence when m s β 2 ε 0 , Γ i . This behaviour is exactly that of Eq.(2) for the enhancement of the naïve inelastic cross section given the β dependence of R plotted in Fig. 3 resulting from numerical solutions and as discussed above. Third, the final limiting value of βσ BP i has a 1/ε 2 0 dependence and shows that for WIMPonium bound states close to threshold the cross sections are further enhanced. (In the above we have made the assumption that ε 0 > Γ i which is true over the majority of parameter space. However, for the case where ε 0 < Γ i the true dependences can be more complicated involving cross terms of the inelastic and elastic widths.) Depending on the size of ε 0 the size of the 1/β 2 enhancements can be significant since for DM annihilations in present day astrophysical systems the relevant range of β is ∼ 10 −3 − 10 −5 .
A useful way of thinking about and calculating approximately the amplification of the elastic and inelastic cross sections is in terms of the diverging scattering lengths that occur when a bound state energy tends to zero. Recall from the elementary theory of non-relativistic elastic scattering that the S-wave phase shift satisfies where r e ∼ 1/m n is the effective range of the potential and l s is the scattering length which is related to the near-threshold bound state energy ε 0 by These equations imply cot δ 0 = − |ε 0 |/E where E = k 2 /m s is the CM scattering energy, and we have assumed k → 0 before ε 0 → 0. From the standard expression of the elastic cross-section in terms of δ 0 one then finds agreement with the BP formula Eq.(5) in the same limit of k → 0 before ε 0 → 0. Comparing with the BP form we learn that γ 2 e = 4|ε 0 | → 0 as the bound state approaches exactly the zero-energy threshold. The advantage of this approach is that the S-wave component of the distorted incoming plane waves can be simply expressed in terms of δ 0 and thus ε 0 . To a good approximation the wavefunction is This form and its dependence on the scattering length s will be useful to us when we discuss radiative capture.

III. DECAYS, CAPTURES, & TRANSITIONS
WIMPonium bound states possess a rich phenomenology of radiative captures to various bound states, and decays and transitions from or among the bound states.

A. Radiative Capture
For interesting bound-state to bound-state transitions to be relevant, the radiative capture cross section must be significant. Similar to elastic scattering the radiative capture cross section is enhanced when there is a nearthreshold bound state, as is approximately the situation in neutron-proton scattering where enhanced radiative capture to a bound deuteron is possible.
We will consider transitions into both the = 0 and = 1 bound state energy levels. The most economical way to perform such a calculation is to first expand the continuum state in terms of partial waves where we will keep only the S-partial wave, Eq. (14) (the incoming P-wave gives terms that are suppressed by β 2 ). Labeling the bound state wavefunctions into which we will be capturing as ψ nllz , the two most important states are ψ 100 ≈ (πa 3 0 ) −1/2 exp(−r/a 0 ) and ψ 210 ≈ (32πa 3 0 ) −1/2 r cos θ exp(−r/2a 0 )/a 0 , where we've assumed that the Ω (n, ) bound state wavefunctions are similar to hydrogen with a 0 ∼ 2/m s α. (The bound states ψ 200 and ψ 21±1 do contribute, with the total rates for radiative capture changing by, at most, an O(1) factor. As we are interested in the general parametric dependence we neglect the contributions from capture into these bound states.) The radiative capture cross section depends on matrix elements of the form I = ψ f |O|ψ i , where ψ i (r) is the initial distorted partial wave of the continuum state, ψ f (r) is the wavefunction of the bound state into which we are being captured, and O is the interaction Hamiltonian. In our case O = λ 2 e ip.r , where p is the momentum of the radiated scalar n state.
Expanding the exponential in powers of p.r and considering radiative capture into the 1s and 2p states, the overlap integrals take the forms where the first non-zero integral comes at second order in p.r for capture into the 1s state and at linear order for capture into the 2p state. From this the rates are found to be where p 1 ≈ E 2 1 − m 2 n and p 2 ≈ E 2 2 − m 2 n are the momenta of the n scalars due to transitions into the 1s and 2p states respectively. For simplicity, from now on we will assume that m n is small relative to transition energies and so can be neglected. Clearly if m n is not small then there are trivial kinematical suppression factors.
Using the S-wave phase shift, Eq. (12), and taking the small k limit gives δ 0 ≈ k s . For large scattering length (compared with 8a 0 ) the rates become where in the last expression we have used the relation, Eq. (13), between the scattering length s and the near threshold bound state energy and continuing the analogy with hydrogen we have taken the bound state energies to be E n = −m s α 2 /4n 2 . This implies radiative capture cross sections This shows the usual 1/β Bethe-dependence of an inelastic cross section near β → 0, and most importantly an additional α 2 m s /ε 0 enhancement relative to the radiative capture cross section if there was not a near threshold bound state (in other words if the a 0 -dependent terms dominated over δ 0 /k in Eq. (16)). The factor α 2 m s /ε 0 is just the ratio of the typical bound state energy compared to the near-threshold energy. In addition to this, it should be noted that for α ∼ O(1) (which is still well within the perturbative regime) the radiative capture cross sections are further increased, although the simple hydrogen-like scaling that we have employed starts to break down and numerical methods must be used.

B. Decays
Consider annihilation of the ss bound state Ω (n, ) to light (mass m s ) degrees of freedom. Let the amplitude for the free 2 → 2 scattering be A f i , then from the standard theory of, eg., positronium decay the amplitude M f i for the bound state decay is where the momentum-space bound state Schrödinger wavefunction is normalized as d 3 p|ψ(p)| 2 = 1. Here p is defined by (0, 2p) = q 1 − q 2 where q i are the 4-momenta of the two s constituents. In the limit of relatively weak f i +... and using ψ(p) = d 3 xe ip.x ψ(x)/(2π) 3/2 gives which enables the calculation of the decay width of various orbital angular momentum bound states. For m s /m n 1, and the ssn coupling λ not very large, the wavefunctions of the bound states are approximately hydrogen-like, which upon applying Eq. (21) gives the annihilation width = 0 This expression is most accurate for the more tightly bound WIMPonium states with small principal quantum number, n = 1, 2, ..., while for higher bound states a precise decay width requires numerical evaluation of the bound state wavefunction. Parametrically Eq. (22) gives a good estimate of the decay width in all cases. From Eq. (21) and the expansion of the free 2 → 2 scattering amplitude, the decay widths of the higher orbital angular momentum states are parametrically suppressed by powers of 1/(m s a 0 ) 2 . Using hydrogen-like scaling once again shows that, eg, Γ =1 ann /Γ =0 ann α 2 which, if we take α < 1, is parametrically small. However, if we take α ∼ O(1) all decay widths can be large and comparable.

C. Transitions
Transitions between the various bound states are possible with either the emission of the CP = ±1 components of n, or if on-shell n-production is kinematically disallowed, by decay to light SM states through virtual n-emission. First assuming on-shell n emission (possible when α 2 m s > ∼ m n ), the relevant WIMPonium matrix element reduces to T f i = R n Y m | λ exp(ipr) |R n Y m . Similar to transitions in hydrogen-like systems we may expand the exponential in powers kr ∼ α (for relativistic n). The first transitions occur at O(kr) with the emission of the CP = −1 state of n with orbital angular momentum l = 1, changing the bound state from φ s φ s or a s a s to φ s a s or vice versa.
A good estimate of the various transition rates follows from a straightforward application of Fermi's Golden Rule. Parametrically the rate for ∆ = 1 transitions emitting a relativistic a n scales as where ∆E ( m n ) is the energy splitting between the bound states. We see the (∆E) 3 behaviour familiar from hydrogen-like systems which favours deep transitions. All other transitions emitting an on-shell n are down by powers of (∆Ea 0 ) 2 ∼ α 2 . For instance ∆ = 0 or 2 transitions emitting a CP = +1 φ n state are suppressed as Γ ∆ =0,2 trans /Γ ∆ =1 trans ∼ α 2 as well as by final state phase space factors (∆E/∆E ) 3 . Transitions to light SM states via virtual n's are further suppressed by both couplings such as α or α em and (at least) three-body phase space factors.
The most important feature of Eq. (23) is that it shows that the transition rate between different bound states can be competitive to the direct decay rate, Eq. (22), of the = 0 tower as long as on-shell n production is kinematically possible. If the s DM particles are captured in a P-wave orbital angular momentum state, then the extra suppression of P-wave annihilations implies that transitions via on-shall n's can dominate. Furthermore transitions to light SM states via virtual n's might be only mildly suppressed relative to P-wave annihilations depending on the strength of the coupling of n to Higgs (and thus other SM) fields. Fig.4 schematically illustrates the dominant decays and transitions and their respective rates while Fig.5 depicts diagrams responsible for discrete γ lines.

IV. CONSEQUENCES FOR INDIRECT DETECTION AND LHC
A wide range of consequences follow from the existence of WIMPonium bound states. As we have already mentioned they provide a new, dynamical mechanism for the "boost factors" that are often introduced to explain anomalies in indirect detection observations. Unlike traditional ρ(x) 2 enhancements they depend on the velocity distribution of the DM particles, and since between the freeze-out epoch and today the velocity changes from β f o ∼ 1/5 to β now ∼ 10 −3 − 10 −5 , the near-threshold bound state enhancement decouples the value of the annihilation cross section determined by successful thermal freeze-out from that observed today in indirect annihilation observations. (See [14] for a discussion of thermal freeze-out production of DM in models closely related to our toy theory.) In fact, because of thermal corrections to the couplings and masses during freeze-out the "path" taken in (α/β, α/ ) parameter space as the universe cools is not exactly a α/ =constant line, but instead can move towards or away from the value at which a resonance occurs exactly at threshold. This further decouples the value of the cross-section at freeze-out from that observed now. Moreover, the enhancement of the cross-section as the universe cools leads to the possibility that there is interesting residual post-freeze-out annihilation of the DM particles, for instance leading to changes in BBN predictions of such elements as 6 Li and 7 Li [26]. Depending on the maximum size of the threshold enhancement, which is determined both by the size of the inelastic widths, and the degree to which the bound state approaches zero energy, this dynamical boost factor can be large enough to potentially favour environments such as dwarf galaxies, since their velocity dispersions go down as low as 3-5kms −1 [27] compared to the typical galactic value of ∼ 200 − 300kms −1 .
Turning to particle physics model-building issues, the existence of WIMPonium bound states implies that the standard supersymmetric neutralino DM picture must be modified somewhat. Although we have explained the phenomenology of WIMPonium in the context of a very simple, purely scalar toy model, we emphasize that similar phenomena are possible if the DM is fermionic, or even neutralino DM. From the condition for the existence of at least one bound state α ≥ 0.84m n /m s we learn that if the DM is a neutralino interacting via the exchange of W ± and Z gauge bosons (and Higgs states) then the neutralino must be heavy m neutralino > ∼ m W /α 2 ∼ 2 TeV. On the other hand if the DM particle interacts with a new strong-interacting sector, say a hidden valley sector [28], then the DM particle can potentially have close-to-weakscale mass.
A particularly attractive possibility is to have the DM state associated with electro-weak symmetry breaking dynamics in some way, so that it interacts with the Standard Model via the so-called Higgs Portal [14], and has strong Higgs-mediated self interactions. This is in the same class of models as our toy theory, although the DM particle s could be fermionic in which case the spectrum of bound states and associated decays and transitions is even richer. In all three cases the LHC search strategy for the DM state is greatly modified compared to the standard expectation of weak-scale neutralino DM.

V. CONCLUSIONS
We have shown that theories of TeV-scale physics can have dark matter candidates whose annihilation proceeds via the formation of near-threshold WIMPonium bound states. Depending upon the closeness-to-threshold of the weakest-bound state, these can lead to a substantial velocity-dependent amplification of the dark matter annihilation cross-section, preferring the lowest velocity dispersion environments all other factors being equal, and providing a new dynamical source of "boost factors" for indirect detection signals. In addition, the amplified radiative capture to more deeply bound states, and the transitions among such bound states can both lead to a rich spectrum of discrete γ lines which, if observed, would give striking confirmation of the mechanism.
During the preparation of this work, [11] appeared. This paper also considers aspects of the phenomenology of WIMP bound states in the context of their annihilations and the consequences for dark matter indirect detection.