Calculations of the Sommerfeld Effect in a Unified Wave Function Framework

In this paper, we suggest a method of a complete calculation of the Sommerfeld effect in both the s-wave and p-wave annihilation situations. The basic idea is to uniformly consider the equal-time Beth-Salpeter wave functions of both the dark matter and the final state particles, then short-distance interactions can be parametrized as a cross-term considering the equations. Solving these equations will result in the complete formulas of the cross sections with the Sommerfeld effects.


I. INTRODUCTION
If the dark matter particles interact with themselves through exchanging some bosonic light mediator, the annihilation behaviour of a dark matter particle pair as a function of their relative velocity can be significantly altered. If the force that the mediator carries is attractive, then the annihilation cross section σv will be amplified. This is usually called the Sommerfeld enhancement in the literature (For the original work by A. Sommerfeld, see [1]. For some early applications in the dark matter, see [2][3][4][5][6][7][8][9][10][11][12][13][14][15]). Usually, it is the resonance effect of the "zero-energy bound state".
Traditional calculations of the Sommerfeld effects usually involve solving the Schrödinger equation only including the "long-distance effects" at first, and then calculating the perturbative "short-distance" cross section σv 0 . The final result of the cross section becomes σv = σv 0 × S(v), where S(v) is the boost factor. In the s-wave annihilation processes, S(v) is calculated to be |ψ( 0)| 2 , where the ψ( 0) is the zero-point wave function of a stationary scattering state. This is equivalent to calculating the Feynmann diagrams as shown in the Fig. 1. In this paper, we try to regard both these interactions in a unified framework of wave functions in which both the Sommerfeld effect and the annihilation processes are all calculated by solving the equations of wave function. This is equivalent to resumming all the "trail"-diagrams as in the Fig. 2. In the literature, the Ref. [16] had followed this idea by introducing a complex δ-function in the dark matter potential (Ref. [3,17] have also adopted this potential), indicating a short-distance coupling making the Hamiltonian non-Hermitian.
The δ-function in the quantum mechanics requires renormalization [18]. The Ref. [16] had also pointed out that the traditional calculation method can break the unitarity in some particular parameter space [19], and their calculations can cure this problem. This technique is effective in calculating the s-wave annihilation processes.
The p-wave annihilation for the indirect detection signals (sometimes with the Sommerfeld enhancement) have also been discussed in the literature [20][21][22][23]. Generalizing the method in the Ref. [16] to the p-wave case confronts with difficulties. δ function potential only interact with the l = 0 partial waves. If we simply add some complex potential, e.g., ∂z , that interacts with the l = 1 partial wave, since usually l = 1| ∂δ (3) ( r) ∂z |l = 0 = 0, such a term will disturb the l = 0 partial waves. This does not happen in some practical dark matter model. In this paper, we will make another approach. If we uniformly describe the annihilating dark matter and the final state particles in the form of Bethe-Salpeter wave functions, we can then parametrize and write down the short-distance interaction terms between them in the wave function equations. In this situation, an interaction term like ∂δ (3) ( r) ∂z couples with the different particle states, thus avoiding the problems described before.
The Sommerfeld effect arises in the non-relativistic limit of the dark matter, so the dark matter can be described by the simple Schrödinger equation. Usually, the annihilation products are usually light standard model (SM) particles compared with the dark matter particles, and can be regarded as massless particles. It is convenient to describe the finalstate particles with the massless Klein-Gordan equations, or d'Alembert equation with an interaction term. We will prove in our paper that the Beth-Salpeter wave function of the relative motion of a pair of particles in many cases satisfy the d'Alembert equation. We will solve these equations and extract the cross sections from the results. We will show the final expressions as well as the derivation processes in detail.

II. SOME GENERAL DISCUSSIONS
Let χ( r) be the dark matter's stationary wave function of the relative motion, and φ( r) be the stationary wave function of relative motion for the final state, where r is the relative displacement between the particles inside the pair. The short-distance scattering effects can be parametrized by an interaction term which is non-zero only at the origin. For an swave scattering process, only the δ-function arises. Considering the long-distance effects by exchanging some light mediators between the dark matter particles, the stationary equations of the dark matter-annihilation product system can be given by for s-wave interactions.
For a p-wave process, the derivative of a δ-function is utilized. In a practical model, the index of the derivative should be contracted with another index. In this paper, we consider a special case that a L = 1, S = 1, J = 0 dark matter pair annihilates into the L = 0, S = 0, J = 0 final states (e.g., χχ → Φ → f f , where χ is the majorana fermionic dark matter, Φ is an s-channel CP-even scalar mediator, and f is the final product). In the section V, we will prove that the L = 1, S = 1, J = 0 final states can also be described by the same equation. In this case, the index of the derivative of a δ-function is contracted with the triplet spin-index of the initial state. Extract the third component of the initial-state triplet wave functions, we acquire One can consider some other cases such as L = 1, S = 0, J = 0 dark matter annihilating to the L = 1, S = 0, J = 0 final states (e.g., χχ → Z → f f through a vector mediator Z ). We should note that all the differences during the calculations can finally be absorbed into the renormalization parameter k Λ2 , although we are not going to present the detailed derivation in this paper. In the above equations, µ is the reduced mass of the dark matter pair. u s,p are the coupling constants in the s-wave and p-wave scattering cases, respectively.
Because of the conservation of the energy during the annihilation processes, the total energy where E k is the total kinematic energy of the dark matter particle pair. V ( r) is the potential indicating the long-distance interaction between the dark matter particles.
Here, we only concern the spherically symmetric potential, so that V ( r) = V (r). Usually lim r→∞ V (r) = 0 and V (r) can have a pole with no more than order one at the origin.
Temporarily omitting the interaction terms, after separating the angular variables, the radial equations are given by where l 1,2 indicate the quantum number of the total orbital angular momentum of the χ and φ partial waves respectively. The definitions of the χ l 1 and φ l 2 are given by where m is the orbital angular momentum quantum number along the z-direction, Y lm (θ, φ) are the spherical harmonics, and θ, φ are the angles in the polar coordinate system. Now we consider the eqn. (3), and define k 2 = 2µE k , (3) becomes −χ l 1 (r) + l 1 (l 1 + 1)χ l 1 (r) + 2µV (r)χ l 1 (r) = k 2 χ l 1 (r).
For each k ∈ R and k ≥ 0, there are two linearly independent solutions u k (r), v k (r), with the two different asymptotic boundary conditions In order to calculate the cross section, we need to decompose the wave functions into the asymptotic incoming and outgoing wave functions e ∓i(kx+δ l 1 ) . Generally, the incoming and outgoing wave functions of a particular momentum k can be written in the form If V (r) is an attractive potential, there will be some further bound states with a scattered spectrum of E i = k 2 i < 0, where i is some index to distinct the different bound states. We still express these eigenfunctions as u l 1 (k i , r), and the normalization convention is defined The asymptotic behaviour in the r → 0 is the same as the k 2 > 0 cases, but the asymptotic behaviour in the r → ∞ is calculated to be According to the Sturm-Liouville theory, the {u l 1 (k, r)}, including all of the u l 1 (k i , r)'s, form a set of complete orthogonal basis on the function space defined on the r > 0 domain.
We should note that the v l 1 can also be expanded on the basis of u l 1 .
However, if we restore the above incoming and outgoing radial functions to the threedimensional wave function χ l 1 and φ l 1 , there will be a divergence at the origin, making them difficult to be coupled with the interaction term in the (1) and (2). Since {u l 1 (k, r)} is a set of complete orthogonal basis, and is regular at the origin, we can expand u l 1 ,in/out with the basis of u l 1 , and indirectly calculate the couplings at the origin with their aids. In order to calculate the inner product defined by the integration, consider two functions u 1 (r) and u 2 (r) with the eigenvalue k 1 and k 2 , satisfying −u 1 (r) + l 1 (l 1 + 1) r 2 u 1 (r) + 2µV (r)u 1 (r) = k 2 1 u 1 (r), −u 2 (r) + l 1 (l 1 + 1) (11)×u 2 -(12)×u 1 , and after a few steps, we acquire Now we are ready to calculate the ∞ 0 u l 1 ,in (k, r)u l 1 (k , r)dr and ∞ 0 u l 1 ,out (k, r)u l 1 (k , r)dr. Notice that if we add a small imaginary part ∓iλ (λ > 0) to the k when solving the (6), we will acquire the damping wave functions u l 1 ,in (k − iλ, r) ≈ u l 1 ,in (k, r)e −λr and u l 1 ,out (k + iλ, r) ≈ u l 1 ,out (k, r)e −λr .
Therefore, we acquire lim r→∞ u l 1 ,in, out (k ∓ iλ, r)u l 1 (k , r) − u l 1 ,in, out (k ∓ iλ, r)u l 1 (k , r) = 0. Finally, the inner products are given by If k = k , u l 1 ,in (k, r)u l 1 (k, r)−u l 1 ,in (k, r)u l 1 (k, r) and u l 1 ,out (k, r)u l 1 (k, r)−u l 1 ,out (k, r)u l 1 (k, r) become constants (these are called the "Wronskian") independent of r. From the asymptotic definition (7), we can acquire Here, A l 1 ,k is a real number according to our convention in (8). For simplicity, we can adjust the phases in the (7) to guarantee A l 1 ,k > 0. Then, Immediately, we acquire We can then calculate the normalization coefficient of the basis {u l 1 (k, r)}. Calculate , and take the limit λ → 0, we can prove that With this result, we can then expand u l 1 ,in, out on the basis {u l 1 (k, r)}. Considering both the contributions from the bound states and the continuous spectrum states, the result is then The φ l 2 's have the similar results, and the above discussions can be directly transferred to the φ l 2 .

A. No Sommerfeld Effect
In order for the comparison and preparation, we at first calculate the S-wave scattering process without Sommerfeld effect, i.e., V (r) = 0, and l 1 = l 2 = 0. In this case, u 0 (k, r) = sin(kr), v 0 (k, r) = cos(kr); The Fourier expansion then becomes Just like in the quantum mechanics case, we write the χ 0 and the φ 0 in the composition of the incoming, reflecting and scattering wave functions where due to the conservation of the energy. Therefore, To solve the (1), we also need to expand the three-dimensional δ-function, and the zero-point value of the wave function is also calculated to be Note that the integral in the (32, 33) is divergent, some renormalization scheme is re- The last step is because the second integrand is an even function. The first divergent integral can be renormalized by a cut-off k Λ , and the second integral is calculated to be ∓ iπk 2 . In this process, we have applied a different form of, but fundamentally equivalent renormalization scheme compared with the Ref. [18,24,25]. Therefore, Then, with a few steps, we can see that if (1) can be satisfied for the eigenvalues E k = k 2 1 2µ and E 2 t = k 2 2 . Solving (36), we have Notice that for a nonzero u s , C < 1, and for nonzero k Λ1,2 , a phase shift will be induced, too. The cross sections for χ → χ scattering and χ → φ annihilation is Although a matching with the quantum field theory perturbative calculations requires both the results from (38) and (39), we only write down the annihilation cross section explicitly because of our main topic, . (40) For the purpose of matching with the perturbative calculations of the quantum field theory, we also need to calculate the φ → φ cross sections. This is calculated to be For the brevity of this paper, we have skipped all of the lengthy but familiar derivation of the (41).

B. Sommerfeld Effect Calculations
If the Sommerfeld effect is considered, i.e., in the V (r) = 0 case, the A 0,k should be calculated. Traditionally, a scattering wave function in the V (r) potential ψ( r) with the asymptotic condition ψ( r) ∼ e ikz (r → ∞) is calculated, and its origin value |ψ k ( 0)| is acquired for further calculations. From the definitions in the (15), we learn that A 0,k is actually the the |ψ k ( 0)|. Therefore, the (29) should me modified to The δ-function in the χ-equation in the (1) can be written in the form (32) also changes the integral in (32) to be The integral can again be separated into two parts, The (45) has two problems. To integrate out the first term and compare it with the (34), we need to know the ψ k ( 0) for every value of k, and the second term should not be calculated by simply picking up the residue only near the real axis because the full analytical properties of ψ k ( 0) in the complex plane are unknown. However, the k = ±(k±iλ) residue contributions in the second term are proportional to i(1−C), while the all the other parts are proportional to (1+C) and can be attributed to one term t ψ k ( 0) ·k Λ . t can be calculated in any value of the cut off k Λ by principle. However, notice that if k Λ k 1 , the integrand in the high-energy area mainly contribute to the results, and notice that because lim k→∞ ψ k ( 0) = 1, in this case t ∼ 1. On the other hand, If k Λ µ, we will see the final result will be not sensitive to the detailed value of t. In the following discussions, we will still explicitly preserve t, and the (35) can be modified to After a similar calculation, the (36) now becomes Solving this equation, we have It can be proved easily that |C| ≤ 1 so the unitarity will never break up. and finally, σ ann, Sommerfeld = (49) To observe the dependence of σ ann on the |ψ k 1 ( 0)| 2 , we take the limit k Λ1 = k Λ2 = 0, which is reasonable in some special cases to be described, then It is now clear that when k 1 k 2 u 2 s µ 2 8π 2 , σ ann, Sommerfeld ≈ |ψ k 1 ( 0)| 2 σ ann if |ψ k 1 ( 0)| 2 is not so large. For extremely large |ψ k 1 ( 0)| 2 , σ ann, Sommerfeld then becomes saturated before reaching the unitarity bound.

C. Some Discussions
In this subsection, we will take some discussions on the renormalization scales. From the aspect of perturbative quantum field theory, The divergences arise from these two diagrams in the Fig. 3. In the quantum field theory, the t-and u-channel loop diagrams are also considered. In the wave function calculations, the t-and u-channel loop diagrams are omitted, but the diagrams in the Fig. 3 are connected from head to tail and are resummed.
The divergences induced by these diagrams are cancelled out by the two counter-terms in the Fig. 4. One can try to add the corresponding self-interaction terms in the (1), however, this is equivalent to adjusting the renormalization scale k Λ1 and k Λ2 .
It is easy to recognize that in the (1, 2), we did not introduce the elastic short-distance terms indicating the χ → χ and φ → φ processes. However, such elastic scattering effects can still arise due to the Fig. 3 loops together with the corresponding counter terms in the • Calculate the short-distance perturbative s-wave dark matter annihilating cross section, dark matter s-wave self-interaction cross section, and the final state selfinteraction cross section by the quantum field theory.
• Compare the perturbative results with the (38, 40, 41) to determine the u s , k Λ1 , k Λ2 parameters.
• Calculate the φ k 1 ( 0) as usual, and then determine the value of the t. Take all of these values to (50) to calculate the Sommerfeld corrected cross section.
To determine the t, we need to know what we are comparing when we are discussing the "boost factor". Let us investigate the meaning of the cross section without the "long-range" interaction, i.e., the "σv 0 " from the perturbative quantum field theory aspect. Usually, if we shut down the long-range interactions to calculate the "σv 0 ", all the divergences will structures. However, if we want to compare the two models with the same "renormalized parameters", which means that in this case, the "renormalized masses" of the Φ mediator remains unchanged, which is proportional to the cut-off scale k Λ1 in the effective quantum mechanics models, we immediately know that t = 1.
Another approach to determine the t can be found in the Ref. [16]. The authors of the Ref. [16] define the Sommerfeld factor by comparing the low-energy and the high-energy cross sections of exactly the same model. This prevent the problem of comparing two different models, and is actually equivalent to the t = 1.
If the counter terms "coincidently" cancels nearly all the divergence terms, the elastic scattering can be nearly prohibited. In this case, k Λ1,2 µ. This can be an approximation for the models such like the dark matter annihilates into some light mediators through t-channel diagrams. In this case, the dark matter's self-interactions are mainly due to the long-distance effects. Therefore, our discussions in the (51) are available in such cases.

IV. P-WAVE SCATTERING
Now, we are going to calculate the stationary scattering states of the (2). Since the derivative of a three-dimensional δ-function at the origin only has a l = 1 expansion along the θ − φ direction, and particularly, ∂δ (3) ( r) ∂z only has a l = 1, m = 0 expansion, (2) can only couple the p-wave in the χ-field with the s-wave φ-field. Therefore, l 1 = 1 and l 2 = 0 in this case.

A. No Sommerfeld Effect
Again, let us try V (r) = 0 at first. In this case, where j n (x) = (−x) n 1 x d dx n sin x x are the spherical Bessel functions. The generalized Fourier expansion of the incoming and outgoing wave functions are given by u 0 (k , r), Define and then, The derivative of the wave function χ( r) at the origin is calculated to be To separate the divergence, The first term is renormalized to be k 2 Λ , so ∂χ r ( r) The equations corresponding to the φ are similar to the s-wave cases, therefore we do not repeat them here.
Again, from the (15), A 1,k = |ψ k0 |. The expansions of the incoming and outgoing wave functions have changed to Here is the expansion of the ∂δ (3) ( r) ∂z , and the derivative of the wave function χ( r) is modified to . (68) After renormalizing the divergent part into tk Λ1 , we have Now, the algebraic equation is given by The solution of this equation becomes Correspondingly, σ ann, Sommerfeld = (72) 9 × 12k 1 k 2 π 7 |ψ k 1 0 | 2 u 2 p µ 2 The discussions of the renormalization is similar to the s-wave case, so we do not repeat it in this section. However, we should note that in a real dark matter model, loop diagrams in the Fig. 3 will also contribute to the |l = 1, m = 0 → |l = 1, m = ±1 amplitudes. This is done through the ∂δ (3) ( r) ∂z terms which should by principle appear in the (2). The cross section of such a channel is at least proportional to v 4 , where v is the relative velocity of the two dark matter particles, so we do not concern this suppressed effect in this paper. Again, if the counter term cancels most of the divergences, we have Again, it is obvious that the C in the (71) will be saturated before it reaches the unitarity bound.

V. WAVE FUNCTION OF THE MASSLESS FINAL STATE PARTICLE PAIR
In this section, we are going to patch a leak. The first equations of the (1-2) are Schrödinger equations and there have been sufficient discussions about this in the literature (For an example and an application in the dark matter, see Ref. [4]). However, the final state particles can usually be considered as massless in most of the dark matter models, and whether the "wave function of relative motion" of such a relativistic system should satisfy the d'Alembert equation remains a quest. From the aspect of quantum field theory, the stationary state |N is characterized by the Bethe-Salpeter wave function (For a review, see Ref. [26]. For the original paper, see Ref. [27]. For an application in the dark matter bound state, see Ref. [28]) defined by where N | is the corresponding n-particle stationary state, and ψ i (x i ) can be the same or different particle field operators. Generally, the time component x 0 1,2,...,n can be different. However, if all the interactions among the particles can be regarded as instantaneous interactions, and the effective vertices are independent on the q 0 , which is the time-like component of the transfer momentum, as just the case in our paper, the Bethe-Salpeter wave function can be reduced to the equal-time wave function ψ N (t, r 1 , r 2 , ..., r n ) = N |ψ 1 (t, r 1 )ψ 2 (t, r 2 )...ψ n (r, r n )|0 .
Solving these wave functions involve the Beth-Salpeter equation. However, in this paper we adopt another convenient method. We can directly take the derivation on t at both sides of the (75), and then replace the operators on the right-handed side by the equation of motions.
As for the p-wave scattering processes into fermionic particle pairs, the final state can be either S = 0 or S = 1. The majorana dark matter particle pair's p-wave annihilation through an s-channel scalar particle usually result in S = 0, or longitudinal S = 1 final states. The following discussions are mainly based on this situation. Another possibility is that the complex scalar dark matter annihilating into transverse S = 1 final states though an s-channel massive boson. We do not present the detailed derivation for brevity.
The massless fermionic final state particle pairs can be described by two Weyl spinors.
Let ψ 1 (t, r) and ψ 2 (t, r) be the two 2-dimensional Weyl spinors, each of them satisfy the equation of motion where σ i are the Pauli matrices, and O( r) is the interaction term containing the dark matter fields, e.g., for a s-wave interaction effective theory, O(t, r) can be gΦ 1 (t, r)Φ 2 (t, r), where Φ(t, r) is the dark matter's field operator, g is the coupling constant which can be complex, and for a p-wave interaction effective theory, O(t, r) = gΦ 1 (t, r)∂ i Φ 2 (t, r). The definition of the scalar-longitudinal-vector Bethe-Salpeter wave function is given by where |N not only contains the ψ 1 -ψ 2 state, the two-dark matter states Ψ 1 Ψ 2 are also included. Because we only care about the wave function of the relative motion of a particle pair, with the total momenta p tot = 0 at rest, we can learn that such a wave function only depends on the relative coordinate r 1 − r 2 , so we define This wave function can be decomposed into a vector part and a scalar part ψ w,2 (t, r) = V (t, r) · σ + A(t, r)I, where I is the 2 × 2 unitary matrix. Take the derivative the (8) according to the t parameter and consider the (76), we have Here, O(t, 0) = N |O(t, 0)|0 , and D w ( r) = 0|(iσ 2 )ψ * 2 (t, r)ψ T 2 (t, 0)(iσ 2 )|0 + 0|ψ 1 (t, r)ψ † 1 (t, 0)|0 . We should note that although we are calculating the (80) in the Heisenberg picture, we can still contract the operators of the same particle by the Wick's theorem since we have ignored all the n > 2 multi-particle states in this paper. Finally, D w is calculated to be The above derivation have applied the properties (iσ 2 )σ µT (iσ 2 ) =σ µ .
Then we can separate the different SO(3) representations in the (80), into spin-0 part A and spin-1 part V . We have Eliminating the V , we acquire an equation about A, Note that we are discussing a stationary state N |. According to the Heisenberg's equation, where H is the Hamiltonian, and E N , E 0 are the energy eigenvalue of the state |N and |0 respectively. Usually, the vacuum energy E 0 is renormalized to be zero, so which proved that the scalar part of the relative wave function of a massless fermionic particle pair does satisfy the d'Alembert equation with an interaction term.
As for the vector part V , from the (83), we can see that then where f ( r) is some vector field which is independent of t. Since we are only interested in stationary wave function with a positive energy E, there will be an unavoidable e iEt factor in the V (t, r). For a non-zero V (t, r), the only choice is that f ( r) = 0. Therefore, This means that V should be a longitudinal vector field.
Notice that A and V are not independent fields. Through (83) we can see that they can be derived to each other. If A is a s-wave (L = 0) field, from (83)  For simplicity, we can just choose the scalar A as the final-state wave function in the (1,2).
In order to describe the p-wave annihilation to the transverse S = 1 final states, we should construct the following Bethe-Salpeter wave function, These can again be decomposed into a scalar part A 11, or 22 (t, r) as well as a vector part field situation. Other situations should follow a similar process. Let φ 1 (t, r), φ 2 (t, r) be the two final state particle fields, and they satisfy the equation of motion The wave function is given by Again, we only discuss the wave function of relative motion so the wave function only depends on r 1 − r 2 , ψ s,2 (t, r 1 − r 2 ) = ψ s,2 (t, r 1 , r 2 ) = N |φ 1 (t, r 1 )φ 2 (t, r 2 )|0 .
Here we have proved that the wave function of the massless scalar particle pair satisfies d'Alembert equation.
Finally, we need to explain the different orders of the divergences between the diagramatic and wave-function calculations. For example, if we calculate the loops in the Fig. 3, we acquire a logarithmic divergence in the scalar case, and a quadratic result in the fermionic case. However, a wave-function based calculation gives a linear divergence. This is due to the following two reasons, • In the wave-function approach, we describe the "short-distance" interactions by some "point-like" δ-functions. This is usually not the truth, and will modify the orders of the divergences.
• The "ladder approximations" of the wave-function approach usually pick up some of the poles of the particle propagators, while ignores some others. Although these ignored poles does not affect the nearly on-shell performances described by the wave functions, they will actually contribute to the divergence orders.
Fortunately, all these differences can be absorbed to the renormalization scale k Λ2 , and since usually, the k 2 ∼ 2µ remains nearly unchanged when the dark matter relative velocity v 1, we also do not need to worry about the running couplings which behaves quite differently in the bosonic and fermionic cases. Therefore, our calculations in this paper are still reasonable.

VI. SUMMARY
In this paper, we have made the general discussions and calculations on the s-wave or the p-wave dark matter particle pair annihilating into a pair of s-wave massless particles.
Practically, one can compare the results from the perturbative quantum field theory with the information of dark matter annihilation, dark matter elastic scattering and final state particle elastic scattering informations with the formulas in the section III A, or IV A, to determine the coupling constant u, and the two renormalization scales k Λ1,2 , and then calculate the ψ k ( 0) and ψ k0 by following the formulas in the section III B, or IV B for a complete calculation. We have also proved that the wave function of the final state can actually be described by the wave function equation (1) and (2).