Systematizing the Effective Theory of Self-Interacting Dark Matter

If dark matter has strong self-interactions, future astrophysical and cosmological observations, together with a clearer understanding of baryonic feedback effects, might be used to extract the velocity dependence of the dark matter scattering rate. To interpret such data, we should understand what predictions for this quantity are made by various models of the underlying particle nature of dark matter. In this paper, we systematically compute this function for fermionic dark matter with light bosonic mediators of vector, scalar, axial vector, and pseudoscalar type. We do this by matching to the nonrelativistic effective theory of self-interacting dark matter and then computing the spin-averaged viscosity cross section nonperturbatively by solving the Schrodinger equation, thus accounting for any possible Sommerfeld enhancement of the low-velocity cross section. In the pseudoscalar case, this requires a coupled-channel analysis of different angular momentum modes. We find, contrary to some earlier analyses, that nonrelativistic effects only provide a significant enhancement for the cases of light scalar and vector mediators. Scattering from light pseudoscalar and axial vector mediators is well described by tree-level quantum field theory.


Introduction
The majority of mass in our universe is in the form of dark matter (DM), but the underlying nature of DM and its interactions remains elusive. Many attempts to understand the nature of dark matter rely on searching for its interactions with ordinary matter, through direct or indirect detection experiments or attempts to directly produce dark matter particles at colliders or fixed-target experiments. We may also be able to learn about the particle nature of dark matter if it has significant self-interactions, which can alter the astrophysical and cosmological signals of dark matter in ways that may be detectable [1,2]. In recent years, many aspects of the astrophysics and cosmology of self-interacting dark matter have been extensively studied, so we will of necessity refer to only a small fraction of the literature. An excellent recent review article with extensive references is [3].
Dark matter self-interactions have been posited as a possible explanation for a number of discrepancies on small scales between observation and the results of classic N -body dark matter simulations of standard ΛCDM cosmology, such as the core/cusp problem [4,5], the missing satellites problem [6,7], the too-big-to-fail problem [8], and the diversity problem [9,10]. A weakness of such arguments is that baryonic physics, such as supernovae or AGN activity, can alter the distribution of dark matter in galaxies and clusters in ways that are unaccounted for in dark matter-only simulations. Incorporating baryonic physics in simulations accurately is challenging, and may solve some or even all of the potential problems with ΛCDM on small scales [11][12][13][14][15][16].
In this paper, we do not take sides in this debate. We expect that in the future our understanding of dark matter and baryonic effects could be refined to allow for an averaged DM scattering rate to be extracted from data in systems on a variety of scales, from dwarf galaxies to clusters. Because the typical velocity of a dark matter particle is much lower in a dwarf galaxy than in a galaxy cluster, such observations could map out the velocity dependence of dark matter scattering, which encodes information about the underlying particle physics [17,18]. An attempt to do this with current data has been made in [19] and argued to fit a Yukawa potential with a light mediator. Again, one can debate whether these conclusions are robust to uncertainties in baryonic physics, but in any case an important message to take away from [19] is that observations may eventually tell us the quantity σv /m DM as a function of the rms velocity v.
With such measurements in hand, we would naturally want to answer the inverse problem: what underlying model of dark matter self-interactions produces the observed velocity dependence of the cross section? This question is nontrivial because dark matter is nonrelativistic, and in the nonrelativistic limit low-velocity scattering can be nonperturbatively enhanced by the Sommerfeld effect [20]. To compute the cross section in the low-velocity regime, we must match the relativistic QFT of interest onto a nonrelativistic effective theory [21,22] in which we can compute the cross section nonperturbatively by (numerically) solving the Schrödinger equation. In the dark matter context, the need for such calculations was first appreciated in the annihilation of heavy WIMPs [23][24][25]. More recently, it has been applied to scattering of dark matter through light mediators that generate an effective Yukawa potential [17,18,26,27]. However, relativistic QFTs can match to a much wider range of nonrelativistic effective interactions than a simple Yukawa potential. The classification of such interactions [28] has been used to study a variety of possible recoil spectra in direct detection experiments [29][30][31] (see also [32,33]). Most relevantly for our current work, it has also been used to classify dark matter self-interactions in [34], the "effective theory of self-interacting dark matter." This paper is closely related to that one, although some of our conclusions differ.

Goals of this work and relation to the previous literature
Our goal in this paper is to calculate the velocity dependence of dark matter self-interactions for the case of spin-1/2 dark matter interacting via a light boson, which may be a scalar, pseudoscalar, vector, or axial vector. In each case, we match to a nonrelativistic effective theory, then solve the Schrödinger equation numerically to obtain the velocity dependence of the cross section.
Our work differs from, and extends, earlier work on nonrelativistic dark matter scattering in several respects. The first is our matching procedure: we match the Born approximation to shortdistance scattering in the quantum mechanical effective theory to the tree-level perturbative QFT approximation to short-distance scattering. This provides a boundary condition at small radius, from which we can integrate outwards to solve the Schrödinger equation and capture long-distance effects of light mediators. The specifics of our matching procedure are inspired by earlier work [21,22], but the details are novel: our procedure is streamlined and easy to apply. A particular difference from some previous work on self-interacting dark matter is that we do not just consider the effective potential generated by t-channel exchange, but include the contact interactions arising in other channels. In the case of pseudoscalar mediators, this is crucial to obtain the correct cross sections at small velocities.
Our approach clarifies a number of issues related to pseudoscalar mediators. We sum over all angular momentum partial waves. While a simple Yukawa potential conserves orbital angular momentum, the exchange of a pseudoscalar or axial vector leads to interactions that couple different l modes (while conserving the total angular momentum j). Some earlier work has considered a toy model aimed at approximating pseudoscalar exchange without treating the coupled channels carefully. This toy model includes a singular 1/r 3 potential. We will show that this misses an important aspect of the physics, namely that the Born approximation for the pseudoscalar potential is non-singular at short distances. Furthermore, our matching procedure is robust to variations in the matching scale. Unlike some results in the literature, our cross section is entirely determined by the underlying QFT and does not depend on ad hoc constants introduced in the nonrelativistic effective theory.
In the end, we find that the correct matching of a weakly-coupled effective field theory with a light pseudoscalar or axial vector mediator leads to a nonrelativistic effective theory in which there is no enhancement of the cross section at low velocities. For these theories, unlike the case of light scalar or vector mediators, tree-level QFT is reliable. This simple result is in contrast with some earlier claims in the literature, for instance, in studies of annihilating DM with a pseudoscalar mediator [35,36].
We proceed as follows: in §2, we introduce the basic models of mediators in QFT. We review why the axial vector mediator is special, in that its couplings to dark matter must vanish as the mediator mass goes to zero. In §3, we describe our matching procedure, the way in which we set boundary conditions, and the process of extracting the S-matrix from numerical solutions of the Schrödinger equation. In §4, we present the numerical results of solving the Schrödinger equation and computing the cross section for the various models. In particular, we explain the absence of a Sommerfeld enhancement for a pseudoscalar mediator, arguing that a 1-loop calculation in the perturbative QFT also suggests that the effect should be absent. We explain the differences between our results and certain claims in the literature. We offer concluding remarks in §5.

SIDM and Our Examples
We consider a weakly coupled dark sector with spin-1/2 fermionic dark matter interacting via scalar or vector mediators. If these mediators are light and generate an attractive potential, then it is possible that Sommerfeld enhancement can significantly boost the scattering cross section. So, it is crucial to analyze the potentials generated by various renormalizable interactions.
The details of the calculation depend on whether the dark matter carries an approximate conserved charge. If dark matter is a Majorana fermion, then the χχ → χχ scattering process receives contributions from s-, t-, and u-channel diagrams. In the Majorana case, there is no vector interaction, i.e., the coupling A µ χ † σ µ χ takes the form of an axial vector interaction when packaged into a 4-component field. If dark matter carries an approximately conserved charge, we should consider a Dirac fermion, and the dark matter abundance may be primarily of one charge (asymmetric dark matter [37][38][39]) or it may contain particles of both charges (symmetric dark matter). Scattering in the asymmetric case, χχ → χχ, receives contributions from only t-and u-channel diagrams. In the symmetric case, there is additionally the process χχ → χχ, which receives s-and t-channel contributions but no u-channel contribution.
Our goal in this work is to clarify conceptual issues in matching theories of self-interacting DM to a nonrelativistic effective theory, and to understand in which cases a Sommerfeld enhancement is present. These aspects of the physics are not sensitive to the Majorana or Dirac nature of the fermion or to the symmetric or asymmetric nature of the dark matter population. Thus, for concreteness, in the remainder of the paper we will only discuss the case of Dirac dark matter and χχ → χχ scattering. All of our results can be straightforwardly generalized to the other cases.
We denote the dark matter field as χ. The consistency of a Dirac effective theories requires a good approximate global symmetry, In the absence of such a symmetry, we could write Majorana mass terms which split the Dirac fermion into two Majorana mass eigenstates. We now list various cases that we will individually consider in this paper. The Lagrangian for dark matter coupling with a real scalar φ is, If instead we wish to couple to a pseudoscalar φ, we obtain the following, Note that we have assumed a CP symmetry to a good approximation, so that φ has well-defined CP quantum numbers. We can also couple the dark matter to a U(1) gauge field φ µ . If the charge assignments are vectorlike (i.e. we simply gauge the Dirac U(1) symmetry), then we end up with the vector interaction of the gauge field with the dark matter. In order to generate a mass for the U(1) gauge field, we can couple it to a charged scalar that gets a vev. Imposing a separate global U(1) Dirac symmetry under which the scalar is a singlet forbids a coupling of this scalar to the dark matter.
To obtain a purely axial vector interaction of Dirac fermion dark matter is a little intricate. For clarity, we start with two left-handed Weyl fermions, χ 1 and χ 2 , with the same charge under the U(1) gauge group. In this case, the Dirac mass term is not invariant under a U(1) gauge transformation, and must arise from the U(1) breaking. This also implies that the coupling of an axial vector mediator to a massive fermion must vanish in the limit that the vector boson mass goes to zero. This distinguishes the axial vector case from the other cases, in which we are free to take the couplings to be order-one numbers. Explicitly, we write an abelian Higgs model with a Higgs boson of charge 2: 1 Again, an additional global U(1) Dirac symmetry needs to be imposed to ensure that terms like Hψ 1 ψ 1 are absent. After SSB, we can expand the Higgs around its expectation value: . Keeping the relevant terms, we find in unitary gauge, The masses and couplings of the massive vector, radial Higgs mode, and fermion are given in terms of the fundamental parameters as, In particular, this scenario predicts that the coupling of the axial vector mediator to the dark matter behaves as so light axial-vector mediators are necessarily weakly coupled. From an effective field theory point of view, all these interactions are equally well motivated to analyze. Furthermore, if we keep an eye towards UV completions, some of these scenarios arise more naturally than others. Light scalars that are not pseudo-Nambu-Goldstone bosons are unnatural. Light vectors are also unnatural, if their mass comes from Higgsing, because the Higgs boson mass must also be protected. On the other hand, pseudoscalars are particularly well motivated because we know examples of underlying dynamics which can generate pseudoscalar couplings with an associated light boson, such as pions in QCD. So it is much easier to embed a light pseudoscalar into a UV complete theory.

General Procedure
In this section, we outline the general procedure of obtaining the scattering cross sections, starting from a weakly-coupled Lagrangian. For light mediators, there can be substantial effects from multiple exchanges of the mediator in non-relativistic scattering processes, so that the tree-level approximation does not reflect the true answer. A convenient way to resum these contributions in this case is to map the problem on to the equivalent quantum mechanical scattering problem. We can then solve the QM problem non-perturbatively and extract the scattering matrix elements. These amplitudes include the putative Sommerfeld enhancement effects.
Before we outline the procedure, it is worth highlighting the validity of the procedure carried out below. At small enough values of the coupling, the scattering amplitude at some fixed velocity is wellapproximated by the the tree-level amplitude in the QFT. If the relativistic corrections are small, the same amplitude is also well described by the Born approximation in a quantum mechanical system. The scattering potential is calculated by matching a QFT amplitude with the corresponding QM amplitude in the Born approximation. However, a subtlety that we will encounter involves potentials where the higher-order Born terms are "divergent". This is the case for singular potentials that grow faster than 1/r 2 as r → 0. In such a case, to calculate the higher-order Born terms we would need to regulate and renormalize the quantum mechanical problem, and consistently match with the QFT at loop level. Indeed, this is true for the potentials we consider, but these potentials do have a well-defined first order Born limit as r → 0 which will be sufficient for us to avoid the issue of divergences.
Our problem neatly factorizes into two pieces. The short-distance r → 0 piece is the part where the potential divergence shows up, but from the QFT point of view this corresponds to high energy scattering, where Sommerfeld enhancement should not play a role and the amplitude should be well approximated by the tree-level diagram, or equivalently the first Born approximation in the QM picture. The potential Sommerfeld enhancement from multiple exchanges of the mediator appears in the larger r region, where the QM potential is well-behaved.
Thus, we use the following procedure, described in further detail below. We match the QFT tree-level amplitude with the QM amplitude in the first Born approximation to calculate the matching condition close to the origin r = a. Here a 1/m χ serves as a UV cutoff. Using this as the boundary condition, we solve the Schrödinger equation numerically, and thereby derive the scattering amplitude. We illustrate this procedure schematically in Figure 1. We then proceed extract the S-matrix from r u(r) Deep Figure 1: A schematic of our matching procedure. In the deep UV, we expect the Born approximation to hold, which we use to set our boundary conditions at r = a. Beyond this is the IR of our theory where Sommerfeld enhancement might be important and we need to solve the Schrödinger equation with the appropriate potential.
the QM solution, highlighting the coupled channel case. We will then be in a position to critically evaluate which interactions give rise to a Sommerfeld enhancement.

Computing the Tree Level Potential
To compute the tree level potential, we first compute the tree-level perturbative QFT amplitude for the process we are interested in. To illustrate our procedure, we consider χχ → χχ. At tree level, this has contributions from an s-and a t-channel Feynman diagram. Next, we take the non-relativistic limit of this amplitude and keep the leading terms. The scattering amplitude in the Born approximation in non-relativistic quantum mechanics is given by Comparing this expression to the non-relativistic limit of our QFT amplitude gives us where the extra factors of 2E p come from the difference in the conventional normalizations of relativistic and non-relativistic single particle states. Finally, to compute the real space potential V ( r), we have to Fourier transform V ( q) with respect to q.
We obtain the following potentials The terms in these potentials arising from t-channel contributions have been previously computed (e.g., [29,34,40,41]). The s-channel contributions provide contact terms that we have written in the form of a nonrelativistic potential via Fierz rearrangement (see, e.g., [32]). As detailed in Appendix B, the spin matrices for antifermions come with an additional minus sign.
In the pseudoscalar potential, we observe two delta-function terms. The first term comes from the s-channel, but the second term comes from the t-channel, which can be seen in the following manner:

The Schrödinger Equation
We solve the Schrödinger equation in the partial wave expansion using the potential derived above. The rotational invariance of the Hamiltonian implies that j is a conserved quantum number. For the 2 → 2 scattering process involving spin-1 2 particles that we consider, the total spin is s = 0, 1. Consequently, it will be convenient to work in a basis of states |r, j, σ, , s , where and s are the total orbital and spin angular momenta respectively, and σ is the z-component of the total angular momentum. Some of the terms in the potentials we consider mix terms with different values, so for a given j, we get a 4 × 4 block diagonal Hamiltonian. The wavefunction ψ( r) can be separated, The potential preserves j, σ, so without loss of generality we can set σ = 0. The Schrödinger equation for a fixed value of j becomes The matrix elements of the potentialV (r) are given in Appendix A. The subscripts ( , s) take on 4 possible values {(j, 0), (j − 1, 1), (j, 1), (j + 1, 1)}. We have suppressed the j and σ quantum numbers for notational clarity. Conveniently, the four ( , s) states separate into two states that mix with each other (|j ± 1, 1 ) and two that evolve independently (|j, 0 and |j, 1 ).
Assuming that rV (r) → 0 as r → ∞, the asymptotic wavefunction should be a solution to the free particle equation. The solutions of the radial free particle equation, denoted s and c , are given in terms of the spherical Bessel functions. The K-matrix is the generalization of the more familiar partial wave phase-shift tan δ l . By solving the coupled differential equations numerically, we can extract the K-matrix, and then calculate the scattering cross section. We describe these steps in detail next.

Setting Boundary Conditions
The boundary conditions for the Schrödinger equation above are set using the Born approximation in the region r < a = 1/m χ as outlined in Section 3. This region of small radii probe the UV of our effective quantum mechanical description and we expect this to match onto the corresponding QFT. Therefore, we expect the first Born approximation to reproduce the tree level perturbative QFT approximation. We show this matching in Figure 2. The K-matrix in the Born approximation is simple to calculate, where we have restricted the integral to r < a. This gives us our boundary conditions for the numerical solutions (3.14) which we can use to evaluate the wavefunction and its derivative at r = a = m −1 χ . In our numerical calculations below, we have checked that varying the matching radius has little effect on our results, provided that we choose a of order m −1 χ . In this whole discussion so far, we have not specified the structure of the potential. In particular, potentials in quantum mechanics can diverge at the origin, and the Born limit may not be welldefined. A potential is considered non-singular if r 2 V (r) → 0 as r → 0. Since s (kr) ∼ r l+1 for small r, for such potentials clearly the integral above is convergent. The pseudoscalar potential does not fit this criterion due to the r −3 piece and naively looks singular when = = 0. However, as we see from the results in Appendix A, these dangerous terms vanish under the action of the operator O T ≡ 3( S 1 ·r)( S 2 ·r) − S 1 · S 2 appearing in the numerator.
Since we consider all potential terms that are generated from tree-level exchanges, it is interesting to note that the QFT seems to produce highly non-generic potentials which may seem singular, but possess operator structures that remove these divergent pieces at the leading Born approximation.

Extracting the S-Matrix from Numerical Solutions to the Schrödinger Equation
The K-matrix can be extracted from the set of asymptotic solutions u The S-matrix follows, We outline the steps involved in obtaining the cross section from the S-matrix in this basis following the notation in Ref. [42]. The differential cross section is given in terms of the scattering amplitude in the familiar way, dσ k , σ 1 , σ 2 →k , σ 1 , σ 2 dΩ = f k , σ 1 , σ 2 →k , σ 1 , σ 2 2 .

(3.19)
We can write f in terms of the (jσ s) basis S-matrix using Wigner-Eckart theorem, f k , σ 1 , σ 2 →k , σ 1 , σ 2 = − 2πi k jσl s p lsp C 1 2 1 2 (s, p; σ 1 σ 2 )C ls (j, σ; m, p) where the C's are Clebsch-Gordan coefficients. We have specialized to the case of elastic 2 → 2 scattering between spin-1 2 particles. The spin-averaged cross section in this basis is, (3.21) In SIDM phenomenology, there are other related quantities of interest, such as the momentum transfer cross section or the viscosity cross section. Their utility arises from their well-behaved soft and forward limits, and they have a direct bearing on the evolution of dark matter phase space in halos. They are defined in terms of the differential cross section. The transfer cross section is defined as [17,27] σ T = dσ dΩ (1 − cos θ)dΩ (3.22) while the viscosity cross section is defined as

Results
In this section, we show the results from the procedure laid out in Section 3.3. The quantity of merit is the viscosity cross section σ V . The angular weighting regulates both forward and backward scattering, which is important since singularities in the forward and backward scattering limit, although physical, do not change the DM velocity distribution and hence have no observable effect. The physical quantity we can extract from measurements is σ V v . The velocity averaging assumes a Maxwell-Boltzmann distribution truncated at an escape velocity v esc = √ 2v rms , as for a virialized halo. The results are shown in Figure 3 for the various interactions we considered. The parameters are chosen such that the cross sections are at approximately the correct order of magnitude over the velocity range of interest. A more dedicated exploration of the viable parameter space fitting the self-interaction cross section measurements from astrophysical data in [19] should be performed, but is beyond the scope of this work.
As we discuss in detail in Section 4.1, our numerical results for the pseudoscalar mediator show no Sommerfeld enhancement. So, we can compute σ V v directly from the QFT. We notice a kink occurring at v ∼ O(m φ /m χ ). The kink exhibits a characteristic factor of 3 increase in the cross section. The cross section plateaus before and after the kink. These are robust predictions for the behavior of the pseudoscalar interaction. The features we highlighted above can be seen more clearly in Figure 2, even though it contains only the tree-level QFT results, since these are a good approximation to the answer in the case of pseudoscalar and axial vector interactions.
As discussed in Section 2, the axial vector interaction is non-trivial. It comes along with a radial Higgs mode that also couples to the fermions. This is a scalar particle coupled to our fermionic dark MeV, λ q = 4π, λ = 10 −3 and y = 1. We assume the dark matter follows a Maxwell-Boltzmann distribution and use a hard cutoff at the escape velocity.
matter, which in turn will induce a Yukawa potential between them. As we saw already, Yukawa potentials generate a Sommerfeld enhancement. On the other hand, if we want to study just the axial vector type interaction, we need to decouple the Higgs mode by making it heavier. One consistent way of achieving this is by making the square root of the Higgs quartic larger than the Higgs coupling to the dark matter ( 2λ q > y). Furthermore, to make the axial vector mediator lighter than the dark matter, we need the gauge coupling to be smaller than the Higgs-dark matter coupling (y > √ 8λ). By doing so, we can decouple the Higgs and study a pure axial vector theory. Formally, by integrating out the Higgs, we generate a four-Fermi interaction which manifests as a delta function in the potential since it was generated by a contact interaction without any q dependence.
Having done so, we can now isolate the potential generated by a purely axial vector type interaction. By inspecting Equation 3.6, we see a term proportional to V pseudoscalar . Again, this term does not generate a Sommerfeld enhancement. In addition, there is a Yukawa term in the potential. This term does not generate a Sommerfeld enhancement because in the light mediator regime, the coupling, which is proportional to the mediator mass, is also small. On the other hand, if we want a large coupling, then the mediator mass increases proportionally and Sommerfeld enhancement turns off. The structure of the underlying UV completion restricts the parameter choices we can make and therefore, we can compute σ V v using the QFT when we work in the decoupling limit.

Sommerfeld Enhancement
In this section we study the existence of Sommerfeld enhancement for various types of interactions. We begin with a discussion of our results for the Yukawa potential, which has previously been studied extensively [17,18,26,27]. We then present novel results for the pseudoscalar and axial vector case, which do not show Sommerfeld enhancement. We discuss a diagrammatic argument that reinforces this conclusion, and comment on disagreement with previous work that has found Sommerfeld enhancement in the pseudoscalar case.
In Figure 4, we present our results for the Yukawa potential generated via a scalar interaction. The solid orange curve is the tree-level QFT cross section. Using the procedure outlined in Section 3.3, we set the boundary conditions and solve the Schrödinger equation. The numerical cross section we obtain is plotted as the solid blue curve. It is worth noting here that there is a considerable enhancement in the cross section in the non-relativistic regime. At higher velocities, the QFT tree-level cross section becomes a better approximation for the cross section. As has been discussed previously in the literature, we expect Sommerfeld enhancement to be a significant effect in the non-relativistic limit for light mediators and our results are in good agreement with this expectation.
Previously, in Figure 2 we showed the agreement between the first Born approximation and the QFT expectation for the cross section for a pseudoscalar mediator. Having set the boundary conditions, we compute the numerical cross section, shown in blue, and compare it to the tree-level QFT cross section in dashed orange in Figure 5. At large velocities, the two curves deviate by O(1%). This discrepancy is saturated by the non-relativistic corrections which scale as β 2 and we see that there is no additional enhancement. While we might have expected our intuition from the Yukawa potential to apply here, we do not find any significant Sommerfeld enhancement for any values of m φ .
The absence of enhancements in the small m φ regime can be understood intuitively. The pseudoscalar potential, in Equation 3.4, has two Yukawa like terms which in principle can generate an enhancement, but both of these terms are suppressed by m 2 φ m 2 χ . Therefore, in the limit of small m φ , the suppression of these terms shuts off possible Sommerfeld enhancement. Furthermore, the more divergent terms, such as the r −3 term, behave effectively like a short range potential and hence don't generate enhancements. In the large m φ limit, we don't expect Sommerfeld enhancement in either the Yukawa or pseudoscalar case as the exponential suppression takes over and we have short range potentials.
Feynman diagrammatic arguments lend further credence to this result. We discuss the general idea and results of this argument for the scalar and pseudoscalar case here and relegate the details to Appendix C. Sommerfeld enhancement occurs in a region of phase space where diagrams with two mediator exchanges are comparable or parametrically larger than a tree level diagram with a single mediator exchange. To understand the scaling behavior in this regime, we can analytically compute the resulting box diagram for the cases of scalar and pseudoscalar mediators.
We study the ratio of M 1-loop to M tree . For the pseudoscalar case, the leading behavior goes like g 2 16π 2 log ξ, where we have defined ξ = m 2 χ /m 2 φ . On the other hand, for the scalar case the leading behavior goes like g 2 4π √ ξ. The dependence on the coupling g is as expected. As always for a perturbative calculation, if we make the coupling arbitrarily large, then loops will be important. What is relevant for the existence of nonrelativistic enhancement is the scaling with ξ, which is different for the scalar and pseudoscalar cases. At large ξ, the ratio in the scalar case diverges, whereas in the pseudoscalar case, it behaves like a log-enhanced loop factor. This indicates that when we make our mediator light, the box diagram becomes important and dominates over the tree level diagram in the scalar, or more generally the Yukawa case. This is Sommerfeld enhancement. It requires that we resum an infinite family of diagrams, which is accomplished by solving the nonrelativistic Schrödinger equation. On the other hand, in the pseudoscalar case, the box diagram is not enhanced by nonrelativistic effects as the mediator becomes light, which corroborates the numerical results and heuristic arguments we gave above.
Our conclusions deviate from certain claims in the literature. Our analysis of the EFT and its operators largely agrees with that of [34], though they included only t-channel contributions from light mediators and neglected contact terms. Both [34] and [35] correctly explain how pseudoscalar exchange can couple modes with l = j ± 1. However, [34] then argued that scattering should be dominated by a singular 1/r 3 potential and proceeded to analyze a single-channel equation with such a potential as an idealization of the pseudoscalar exchange potential. We do not believe that this toy single-channel problem has similar physics to pseudoscalar exchange. In particular, our matching procedure works well for the pseudoscalar exchange problem but is not even well-defined for the l = 0 mode in a 1/r 3 potential, because the integral computing the first Born approximation diverges at small r. This reflects the fact that a 1/r 3 potential is singular whereas, as we have argued, a [3(S 1 ·r)(S 2 ·r) − S 1 · S 2 ]/r 3 potential is not (at least at leading order in the Born approximation). The authors of [35] studied the coupled-channel problem for j = 1 modes and concluded that there can be significant Sommerfeld enhancement. Their approach is similar to ours: they have introduced a short-distance square well regulator, and then solve the full coupled-channel problem at longer distances. However, they have not carried out a detailed matching to perturbative QFT, and as a result they treat the depth of the square well potential and the coupling strength as free parameters to vary. The Sommerfeld enhancement that they observe arises in regions of nonperturbatively large short-distance scattering. We believe that there is no conflict between their numerical results and our claim that Sommerfeld enhancement does not arise within the parameter space one can obtain through perturbative matching to a weakly-coupled quantum field theory. Velocity β σ (GeV) -2 Figure 5: Cross section as a function of velocity for dark matter coupled via a pseudoscalar mediator. The numerical cross section (solid blue) computed with the procedure outlined in Section 3 is compared with the tree-level QFT cross section (dashed orange). We set λ = 1 10 , m χ = 1 GeV and m φ = 10 −1 GeV. At low velocities we do not see any Sommerfeld enhancement. We begin to see deviations at larger velocities, as expected since the tree-level QFT answer is a fully relativistic calculation but the numerical cross section is determined from a non-relativistic potential.

Conclusions
In this work, we have studied the velocity dependence of interactions between spin-1/2 dark matter particles mediated by a light boson. In particular, we studied the scenario where the boson is a scalar, vector, pseudoscalar or axial vector. We derived the associated potentials including both s-and tchannel contributions to the scattering process. We outlined a new procedure for setting the boundary conditions where we match a tree-level perturbative QFT estimate of short-distance scattering to the Born approximation to short-distance scattering in the effective non-relativistic quantum mechanical theory. Numerically solving the Schrödinger equation then allows us to capture the effect of Sommerfeld enhancement. While we have only considered simplified models in this work, our procedure generalizes straightforwardly to more complicated models as well.
We presented numerical results for the scalar and pseudoscalar case. The scalar mediator generates significant enhancement for low velocities and light mediators, an effect that has been studied extensively in the literature previously. Our numerical results for the pseudoscalar mediator show an excellent match to the tree-level perturbative QFT approximation. This lack of Sommerfeld enhancement is further supported by a Feynman diagrammatic argument. We also argued that Sommerfeld enhancement is absent for an axial vector mediator. In this scenario, the mediator mass and the gauge coupling are tied together such that as the mass is dialed down, the coupling gets correspondingly weaker. Our results suggest that, if the shape of the cross section discussed in [19] persists with more data and a better understanding of baryonic effects, then pseudoscalar and axial mediators will not fit the data as well as scalar and vector mediators. (These light-mediator models are not the only options, however; see, e.g., [43]).
The matching procedure that we have described can be applied beyond the four simple models we have studied. For example, the axial vector model also in general has scattering mediated by the Higgs boson that provides a mass to the axial vector field. One could match to a theory that includes both the Higgs and axial-vector contributions. In general, our matching procedure will be useful in cases with both long-range Sommerfeld enhanced scattering and short-distance contributions to the amplitude. Once the velocity-dependence of the cross section in a given model is known, it can be fed into simulations or other studies of structure formation, for instance in the ETHOS framework [44][45][46].
Our results also raise a more abstract question: to what extent are singular potentials in quantum mechanics relevant when matching to an underlying perturbative QFT? There is a large literature on singular quantum-mechanical potentials like 1/r 3 . We have observed that when matching to tree-level QFT, such terms are accompanied by spin-dependent factors that eliminate the dangerous terms in the leading-order Born approximation. In the future, it would be interesting to better understand the general properties of quantum mechanical models arising from weakly coupled QFTs.

B Fermion/Antifermion Spin Matrices and Minus Signs
Following the conventions of [47], we have the following definitions: This means that for a particle, we have ξ s † σ z ξ s = ±ξ s † ξ s = ±δ ss (B .3) and for an antiparticle, we have ξ s † σ z ξ s = ∓ξ s † ξ s = ∓δ ss (B.4) More generically, let's identify η = ξ * χ . is the antisymmetric tensor and σ satisfies the relation: Now this is a scalar quantity so we are free to transpose it and this gives us −ξ † χ σξχ, where this σ gets identified with the spin matrix S. From this, we see that the minus sign naturally arises when looking at η † ση for any generic state, and not just the z eigenstates. Figure 6: The box diagram corresponds to the first diagram in the infinite set of ladder diagrams being resummed by our procedure. Sommerfeld enhancement arises when this diagram gives a contribution that is comparable to or larger than the tree level contribution to the scattering process.

C Feynman Diagrammatic Argument for Sommerfeld Enhancement
There are a few lines of evidence supporting our claim that Sommerfeld enhancement is absent in the pseudoscalar case. We can show this analytically by computing the box diagram in Figure 6. The general amplitude for this diagram is given by where the numerator N is The Γ matrices represent the matrix structure arising from the vertices. We will focus on two cases: the Yukawa interaction, where Γ is the identity matrix, and the pseudoscalar case, where Γ = γ 5 . On the equations of motion, for the pseudoscalar case, the numerator simplifies to u(p 3 )/ l u(p 1 )v(p 2 )/ l v(p 4 ).

(C.3)
We introduce Feynman parameters and perform the integral over the loop momentum. Since we are interested in the non-relativistic regime, we take the v → 0 limit of the amplitude. This allows us to perform the integration over two of the Feynman parameters and we are left with the following expression.

(C.6)
A series expansion around large ξ yields The box diagram then gives a contribution to the matrix element scaling as This should be compared with the tree-level amplitude, which scales as M tree ∼ g 2 . Note that although the t-channel contribution to the tree-level amplitude is momentum suppressed in the nonrelativistic limit, the s-channel contribution is not. Hence M 1-loop /M tree is on the order of a naive, log-enhanced loop factor, and nonrelativistic effects do not enhance the cross section predicted by perturbative QFT. By contrast, for the scalar case, we expect to find a Sommerfeld enhancement. The numerator simplifies to 2 N = u(p 3 )(/ l + 2m χ )u(p 1 )v(p 2 )(/ l + 2m χ )v(p 4 ). (C.9) Computing the integral over Feynman parameters analytically for the scalar case yields − 1 2ξ 2 √ −1 + 4ξ (−2 + 6ξ − 32ξ 3 ) arctan 1 √ −1 + 4ξ + (2 − 6ξ + 32ξ 3 ) arctan 1 − 2ξ √ −1 + 4ξ + −1 + 4ξ(−2ξ − log ξ + ξ log ξ) . (C.10) For large ξ, this behaves like 4π √ ξ. This means that the Sommerfeld enhancement can be important for low mediator masses, when This is consistent with the standard claim about the regime of nonrelativistic enhancement for a Yukawa potential (e.g., [27]). One can understand the origin of the enhancement as follows: the numerator of the integral scales as ξ and the denominator as ξ 2 , so in most of the integration region one expects a suppressed contribution at large ξ. However, when (w − x) 2 ∼ ξ −1 , the denominator takes order-one values and the integrand is of order ξ. This occurs only in a region of size ξ −1/2 within the overall integration region, and accounts for an integral of size √ ξ. This did not occur in the pseudoscalar case, where one can check that, precisely when the denominator becomes of order one, the numerator is suppressed as well.