Quantum theory of dark matter scattering

Dark matter self-scattering is one of key ingredients for small-scale structure of the Universe, while dark matter annihilation is important for the indirect measurements. There is a strong correlation between the velocity-dependent self-scattering cross section and the Sommerfeld enhancement factor for the dark matter annihilation cross section. In this study, we formulate a direct relation between them by the use of Watson’s (initial state/final state) theorem and Omnès solution, and our formulation reproduces the Sommerfeld enhancement factor, which directly computed by solving the Schrödinger equation, from the scattering phase shift.


Introduction
Interaction is one of the most important properties of particles.The interaction of dark matter (DM) is still unknown even though the existence of DM has been firmly confirmed by astrophysical observations.In this study, we focus on two kinds of interactions of DM: annihilation and self-scattering processes.
Recently, there has been growing interests in DM self-interaction.We can prove the selfinteraction of DM in cosmological observations though it is extremely challenging to find the self-interaction of DM in terrestrial experiments.DM self-scattering may explain the smallscale structure problems of the Universe [1] (see Ref. [2] a review).It has been extensively studied that these problems may be resolved even within the collisionless cold DM paradigm once we take into account baryon feedback processes in simulations.However, it is still under debate to what extent the baryon feedback processes accommodate the problems, and hence it is intriguing possibility that DM has self-interaction.Combining observations of DM halos in a wide range of masses may indicate the velocity dependence of the DM self-scattering cross section [3].The self-scattering cross section per mass is preferred to be σ/m ∼ 0.1 cm 2 /g at the galaxy clusters, where the collision velocity is of O( 103 ) km/s, in order to explain the core formation of the galaxy clusters [4] and to avoid the constraints of the merging clusters [5,6].As for the galaxy scale, where the velocity is of O( 102 ) km/s, σ/m ∼ 1 cm 2 /g is preferred to explain the galactic rotation curve [7,8].The DM self-scattering in the Milky Way (MW) dwarf spheroidal galaxies, where a typical velocity is v ∼ 30 km/s, is still under debate.The preferred scattering cross section per mass is widely distributed as 0.3-40 cm 2 /g [9].The expected cross section can be enlarged if DM halos are in the core-contraction phase: 30-120 cm 2 /g [10].On the other hand, the ultra-faint dwarf galaxies, which are considered to be more DM-dominated, provide constraints on the cross section to be σ/m ≲ 0.1 cm 2 /g [11].If the preferred scattering cross section per mass is of O(1-10 2 ) cm 2 /g for the typical velocity of v ∼ 30 km/s, DM models with a light mediator are one of the simplest models to explain the velocity-dependence of the self-scattering cross section.
Meanwhile, the DM annihilation into the standard-model (SM) particles (or the mediator particles) is important for determining the DM abundance today in some scenarios with the freeze-out mechanism.The annihilation rate today is also important for the indirect searches for DM.The annihilation cross sections for these processes can significantly differ due to different velocities even if the microscopic interactions controlling these processes are same.A typical velocity of DM is v/c ∼ O(0.1) at freeze-out, while the maximum velocity of DM is v/c ∼ O(10 −2 ) at the cosmic structures.When the DM particle couples to a light force mediator, the annihilation cross section can be significantly enhanced in the low velocity regime [12][13][14][15][16], which is known as the Sommerfeld enhancement [17].The indirect measurements of DM put strong constraints on the DM annihilation into the SM particles, such as e ± and γ.On the contrary, the velocity-dependent annihilation cross section has been often utilized to explain the observed cosmic-ray excesses for a decade after the positron excess (such as PAMELA [18], Fermi-LAT [19], and AMS-02 [20]) was reported [16].
In the context of the self-interacting DM (SIDM), the presence of the light mediator leads not only to the velocity-dependent self-scattering but also to a large annihilation cross section.When the Sommerfeld-enhanced DM annihilation ends up with the energy injection into the electromagnetic channels, the indirect measurements put strong constraints on the annihilation cross section [21].The Sommerfeld enhancement factor can also be sizable for the SIDM model providing the self-scattering cross section that almost saturates the Unitarity bound, as pointed out in Ref. [22].There can be a connection through quantum mechanical properties between the annihilation process and the scattering process.The Watson's (final state/initial state) theorem [23] states that the discontinuity of a weak-coupling annihilation amplitude is given by the phase shift of the relevant scattering process.The annihilation amplitude is given by socalled Omnès solution as far as the scattering process is in the elastic regime.In this study, we apply the Omnès solution to a system that contains the DM particle χ and the light mediator particle ϕ.The amplitude for DM annihilation into ϕ particles is related with the phase shift of the self-scattering process via the Watson's theorem.
The organization of this paper is as follows.In the next section, we introduce the standard methods to compute the phase shift for the DM self-scattering and the Sommerfeld enhancement factor for the DM annihilation under the attractive Yukawa potential, and we confirm that there is a strong correlation between the resonant self-scattering and the large Sommerfeld enhancement factor.In Section 3, we introduce the Watson's final-state/initial-state theorem and the Omnès solution, which relate the behavior of the scattering phase shift and the annihilation amplitude.We also provide a new computational method for the Sommerfeld enhancement factor via the Omnès solution.One of key ingredients for our formulation is the number of bound states and corresponding poles of the amplitudes.We summarize the phase shift and the bound states by introducing the effective range theory in Section 4. We can find the pole of the annihilation amplitude corresponding to the shallow bound state by the use of the effective range theory.In Section 5, we demonstrate the agreement of the Sommerfeld enhancement factors in two methods.Even though we focus on the Yukawa potential in the main text, we also discuss the attractive Hulthén potential in Appendix A. Section 6 is devoted to conclude our study.

Self-Scattering and Annihilation
In this section, we introduce the standard computational methods for the phase shift of the DM self-scattering and for the Sommerfeld enhancement factor of the DM annihilation in the presence of the long-range force.We consider a system composed of DM particle χ (with the mass of m χ ) and mediator particle ϕ (with the mass of m ϕ ) interacting through the attractive Yukawa potential1 : Here, r denotes the relative distance between two DM particles, and α is a fine-structure constant of the long-range force.We define the dimensionless parameters, a ≡ v/2α and b ≡ αm χ /m ϕ with v being the relative velocity.First, we discuss the DM self-scattering, χχ → χχ, described by the Yukawa potential.For DM self-scattering with low velocity mediated by the light particle, neither the Born approximation (b < 1, for which the cross section is computed perturbatively) nor the classical approximation (ab > 1, for which a classical cross section describes the DM scattering) is valid.The non-perturbative effect by a multiple exchange of the light mediator becomes important for the parameter space, b > 1 and a < 1.Since we do not have an approximate formula for the cross section σ for such parameter space, we must solve the Schrödinger equation for the DM two-body system, with the partial-wave analysis in order to take into account the non-perturbative effect.Here, µ = m χ /2 is the reduced mass, E = k 2 /2µ, and k = µv is the relative momentum.The wave-function has an asymptotic form at large r as Here, f (k, θ) denotes the scattering amplitude, and the differential cross section for the scattering process is given by dσ The scattering state is decomposed into its partial-wave contributions (with multipole ℓ) as and the one-dimensional Schrödinger equation for the radial wave-function R k,ℓ is written as Since the potential term rapidly vanishes as r → ∞, the asymptotic form of the radial wavefunction at large r is where the scattering phase shift δ ℓ is real.We choose the normalization of δ ℓ to be δ ℓ (k) → 0 as k 2 → ∞ since the potential term is negligible at high velocity.From the asymptotic form of ψ k (x) − e ikr cos θ , we obtain the partial-wave decomposition of the scattering amplitude as follows: The elastic scattering cross section is decomposed into its partial-wave contributions.
The radial wave-function R k,ℓ (r) for the two-body DM state also determines an important effect for the DM annihilation, χχ → ϕϕ, with low velocity, namely the Sommerfeld enhancement factor.The annihilation occurs near r = 0, and the wave-function distorted from plane-wave near the origin significantly changes the size of the annihilation cross section.The enhancement factor with multipole ℓ, S ℓ , is given by the ratio of wave-functions with and without the potential at the origin [24,25] Here, R k,ℓ denotes the radial wave-function without the potential term.Now, we introduce numerical methods of computing the phase shift δ ℓ and the Sommerfeld enhancement factor S ℓ .Since there is no analytic solution for the Schrödinger equation with the Yukawa potential, we numerically solve the equation following Ref.[26] for the numerical method for obtaining the phase shift.We solve the one-dimensional Schrödinger equation Eq. ( 5) in a certain domain r i ≤ r ≤ r m , covering the one where the potential term dominates the equation compared to the kinetic term and the angular-momentum term.Meanwhile, the asymptotic solution of the equation at large r is given by where j ℓ (n ℓ ) is the spherical Bessel (Neumann) function.Matching the numerical solution to the asymptotic solution at r = r m , we obtain the phase shift δ ℓ : The normalization of the wave-function is irrelevant to determine the phase shift, then we take a simple initial condition for numerically solving the Schrödinger equation Eq. ( 5).The angularmomentum term dominates Eq. ( 5) near the origin, and hence χ k,ℓ (r i ) ∝ (kr i ) ℓ+1 .Therefore, we use χ k,ℓ ≡ A ℓ χ k,ℓ instead of χ k,ℓ (r) for the numerical computation, which satisfies an initial condition χ k,ℓ (r i ) = 1 and χ ′ k,ℓ (r i ) = (ℓ + 1)/r i .We numerically compute the Sommerfeld enhancement factor by using the same wavefunction χ k,ℓ .Meanwhile, near the origin, the wave-function χ k,ℓ (r) asymptotically behaves as Here, B ℓ represents the distortion of the wave-function near the origin relative to the plane-wave solution due to the presence of the potential.In other words, the Sommerfeld enhancement factor is given by S ℓ = |B ℓ | 2 .Therefore, matching the numerical solution to the asymptotic solution at r = r i (where χ k,ℓ (r i ) = 1), we obtain the Sommerfeld enhancement factor as Here, we determine the normalization A ℓ by using the phase shift given by Eq. ( 11) and the numerical solution χ k,ℓ (r) at r = r m , as We show the resonant self-scattering s-wave cross section and Sommerfeld enhancement factor for s-wave annihilation process in Fig. 1.Both the scattering cross section and the Sommerfeld enhancement factor resonantly enhance in low velocity (i.e.low a), and we can find that the resonant points for both quantities correspond with each other.In the following sections, we clarify the underlying reason why these resonant points are closely associated with each other.

Watson's Theorem and Omnès Solution
Now, we discuss a relation between the self-scattering phase shift and the Sommerfeld enhancement factor for the DM annihilation.A key to demonstrate the relation is the Watson's theorem [23]: the discontinuity of an amplitude is determined by the scattering phase shift of the initial (final) state in the presence of elastic scattering of the states.
In-and out-states are defined as eigenstates of Hamiltonian of interacting system, H = H 0 + V with free Hamiltonian H 0 and the interaction term V .These states become a free state Φ for V → 0, and are formally expressed by the Lippmann-Schwinger equation where ϵ > 0, and ± indicates in-state (+) and out-state (−).α collectively denotes labels of the states, such as momenta, spins, and so on.The DM annihilation with partial-wave ℓ is induced via an operator Θ χ , and we introduce the amplitude of DM annihilation as where |Ψ + ℓ ⟩ denotes an in-state for the initial two-body system (with partial-wave ℓ).There is a branch cut along real k 2 axis for k 2 > 0, and there may be poles for k 2 < 0 corresponding to bound states in the first Riemann sheet of the k2 -plane as shown in Fig. 2.There may also be virtual levels, resonances (unstable states), and anti-resonance in the second Riemann sheet.The amplitude Γ ℓ is a real function when the system has the time inversion invariance.
When the system has time-inversion invariance, by inserting the complete set of out states, we can rewrite the amplitude as follows. 2 Here, e 2iδ ℓ originates from the S-matrix element: and δ ℓ is the phase shift of the χχ → χχ scattering process.This indicates that the discontinuity of the amplitude is given by the scattering phase δ ℓ , which is known as the Watson's theorem [23].The analytic solution of Eq. ( 17) is known as the Omnès solution [27,28].The form of the solution is given by Here, F ℓ (k 2 ) is a real rational function, which is responsible for the pole structure of Γ ℓ (k 2 ) in the first Riemann sheet in Fig. 2, namely the bound states. 3b denotes the label of bound states.In the presence of bound states in the χχ → χχ scattering, F ℓ (k 2 ) has poles along the positive imaginary axis on the complex k plane.For instance, when there is a pole at k = iκ b .We will discuss the numerator after introducing the Levinson's theorem.The factor e ω ℓ (k 2 )+iδ ℓ (k) , called as the Omnès function, reproduces the branch cut of the amplitude. 4he Levinson's theorem [29] relates the behavior of the scattering phase shift to the number of bound states.We will introduce an intuitive proof of the theorem in Appendix B. The theorem states that where n ℓ denotes the number of the bound states, and N = 1 only for s-wave scattering just on the zero-energy resonance and N = 0 for other situations (other partial-waves or off the zero-energy resonance).As we will see in detail in the next section, there are zero-energy states on the resonances, which we refer to the zero-energy resonances: they are non-normalizable (i.e., virtual levels) for s-wave5 , while they are normalizable (i.e., bound states) for ℓ ≥ 1. n ℓ counts only bound states, but not virtual levels.δ ℓ (0)/π is closely related to the number of bound states due to our normalization of δ ℓ , δ ℓ (∞) = 0. We determine the numerator of F ℓ (k 2 ) as follows.We choose the normalization of the amplitude ignoring the potential term to be Γ ℓ (k 2 ) at high velocity.This determines the behavior of F ℓ (k 2 ) in the limit of k 2 → ∞. ω ℓ (k 2 ) goes to zero in this limit since we take the normalization δ ℓ (∞) = 0, and hence the function Meanwhile, according to the Levinson's theorem [29], the phase shift at small k relates with the number of the bound states.The ω ℓ (k 2 ) logarithmically diverges at small k, ω ℓ (k 2 ) ≃ π −1 δ ℓ (0) ln(Λ 2 /k 2 ), since the integral ( 19) is dominated at p 2 = k 2 .Λ characterizes the maximum momentum up to which the phase shift is regarded to be constant.Hence, the Omnès function behaves as at small k with the number of bound states n ℓ .Since the Sommerfeld enhancement factor should be constant at small k except for the zero-energy resonance (we will discuss later), the functional form of ) n ℓ at small k. 6 Hence, we obtain the functional form of F ℓ (k 2 ) as in Eq. (19).Combining the behavior of the Omnès function at small k, we have a constant Γ ℓ (k 2 ) at small k except for the zero-energy resonance (κ b = 0).For our formulation, we need to know the number of the bound states (i.e., the low-energy behavior of the phase shift) and the position of the bound states κ b .We discuss them in detail by introducing the effective range theory in the next section.Finally, the Sommerfeld enhancement factor is obtained as a ratio of the velocity-weighted annihilation cross sections with and without potential, and hence it may be rewritten as the ratio of the squared amplitudes with and without potential.
We may compute S ℓ for a given phase shift by the use of this formula instead of the standard computational method of the Sommerfeld enhancement factor given by Eq. ( 9).In Section 5, we compare the factors computed in these two ways.

Phase Shift and Bound States
Our formulation requires to know the number of bound states and corresponding poles κ b in order to determine F ℓ (k 2 ) as introduced in the previous section.The number of bound states is closely related to the low-energy scattering phase shift via the Levinson's theorem.It is a key for our formulation to know the low-energy behavior of the phase shift and the poles corresponding to the bound states.The low-energy scattering of two-body system is approximately described by the effective range theory [30,31], in which only two model parameters describe the two-body scattering.We determine the poles of the annihilation amplitudes by the use of the effective range theory as far as the poles corresponds to the shallow bound states, as we will see later.
We also use the effective range theory in order to find the reference points for our analysis.
It is a key point for the effective range theory that the wave-function of the two-body system must be analytic at k = 0.This analytic property determines the low-energy behavior of the phase shift, and hence the phase shift δ ℓ for multipole ℓ is expanded in the low-energy as Here, a ℓ and r e,ℓ are called the scattering length and the effective range, respectively.This expansion is valid as far as |k| < |r e,ℓ | −1 , otherwise the k 4 and higher terms are not negligible and it is beyond the effective range theory.This formalism has been recently utilized as a model-independent approach in order to approximate the velocity-dependent self-scattering of DM in the context of self-interacting DM [32] (see Refs. [22,[33][34][35] for concrete models).
Using the effective range theory, we find poles of the scattering amplitude given by Eq. ( 7).The poles of the s-wave amplitude are found along the imaginary axis of the complex k-plane as far as 2r e,0 < a 0 .One of the poles that is closer to the real axis, denoted by k = iκ b , is given by b /2µ is negative, and hence there exists a bound state (a 0 > 0) or a virtual level (a 0 < 0; a bound state with non-normalizable wave-function).These states become a zero-energy state as |a 0 | → ∞ (κ b → 0), and we refer to these points such that |a ℓ | → ∞ as the zero-energy resonances.For the higher partial-waves (ℓ ≥ 1), the poles appear at κ 2 b ≃ −2r 2ℓ−1 e,ℓ /a 2ℓ+1 ℓ , and hence there exist both a bound state and a virtual level (for a ℓ > 0 and r e,ℓ < 0, or for a ℓ < 0 and r e,ℓ > 0).The phase shift has interesting behavior for the effective-range theory parameters in a ℓ < 0 and r e,ℓ < 0. For such parameter range, there is a quick change of the phase shift by π at a certain momentum, k 2 = 2r 2ℓ−1 e,ℓ /a 2ℓ+1 ℓ , and the phase shift goes through π/2 (or a half-integer multiple of π).This indicates that the cross section has a peak at the momentum, namely there exists a Breit-Wigner type resonance.
We show the scattering length and the effective range in Fig. 3 for s-wave (top panels) and for p-wave (bottom panels).We mark several reference points in the left panels, which we will focus on in Section 5.These reference points are classified into three types: points below the first zero-energy resonance (▼ and ▲), on the zero-energy resonance (■), and points above the zero-energy resonance ( and ♦).Since the scattering length a ℓ diverges on the zero-energy resonance, we put marks ■ (b = 1.680 for s-wave and b = 9.082 for p-wave) on b-axis in the figure.Meanwhile, the scattering length vanishes, a ℓ = 0, at certain points depicted as ♦ in the figure (b = 4.466 for s-wave and b = 10.6997 for p-wave), and we call the points as the (first) anti-resonance.The effective range theory is no longer valid at these points since 1/a ℓ is indefinite and 1/|r e,ℓ | → 0. For both partial-waves, a ℓ is negative for b below the zero-energy resonance, while a ℓ is positive for b above the zero-energy resonance up to the anti-resonance.The effective range r e,ℓ is always positive for s-wave, while r e,ℓ takes both positive and negative values for p-wave.We use different colors for depicting connected lines for r e,1 in the bottomright panel.On the point ▲ for p-wave, since both a 1 and r e,1 are negative, there exists a Breit-Wigner type resonance.Now, we discuss the phase shift and the cross section for s-and p-wave scattering at points marked as ▼: b = 1.0 for s-wave and b = 6.0 for p-wave.For s-wave scattering, the scattering length is negative (a 0 < 0) and the effective range is positive (r e,0 > 0).The effective range theory Eq. ( 23) tells us that δ 0 → +0 in two limits k → 0 and large k and that δ 0 is maximized at k ≃ (−a 0 r e,0 ) −1/2 .We show the behavior of the phase shift and the cross section in Fig. 4 as gray-solid lines.We also show the scattering cross section (in unit of m 2 ϕ ) in the right panel.The cross section approaches to a constant value σ ∼ 4πa 2 0 in small k limit.Meanwhile, for p-wave scattering, the scattering length is negative while the effective range is positive, and δ 1 is maximized at k ≃ (−r e,1 /a 3 1 ) 1/2 .The gray-dashed lines in Fig. 4 show the behavior of the phase shift and the cross section for b = 6.0 of p-wave scattering.The phase shift scales as cot δ 1 ≃ (−a 1 k) −3 at small k, and hence the cross section scales as k 4 .
Fig. 4 also shows the phase shift and the cross section for the first anti-resonance points marked as ♦: b = 4.466 for s-wave and b = 10.6997 for p-wave.We note again that the effective range theory is no longer valid at these points.For s-wave scattering, the phase shift approaches to π at small k, shown as a purple-solid line in the left panel of Fig. 4.More precisely, in the small k limit, we find that the phase shift behaves as δ 0 ≃ π − b 3 0 k 3 with b 0 being a length parameter, and hence the cross section scales as k 4 at small k.Meanwhile, p-wave scattering, we find that the phase shift behaves as δ 1 ≃ π − b 5  1 k 5 with b 1 being a length parameter, shown as a purple-dashed line in the left panel of Fig. 4. The cross section scales as k 8 from the k-dependence of the phase shift at small k.As a generalization, for the first anti-resonance of multipole ℓ, we expect that the phase shift behaves as δ ℓ ≃ π − (b ℓ k) 2ℓ+3 at small k with b ℓ being a length parameter.Therefore, the effective range theory does not work since k 2ℓ+1 cot δ ℓ ∼ b −2ℓ−3 ℓ k −2 , and the scattering cross section scales as k 4ℓ+4 at small k.In Fig. 5, we show the phase shift and the scattering cross section for s-wave (top panels) and p-wave (bottom panels) near the first zero-energy resonance.Below the zero-energy resonance, the phase shift (depicted as green lines) goes to zero in two limits k → 0 and large k, and hence the phase shift takes a maximum value at a certain k.Therefore, the cross section (depicted as green lines) approaches to a constant value σ ∼ 4πa 2 0 at small k for the s-wave scattering, and the cross section scales as k 4 at small k for the p-wave scattering.The p-wave scattering cross section has a spiky behavior at k ∼ 10 −2 m ϕ for b = 9.081.Since both scattering length and effective range are negative for this b, the phase shift takes π/2 at k = (2r e,1 /a 3 1 ) 1/2 according to the effective range theory.Therefore, the scattering cross section takes its maximum value σ ≃ 3 × 4π/k 2 at this point.This behavior is also understood by the use of the resonant scattering cross section with running decay width [32,36], with . We find a relation between the effective-range theory parameters and the running decay width by matching the cross sections Eq. ( 8) and Eq. ( 25) near the resonance E ≃ E R as follows: The cross section Eq. ( 25) behaves as σ ℓ ≃ 4π(2ℓ + 1)a 4ℓ+2 ℓ k 4ℓ at low energy, which agrees with the low-energy behavior of the cross section Eq. ( 8).We confirm that the cross section Eq. ( 25) is in agreement with the green line in the right-bottom panel of Fig. 5 as far as k ≲ |r e,1 | −1 .
The phase shift goes to π at small k for b above the first zero-energy resonance for s-and p-waves, which is depicted as cyan lines in Fig. 5.More precisely, the phase shift behaves as δ 0 ≃ π − a 0 k for s-wave, and δ 1 ≃ π − a 3  1 k 3 for p-wave.There is a bound state; this is because the scattering length is positive at this point for s-wave, and a 1 > 0 and r e,1 < 0 for p-wave.This also confirms the Levinson's theorem (with n ℓ = 1), as we will introduce in the next section.The scattering cross section approaches to the constant value, σ ≃ 4πa 2 0 , for s-wave, and the scattering cross section scales as k 4 for p-wave.
Meanwhile, the phase shift approaches to π/2 for s-wave and π for p-wave on the first zeroenergy resonance (depicted as red lines in Fig. 5).For these points, the low-energy behavior of the phase shift is controlled only by the effective range since the scattering length diverges.Namely, at small k, the phase shift behaves as 2δ 0 ≃ π−r e,0 k for s-wave and as δ 1 ≃ π−2r e,1 k for p-wave.The scattering cross section scales as k −2 for s-wave, and the cross section approaches to a constant value for p-wave.

Comparison of Sommerfeld Enhancement Factors
Let us now compare the Sommerfeld enhancement factors computed in two ways: directly computed via Eq.( 13) and computed using the Omnès solution [Eq.(22)] for a given phase shift Eq. ( 11).We consider s-and p-waves here, but we also give a discussion of the first zero-energy resonance for d-wave in Appendix C. First, we discuss the reference points far from the first zero-energy resonance.Concerning the points marked as ▼ in Fig. 3, the function F ℓ (k 2 ) = 1 due to the absence of bound states.Since p ∈ [(−1/a ℓ ) , 1/r e,ℓ ] dominates the integral given by Eq. ( 19), ω ℓ (k 2 ) goes to a non-zero finite value at small k, and the Omnès solution Γ ℓ also goes to a finite value in the limit.As for the anti-resonance points marked as ♦ in Fig. 3, the phase shift approaches to π at small k.The small-k behavior of ω ℓ (k 2 ) is determined by the phase shift at low energy, and the Omnès function behaves as, in the small k 2 limit with a mass-dimension one parameter Λ as shown in Eq. ( 21).Λ characterizes the maximum scale where the phase shift is considered to be constant.Λ ≃ b −1 ℓ since we focus on the anti-resonance points here. 7Meanwhile, since there is a bound state for the parameter, F ℓ (k 2 ) is proportional to k 2 .Hence, the Omnès solution Γ ℓ goes to a finite value in the limit as required by the partial-wave unitarity.
Fig. 6 shows the Sommerfeld enhancement factors computed in two different ways: S ℓ computed directly via Eq.( 13) shown as black-solid lines, and S ℓ computed with the Omnès function shown as color-dashed lines.We choose the position of a bound state to be κ b = 1.32m ϕ for b = 4.466 (s-wave) and κ b = 0.66m ϕ for b = 10.6997(p-wave).We find that two S ℓ are in agreement with each other within our numerical accuracy of at most 10 % even though the computational methods are different.
Next, we discuss the Sommerfeld enhancement factors for b near the first zero-energy resonances, which are marked as ▲ , ■ , and in Fig. 3.The number of bound states is zero for the parameter b below the zero-energy resonance and on the zero-energy resonance only for s-wave, and hence F ℓ (k 2 ) = 1.Meanwhile, there is a bound state for the parameter b above the zeroenergy resonance and on the zero-energy resonance for the partial-wave with multipole ℓ ≥ 1.The rational function is given by ) for these parameters.As discussed in Section 4, we can determine the pole κ b by the use of the effective range theory as far as κ b ≲ |r e,ℓ | −1 .In other words, the pole determination by the use of the effective range theory is valid only when it is a shallow bound state.
We compare the Sommerfeld enhancement factors near the first zero-energy resonance, computed in two different ways, in Fig. 7.These factors are in agreement with each other within our numerical accuracy of about 10 %.For p-wave Sommerfeld enhancement factor with b = 9.081, a spiky behavior appears at k/m ϕ ≃ 10 −2 due to a resonance as with the scattering cross section.As mentioned in the previous section, this peak is correctly reproduced in our method even though we incorporate only the bound states in the rational function F ℓ (k 2 ).In other words, the Omnès function is responsible for reproducing not only the branch cut of the amplitude but also the resonances that are closer to the physical axis in the second As in Eq. ( 21), we determine F ℓ (k 2 ) so that S ℓ is constant at small k except for the zeroenergy resonance.On the other hand, at the first zero-energy resonance, the factor diverges at small k as S 0 ∼ k −2 and S ℓ ∼ k −4 for ℓ ≥ 1 in our formulation.For the first zero-energy resonance of s-wave, there is only a virtual state, and hence F 0 (k 2 ) = 1.Meanwhile, for the first zero-energy resonance of higher partial-wave (ℓ ≥ 1), there exists a zero-energy bound state with κ 2 b = 0, and hence F ℓ (k 2 ) = 1 again.From the Levinson's theorem, we find the phase shift at k = 0 with our normalization δ ℓ (∞) = 0, δ ℓ (0) = π/2 for the first zero-energy resonance of the s-wave and δ ℓ (0) = π for the first zero-energy resonance of the p-wave and higher partial-wave.ω ℓ (k 2 ) logarithmically diverges near k 2 = 0, and the Omnès function for the first zero-energy resonance at small k results in where Λ ∼ r −1 e,ℓ .Combining them, we get the small-k behavior of the Sommerfeld enhancement factor at the zero-energy resonance as shown above.Fig. 7 shows that our formula reproduces the small-k behavior of the factor via the direct computation.Since we choose the function F ℓ (k 2 ) so that the Sommerfeld enhancement factor is constant at small k except for the zeroenergy resonances, we obtain the constant value for the factor for the parameter near the first zero-energy resonance.Below the first zero-energy resonance, F ℓ (k 2 ) = 1 and the Omnès function is constant due to δ ℓ (0) = 0.For b above the first zero-energy resonance, ) and the Omnès function scales as e ω ℓ (k 2 ) ≃ Λ 2 /k 2 .As with the case on the first zero-energy resonance, our formula reproduces the small-k behavior of the factor via the direct computation.In particular, our choice of κ b by the use of effective range theory well reproduces the direct computation for the parameter above the first zero-energy resonance. 8efore closing this section, we comment on the partial-wave unitarity of the Sommerfeldenhanced cross section on the zero-energy resonance.Since the hard annihilation cross section with multipole ℓ scales as (σ ann,ℓ v) 0 ∝ v 2ℓ at small velocity, the velocity-weighted cross section σ ann,ℓ v scales as v −2 for s-wave and v 2ℓ−4 for higher partial-waves ℓ ≥ 1.Therefore, the annihilation cross section violates the partial-wave unitarity σ ann v ≤ 4π(2ℓ + 1)/m 2 χ v at small velocity for s-and p-waves.The annihilation cross section should be regularized for these partial waves to restore the partial-wave unitarity by taking into account the short-range interaction [37].

Conclusion
In this study, we have addressed the correlation between the resonant DM self-scattering and the large Sommerfeld enhancement factor for the DM annihilation.We have directly formulated the relation between the scattering phase shift and the Sommerfeld enhancement factor by the use of Watson's theorem.This is one of powerful theorems of quantum-mechanical scattering process: the discontinuity of the annihilation amplitude is determined by the scattering phase shift.The form of the amplitude is given by Omnès solution, which reproduces the pole structure and the branch cut of the amplitude.We have computed the Sommerfeld enhancement factor by the use of Omnès solution for a given phase shift for attractive Yukawa potential mediated by a light mediator.Our formulation has accurately reproduced the Sommerfeld enhancement factor, which is computed directly by solving the Schrödinger equation, with suitable choices of a function F ℓ (k 2 ) responsible for the pole structure of the amplitude.Moreover, we can use the effective range theory to determine the poles κ b when the self-scattering is close to the zero-energy resonance.If the velocity dependence of the DM self-scattering cross section is confirmed by observations, we can also determine the Sommerfeld enhancement factor by the use of our formalism, which is important for DM indirect searches.
The work of T. K. is supported in part by the National Science Foundation of China under Grant Nos.11675002, 11635001, 11725520, 12235001, and 12250410248.

A Hulthén Potential
In this appendix, we consider a two-body system under the attractive Hulthén potential This potential approximates the Yukawa potential: this scales as 1/r at short distances, and the force mediated by the potential is screened at large distances.m * denotes the screening mass, and hence we assume that m * is related to m ϕ as m * = κm ϕ with an O(1) coefficient κ. 9 The Schrödinger equation with the potential is analytically solvable for s-wave.For the higher partial-waves (ℓ ≥ 1), however modification of the centrifugal term allows us to find the analytical solution [38]: The modified centrifugal term reproduces ℓ(ℓ + 1)/r 2 at short distance m * r ≪ 1.
We can find the analytic form of the phase shift under the Hulthén potential.Defining a ≡ v/2α and c ≡ αm χ /m * and introducing a new variable t ≡ 1 − e −m * r , the analytic solution of the radial wave-function is given by [25] R k,ℓ (r) = t ℓ+1 2r e −ikr Γ(2ℓ + 2) where 2 F 1 denotes the hypergeometric function, and λ ℓ ± ≡ 1 + ℓ + iac ± √ c − a 2 c 2 , and Γ(x) is the gamma function.In the large distance limit (r → ∞, i.e., t → 1), 2 F 1 (λ − , λ + ; 2ℓ + 2; t) behaves as 2 F 1 (λ − , λ + ; 2ℓ + 2; t) Hence, the radial wave-function χ k,ℓ (r) ≡ rR k,ℓ (r) has the asymptotic behavior as 9 The coefficient κ is chosen to be κ = 2ζ(3) = 1.55 . . .by equating the Born scattering cross section with Hulthén and Yukawa potentials [26].As for the Sommerfeld enhancement factor, a similar prescription provides κ = ζ(2) = π 2 /6 = 1.64 . . . .However, the positions of the resonant enhancement for these processes should be agree with each other as seen in the case of Yukawa potential, and hence the coefficient must be identical.Matching the asymptotic behavior of χ k,ℓ (r) to the asymptotic form of the radial wave-function, the phase shift is analytically obtained.The phase shift is computed as the difference of the phases with and without the potential term.
Here, the second term is in agreement with +ℓπ/2 in the limit of k → ∞, which is an inherent phase for the scattering with partial-wave ℓ in the absence of the potential term when the centrifugal term is not modified.This s-wave phase shift gives the scattering length a 0 and the effective range r e,0 as Here, ψ (n) (x) is the polygamma function of order n, and γ E = 0.577 . . .denotes the Euler-Mascheroni constant.The scattering length diverges at c = n 2 with integers n.We illustrate the behavior of a 0 and r e,0 as a function of c in Fig. 8.For the p-and higher partial-waves scattering, we do not obtain the scattering length and the effective range that correctly reproduce those for Yukawa potential due to the modification of the centrifugal term of the Schrödinger equation.Indeed, at small k, the phase shift δ ℓ is expanded as and hence k 2ℓ+1 cot δ ℓ starts from k 2ℓ terms at small k.Specifically, we obtain the effective-range theory parameters for the higher partial-waves as, are numerically in agreement with each other with an accuracy of a few per mille.Compared to the Yukawa potential, we have the analytic formulae for the phase shift and the Sommerfeld enhancement factor, and hence these factors are identical in a high precision.
The behavior of the phase shift and the cross section is understood in a similar way to the Yukawa potential.For all of the reference points, we find good agreement of S ℓ given by the Omnès solution Eq. ( 19) (colored-dashed lines) with the analytic form (black-solid lines).We would emphasize that, for the Hulthén potential, we can find the pole position even for the parameters far from the first zero-energy resonance, in particular the first anti-resonance, thanks to the analytic form of the wave-function.The bound-state solution should be χ k,ℓ ∝ e −κ b r when we take k = iκ b , and thus a coefficient of the second term of Eq. ( 33) vanishes.We find the pole to be κ b = (c − 1)m * /2 for s-wave and κ b = (c − 4)m * /4 for p-wave above the first zero-energy resonance.Even though the effective-range theory parameters are not agree with that of the Yukawa potential for p-and higher-waves, the Sommerfeld enhancement factors for p-wave are accurately in agreement with each other (the right panel of Fig. 10).
In this appendix, we numerically investigate the relation between the phase shift and the Sommerfeld enhancement factor under the Hulthén potential.This indicate that one may have an analytic relation of the gamma function providing a relation between Eq. ( 34) and Eq. ( 41).However, we could not find such relation in the literature.

B Levinson's Theorem
As discussed in Section 3, there is a theorem that relates the behavior of the phase shift and the number of bound states.The original article [29] provides a rigorous proof of the theorem.We give a naïve sketch of a proof for a simple system, which is given in Refs.[41,42].We assume that the system has a central-force potential and that the system is enclosed in a large sphere of radius R. As far as the potential rapidly vanishes at large radius r, the scattering state for given angular momentum ℓ is approximately proportional to sin(kr − ℓπ/2 + δ ℓ (k)).
The boundary condition at r = R requires that the wave-function vanishes and the scattering state has a discretized momentum k n (n denoting an arbitrary integer that gives positive k n R).
Let n ′ be the minimum integer that gives a positive k n ′ R. For a given momentum k n , the difference of the phase shift is written as Here, N ℓ (k n ) denotes the number of the scattering states with the momentum in the range of In the absence of the potential, the phase shift vanishes and (k n − k n ′ )R/π gives the number of the scattering states with the momentum in a range of k n ′ ≤ k ≤ k n .Therefore, the number of the scattering states can change due to the presence of the potential, and its variation is given by Here, we take k n ′ → 0 for a sufficiently large radius R. As k n → ∞, this gives the total change of the number of the scattering states.Meanwhile, the total number of states is unchanged as far as we gradually turn on the potential.Some of the scattering states with the energy E > 0 for V = 0 becomes the bound state with the negative energy in the presence of V .Therefore, a portion of the scattering states is converted into the bound state, and the number of the bound state is given by This is known as the Levinson's theorem that gives a relation between the behavior of the phase shift and the number of bound states.We note that there is a subtlety related to zero-energy bound states and zero-energy virtual levels.As shown in Eq. ( 20), there is an extra term only when we consider the s-wave zero-energy resonance.

C d-wave case
In the text, we discuss the agreement of the Sommerfeld enhancement factors for s-wave and p-wave in two different computational methods.We also demonstrate that the use of the Omnès solution numerically works even for the d-wave case.We show the phase shift δ 2 (k) and the cross section in the unit of m ϕ on the first zero-energy resonance (b = 21.895) in Fig. 11.The phase shift scales as cot δ 2 ≃ (−r e,2 k) −3 at small k, and hence the cross section scales as k 4 at small k.We find that the phase shift crosses π around k ∼ m ϕ , and hence the scattering cross section has dips around there.We do not find the physical reason for that.Fig. 12 shows the Sommerfeld enhancement factors computed in two different ways as discussed in the text.Since the phase shift approaches δ 2 → π at small k, the Omnès function is logarithmically diverges as e ω 2 (k 2 ) ≃ Λ 2 /k 2 with Λ ∼ r −1 e,2 .Even for the d-wave, we find that these factors in two ways are in agreement with each other.

Figure 1 :
Figure 1: Self-scattering cross section σ (in unit of m ϕ ) and the Sommerfeld enhancement factor S for s-wave processes in the Yukawa potential as a function of b.

Figure 3 :
Figure 3: (top): Scattering length a 0 and effective range r e,0 for s-wave scattering in the Yukawa potential as a function of b in unit of δ.Since the scattering length diverges at b = 1.680, we put the mark for the point on the b-axis.(bottom): Same as top panels, but for p-wave processes.We put the mark for the point b = 9.082 on the b-axis.

Figure 4 :
Figure 4: Phase shift and cross section under the Yukawa potential for s-wave (solid) and p-wave (dashed) processes.(left): Phase shift for given b as a function of k/m ϕ .(right): Cross section in unit of m ϕ for given b as a function of k/m ϕ .

Figure 5 :
Figure 5: Phase shift and cross section near the first zero-energy resonance.(top): Phase shifts of s-wave scattering (left) and p-wave scattering (right) for given b as a function of k/m ϕ .(bottom): Cross section of s-wave scattering (left) and p-wave scattering (right) in unit of m ϕ for given b as a function of k/m ϕ .The black-dashed line in the bottom-right panel shows the cross section σ ≃ 3 × 4π/k 2 .

Figure 6 :
Figure 6: Sommerfeld enhancement factors for given b as a function of k/m ϕ : (left) for s-wave, while (right) for p-wave.S from direct computation is shown as black-solid lines, while S from the phase shift and the Omnès function is as colored-dashed lines.

Figure 7 :
Figure 7: Same as Fig. 6, but for the parameter b near and on the first zero-energy resonance.Different line colors correspond to different choices of b for s-wave (p-wave): b = 1.680 (b = 9.082) for on the first zero-energy resonance, b = 1.678/b = 1.682 (b = 9.081/b = 9.083) for below/above the first zero-energy resonance.

Figure 8 :
Figure 8: Scattering length (left) and effective range (right) for scattering in the Hulthén potential as a function of c in unit of m * .Since the scattering length diverges at b = 1.0, we put the mark for the point on the c-axis.

Figure 9 :
Figure 9: Phase shift and cross section under the Hulthén potential.(left): Phase shift of s-wave scattering for given c as a function of k/m * .(right): Cross section in unit of m * for given c as a function of k/m * .

Figure 10 :
Figure 10: Sommerfeld enhancement factors for given c as a function of k/m * .S ℓ from the analytic form is shown as black-solid lines, while S ℓ from the phase shift and Watson's theorem is shown as colored-dashed lines.

Figure 11 :
Figure 11: Same as Fig. 4, but for d-wave phase shift and cross section on the first zero-energy resonance.

Figure 12 :
Figure 12: Same as Fig. 6, but for d-wave Sommerfeld enhancement factors on the first zeroenergy resonance as a function of k/m ϕ .