Self-consistent Calculation of the Sommerfeld Enhancement

A calculation of the Sommerfeld enhancement is presented and applied to the problem of s-wave non-relativistic dark matter annihilation. The difference from previous computations in the literature is that the effect of the underlying short-range scattering process is consistently included together with the long-range force in the effective QM Schr\"odinger problem. Our procedure satisfies partial-wave unitarity where previous calculations fail. We provide analytic results for some potentials of phenomenological relevance.


Introduction
The annihilation cross section of thermal relic weakly-interacting massive particle (WIMP) dark matter controls both the relic abundance through freeze-out in the early Universe as well as the prospects for indirect detection. It is important to establish to what extent the annihilation cross section at freeze-out, when WIMPs are mildly non-relativistic with thermal velocities of order v ∼ 0.2, is related to the annihilation cross section in systems where the WIMP is highly non-relativistic, including the Galaxy (v ∼ 5 × 10 −4 ) or the epoch of recombination (v ∼ 10 −10 for a WIMP with ∼ TeV mass and kinetic decoupling temperature of a few MeV).
A well-known effect that can modify the annihilation cross section in the deep non-relativistic regime is the Sommerfeld effect (SE), arising when the WIMP couples to a light force mediator [1][2][3][4]. The SE was estimated in the literature using a factorized formula (see, e.g. Ref. [5,6]), σv = σv 0 × S(v). (1) Eq. (1) is computed in two steps: 1. σv 0 is the velocity-weighted hard annihilation cross section, computed using the short-range physics alone, ignoring the long-range potential 2. The SE factor S(v) ≡ |ψ 0 (0)| 2 is derived from ψ 0 (x), the wave function for the nonrelativistic two-body dark matter system, determined as a solution to the quantum mechanics (QM) scattering problem using a Schrödinger equation in which the short-range physics is omitted, with asymptotic behaviour ψ 0 (r → ∞) → e ikz + f e ikr r . Eq. (1) was used in many examples in the literature, but in certain situations it breaks down. A noted example is the occurrence of a bound state with zero binding energy in the long-range Schrödinger equation. (See, e.g., §133 in Ref. [7].) If a bound state solution exists with binding energy |E b | much smaller than the typical potential energy V 0 , then the SE factor approaches the form 1 S(v) ∼ V 0 µv 2 /2+|E b | . Since σv 0 approaches a constant for s-wave hardannihilation process (the Bethe 1/v law), Eq. (1) gives velocity-weighted cross section σv ∝ 1/v 2 at small v for |E b | → 0. This violates parametrically the partial-wave unitarity limit [7,10], σ ≤ π/k 2 = π/(µv) 2 . A more obvious example in which Eq. (1) fails, is when σv 0 is itself not very small, such that a finite S(v) can lead to σv violating the unitarity bound at finite velocity without any special resonance of parametric origin.
In the current paper we identify the source of trouble in Eq. (1) and fix it in a simple way. The basic result is easy to explain: probability conservation requires that the naive SE factor S(v) = |ψ 0 (0)| 2 , given by the modulus of the wave function at the origin neglecting the shortrange potential, should be replaced by |ψ(0)| 2 , the modulus of the full wave function including the short-range interaction 2 . Our main task in the paper will be to compute ψ(0) and to show how using it instead of ψ 0 (0) in Eq. (1) restores unitarity.
Our formula for the full annihilation cross section, given in Eq. (40), is only slightly more complicated than Eq. (1). To compute the full SE we require one more quantity [denoted T (v) in Eq. (40)] based on the long-range interaction, in addition to the usual S(v); and one more quantity [the short-range elastic scattering cross section, denoted σ sc,0 in Eq. (40)] based on the short-range interaction, in addition to the hard annihilation cross section σv 0 . We provide analytic expressions of the new T (v) term for specific long-range potentials of interest in typical applications, including the Coulomb and well potentials and the Hulthén potential (approximating the Yukawa potential), for which the usual S(v) term was derived elsewhere in the literature.
In Sec. 2 we derive our main result and present a general formula for the full SE. In Sec. 2.2 we show that close to a resonant point, for perturbative underlying short-range cross section, the full SE can be obtained by simple velocity shift S(v) → S(v + v c ) in Eq. (1). The shift velocity v c contains a product of short-and long-range information that we compute explicitly for the Hulthén potential. In Sec. 3 we provide analytic results for specific potentials. We focus in particular on the Hulthén potential, approximating the Yukawa interaction, where we demonstrate how our SE computation conserves partial-wave unitarity. In Sec. 4 we conclude. 1 Ref. [8] discussed this limit using a modified Breit-Wigner formula introduced by Bethe and Placzek [9]. See also §145 in Ref. [7]. 2 A similar situation was discussed in the context of nuclear physics, e.g., non-relativistic proton-proton fusion [11,12]. There, the emphasis was on resummation of strong short-range elastic scattering. Here we address the problem from a somewhat different perspective, and will mainly focus on the resummation of short-range inelastic scattering (annihilation).

Formulation
The non-relativistic effective action for the dark matter two-body system is obtained by integrating out degrees of freedom with large momentum [2][3][4]. The resulting Schrödinger equation gives an effective QM description for momentum p ≪ µ, Here µ = M/2 is the reduced mass. Long-range light mediator exchange gives rise to the real, rotationally-symmetric potential V (r), r = |x|. The short-range process is encoded in the term uδ 3 (x) with complex parameter u. We note that Eq. (2) is adequate for s-wave shortrange scattering, which is what we address here and is enough to establish the formalism. Generalization to higher-l initial state is left for subsequent work.
Our objective is to solve Eq. (2); the annihilation cross section can then be read from the solution as follows. The probability current associated with ψ(x) is j(x) = 1 µ Imψ * (x)∇ψ(x). Using the Hamiltonian H = − 1 2µ ∇ 2 +V (|x|)+uδ (3) (x), and imposing steady state ∂ t |ψ(x)| 2 = 0, we have The full velocity-weighted annihilation cross section is given by integrating the divergence of the probability current, The boundary condition on ψ(x) is that at large r it should be given by an incident plane wave plus outgoing spherical wave, ψ(x) → e ipz + f e ipr /r. The factor v = p/µ on the LHS of Eq. (4) is the relative velocity between the incoming particles at z → ∞, namely, the flux of the incident wave. We are already in place to compare with earlier literature.
• The full solution ψ(x), where we solve Eq. (2) retaining both the short-and the long-range potentials, gives the full SE as |ψ(0)| 2 . The shortcoming of Eq. (1) is that it fails to conserve probability by using ψ 0 instead of ψ.
We solve Eq. (2) following Jackiw [13] that analyzed the scattering amplitude for the QM delta-function potential. We define the Green's function G(x) and the wave function ψ 0 (x), chosen to satisfy The function ψ 0 (x) is the regular solution to the Schrödinger equation that corresponds to a unit-normalized initial plane wave; that is, ψ 0 (0) is finite, and ψ 0 (x) → e ipz + f 0 e ipr /r at large r, where the second term represents a contribution to the scattered wave. The function G(x) is chosen to give an out-going spherical wave at r → ∞. Namely, the boundary condition for G at infinity is written as where d p is a function of p.
The solution of Eq. (2) is so consistency at the origin implies leading to Using Eq. (9) in Eq. (4) gives us the full SE directly: The denominator corrects Eq. (1) and shows in what way factorization needs to be modified. However, we are not quite there yet, because G(0) is divergent. This divergence is regulated with a redefinition of the parameter u absorbing an infinite constant, as we now discuss.
Rewrite G(x) as From Eqs. (5) and (7), g p (r) can be shown to satisfy 3 We are interested in G(x) near the origin. We limit the discussion to a long-range potential V (r) that does not diverge at the origin faster than 1/r, and can be expanded as 3 To derive Eq. (13), we used Eq. (12), limr→0 rg ′ p (r) = 0 and This type of potential includes the Yukawa, Coulomb, and well potentials. Near the origin g p (r) and G(r) are found as Note that the divergent part of G(0) is independent of momentum p. In addition, the divergence is restricted to ReG(0), while ImG(0) is finite for real V (r).
To subtract the divergence we define a renormalized coupling at reference momentum p 0 , The subscript on G p 0 is there to emphasize that this G is to be evaluated at p 0 . In QFT language, (2µu) −1 is a bare coupling and ReG p 0 (0) is a counter term. Deferring some subtleties to the next subsection, it is convenient to take large reference momentum p 0 and to define the short-range cross section σv 0 as the cross section measured at p 0 . We are left with the task of computing k p 0 from this matching procedure; this we do in the next subsection by appealing to some more basic QM. Before moving on, note that above we used a regularization scheme in which G(0) is replaced by G(ǫ), finally taking ǫ → 0. Jackiw [13] took a different approach, examining both momentum cut-off and dimensional regularization. But then, Jackiw [13] was interested in the scattering amplitude for the delta function potential alone, while we are interested in the interplay between short-and long-range physics. To this end we find our procedure convenient.

Cross section matching
Rewrite Eq. (10) as We take the partial wave expansion as (See, e.g., §123 in Ref. [7]) where the P ℓ are Legendre polynomials and cos θ = z/r. Considering s-wave scattering we define the corresponding radial wave functions Since ψ 0 (0) is finite, the boundary condition for χ s,0 is χ s,0 (0) = 0. This gives The function χ s,0 (r) solves the same 2nd order differential equation as the function Img p (r), with the same (null) boundary condition at the origin. The other boundary condition at r → ∞ only fixes the normalization. Therefore χ s,0 (r) = cImg p (r), where c is an (in general complex) constant multiplicative factor. Thus, for large r, Using Eq. (24), the s-wave scattering amplitude is The elastic scattering, annihilation, and total cross sections are (see, e.g., §142 in Ref. [7]) In terms of g ′ p (0) and k p 0 , σ ann is written as The elastic scattering cross section is To isolate the short-range physics we need to consider the limit of large kinetic energy p 2 /2µ compared to the typical potential energy. To be precise, we define V (r) = V −1 /(2µr) +Ṽ (r), whereṼ (r) is assumed to be regular everywhere. The large momentum limit we are interested in satisfies p 2 ≫ 2µṼ (r) for all r. In this limit the solution for g p (r) that satisfies g p (0) = 1 at the origin, and that matches to an outgoing wave in the leading WKB approximation at large r, is found as where W (k, m, z) is the Whittaker W function. The dots (...) stand for sub-leading corrections to the WKB approximation. For p, p 0 ≫ V −1 , Thus σ ann ≃ 4π p σ sc ≃ 4π Recall that we have defined the renormalized inverse-scattering length parameter k p 0 in Eq. (18) by appealing to the Green's function at reference momentum p 0 . We see that if we set p = p 0 and consider large p 0 w.r.t. the long-range potential, in the sense discussed above, we can use Eqs. (31-32) to fix k p 0 , up to the sign of Rek p 0 , from the cross section defined at p 0 .
In perturbative examples life is easier because 4 Imk p 0 is larger than µ (and so larger than p 0 ) and V −1 is smaller than µ. Assuming Imk p 0 ≫ p 0 , and The sign of Rek −1 p 0 is model-dependent, in the sense that positive (negative) sign corresponds to a repulsive (attractive) short-range force.
In what follows we focus, for simplicity, on the perturbative regime where Eqs. (33-34) are valid and our results take a particularly simple form.
We still need to compute the long-range effect encoded in g ′ p (0) and d p . We first note that where S(v) is the usual naive SE factor from Eq. (1). To prove this, define the Wronskian satisfying W ′ p (r) = 0. At the origin and at r → ∞ we have Together with g p (0) = 1, derived before, this gives Eq. (35). It is useful at this point to refer again to the usual formula Eq. (1). This should arise as the limit of our Eq. (27) for |k p 0 | ≫ |g ′ p (0) − g ′ p 0 (0)|. Using Eq. (35) we have, in this limit, as needed.
Last we define In general, for arbitrary small momentum p, Reg ′ p at the origin must be determined by solving Eq. (12) using the full long-range potential; this situation is similar to that of the usual S(v) factor. Below we will give analytic solutions for specific examples: (i) the Hulthén potential, approximating the Yukawa potential; (ii) the Coulomb potential arising from a massless mediator; (iii) the well potential.
Putting everything together, we can write Eq. (27) as where η = ±1 is the sign of Rek −1 p 0 . An analogous expression can be written for the elastic scattering cross section, but here and in what follows we focus on annihilation.
This concludes our basic calculation: the SE factor for annihilation, including regularization by the short-range process, is given as the coefficient of σv 0 in Eq. (40). In the next section we provide analytic computations of the SE for specific examples. Analytic calculations and approximations for the S(v) term exist in several references; we add to those a calculation of the T (v) term. The short-range information: σv 0 , σ sc,0 , and the sign η of the short-range force, should be provided according to the underlying model. For example, a thermal-freezout WIMP model would have σv 0 ≈ 3 × 10 −26 cm 3 /sec, etc.

Regularising velocity near a zero energy bound state
Particularly strong SE arise when the Schrödinger equation including the long-range potential admits a bound state solution with zero binding energy |E b | → 0 [2][3][4]14]. At small velocity near |E b | → 0, the naive SE factor scales as where V 0 is the characteristic scale of the potential energy. Ref. [2][3][4] noted that the singular 1/v 2 scaling is cut-off when the short-range contact term is accounted for perturbatively in solving the Schrödinger equation, and that the SE saturates at some velocity v c . Using Eq. (40) we can verify these observations directly. For simplicity and for the sake of comparison with Refs. [2][3][4] we neglect the real part of the inverse-scattering length induced by the short-range interaction. This amounts to setting σ sc,0 = µ 2 σv 0 4π σv 0 in Eq. (40). Plugging the asymptotic form Eq. (41) into Eq. (40) we find Thus we derive and confirm the use of a regularising velocity 5 . This can be plugged in the usual formula Eq. (1) using The velocity scale v c is useful for finding the SE saturation point. To assess this scale for phenomenological applications we examine it for the Yukawa potential V Y (r) = −αe −mr /r arising from a light mediator with mass m ≪ M , where M = 2µ, again, is the WIMP mass. In We scaled the high-velocity annihilation cross section to the relevant value for thermal freeze-out dark matter. For parameters in the ballpark chosen in Eq. (44) with α = O(1), the numerical value of v c makes it phenomenologically relevant for indirect searches for annihilating dark matter, implying that the SE would be saturated in Galactic systems. For comparison, the typical virial velocity in the Milky Way halo is about a factor of 3 smaller than this v c , while in dwarf galaxies the WIMP velocity would be much smaller. Models where these parametrics apply include, for example, the dilaton portal of [15].
In the early universe, before the dark matter forms virialised halos, the WIMP velocity may be much lower than in a typical Galaxy and the regularization may be important even for TeVscale or lighter WIMPs with perturbative interactions. However, away from resonances, the saturation velocity of the usual S(v) factor (v ∼ m/αM ) tends to greatly exceed the regularization velocity v c , for α ≪ 1. Accordingly, for perturbative interactions, pronounced effects of the regularization in the early universe will be confined to the small regions of parameter space near resonance peaks.
It has been previously noted that in such resonant regions, annihilation may recouple during the radiation-dominated epoch after the dark matter kinetically decouples from the Standard Model, since the ∼ 1/v 2 scaling of S(v) leads to an increasing annihilation rate per comoving volume [16]. This effect can lead to a second phase of exponential depletion of the dark matter density, meaning that the correct relic density can be obtained with a much smaller short-range annihilation cross section. In this case, it can be very important to take the regularization into account; examples that we have calculated indicate that the short-range annihilation cross section required to yield the correct relic density in such a resonant scenario can increase by several orders of magnitude, relative to the case with no regularization, albeit in all cases it remains well below the "standard" thermal relic value of σv ≈ 3 × 10 −26 cm 3 /s.

Hulthén potential
The Hulthén potential is We consider this potential because, choosing m * = (π 2 /6)m, it provides a good approximation of the commonly encountered Yukawa potential V Y (r) = −αe −mr /r with consistent asymptotic behaviour at large and small r.
The SE with Eq. (45) admits an analytic solution (see Refs. [17,18] for a computation of the S(v) term). Here we add the ingredients needed for the full formula Eq.(40). Define dimensionless variables, Using these, g p (r) is written as the Gauss hypergeometric function. This gives where H(z) is the Harmonic Number. It is easy to verify that |d p | 2 = 1 p Img ′ p (0) = S(v). In Fig. 1 we show the SE for the Hulthén potential. On the left we show the naive SE given by S(v) of Eq. (50), corresponding to the usual computation in the literature. The regulated SE computed from Eq.(40) is superimposed, but the difference from the usual method is invisible on the scale of the plot. On the right we zoom-in on one of the resonant peaks. The regulated SE of Eq.(40) is the smooth lower curve. For this plot, we took α = 1, v = 10 −6 , and shortrange cross sections σv 0 = 1/(32πM 2 ) and Rek p 0 = 0 corresponding to short-range elastic cross section σ sc,0 = (µ 2 σv 0 /4π)σv 0 .

Coulomb potential
The SE for the attractive Coulomb potential, V (r) = −α/r, can be obtained from the result of Sec. 3.1 by taking the mass parameter m * → 0. This leads to where ψ(z) = Γ ′ (z)/Γ(z) is the PolyGamma function. We do not quote here the result for the factor d p ; the normalization of d p is given of course by |d p | 2 = S(v), but obtaining the phase requires a modification to the plane wave asymptotic states that we have assumed in Eq. (14). This is the standard logarithmic asymptotic phase modification due to the slow-falling 1/r potential (see e.g. §7.13 in Ref. [19], and §135 in Ref. [7]). Note that d p is not needed in the full SE formula for annihilation, Eq. (40), for which one only requires S(v) and T (v).

Well potential
Consider the attractive well potential, g p (r) is written as wherep = p 2 + p 2 V . c p and d p are determined from continuity of g p and g ′ p at r = R. g ′ p (0) is obtained as, Then,

Conclusion
Several references in the literature noted that the naive factorized estimate of the Sommerfeld effect (SE), given by σv ≈ σv 0 × S(v) [see Eq.
(1) and discussion around it], fails in specific parametric points where the Schrödinger equation for the long-range force develops a bound state with zero binding energy. At these isolated points, the naive calculation violates partialwave unitarity. We add to this the observation that a large but finite S(v) factor can boost σv, naively estimated as above, such that it violates the unitarity limit even if the short-range process itself is perturbative, σv 0 ≪ π/µ 2 v 0 . Also in the literature, referring specifically to the zero energy bound state problem, it was commented that the apparent unitarity violation should be regularised by the finite width of the zero energy quasi-bound state. In the analyses of Hisano et al [2][3][4] a full Schrödinger equation is derived from the quantum field theory (QFT), containing explicitly the non-hermitian contact term representing short-range annihilation; it is then argued there that the apparent unitarity violation is an artefact of treating the contact term perturbatively in the solution.
We identify the problem from another angle, taking the effective quantum mechanics (QM) perspective where it is clear that the usual formulation of the SE does not conserve probability, because it does not account for the probability flux absorbed at the origin r → 0 by the shortrange annihilation process.
We solve the problem by referring to Jackiw's [13] calculation of the scattering amplitude for a delta function potential. In this respect, our contribution is in that we fulfil the suggestion of Hisano et al [2][3][4] to solve completely the contact term effects. The final expression for the full SE factor for annihilation is given in Eq. (40).
Calculating the full SE with our formalism is as straightforward as it is in the usual naive formula. A technical difference is that the full result requires one more quantity encoding the long-range force [contained in the new term T (v), a natural counterpart to the usual S(v) term]; and one more quantity encoding short-range physics (the short-range elastic scattering cross section σ sc,0 ). In Sec. 3 we gave analytic expressions of the new T (v) term for the Hulthén potential (a close approximation to the Yukawa potential) as well as for the Coulomb and well potentials. We also recap the formulae for S(v), found in existing literature for these potentials.
For perturbative underlying short-range cross section σv 0 ≪ π/µ 2 v 0 , under the effect of long-range Yukawa, Hulthén, or well potentials, a particularly simple prescription is given in Sec. 2.2: the full SE is obtained by the velocity shift S(v) → S(v + v c ) with v c given in Eq. (43).
While we have clarified and solved the basic problem with the SE at large amplification, several issues were left out in the current paper and merit further work. These include (i) generalization to the multi-state problem, relevant for e.g. the electroweak super-partner scenario of [2][3][4]; (ii) generalization to p-wave and higher-l initial states; (iii) complete characterization of the SE for elastic scattering.