Unitarity in the non-relativistic regime and implications for dark matter

Unitarity sets upper limits on partial-wave elastic and inelastic cross-sections, which are often violated by perturbative computations. We discuss the dynamics underlying these limits in the non-relativistic regime, namely long-range interactions, and show how the resummation of the 2-particle-irreducible diagrams arising from squaring inelastic processes unitarizes both elastic and inelastic cross-sections. We provide a simple prescription to obtain the unitarized cross-sections from those that do not include resummation of the squared inelastic processes. Our results are model-independent, apply to all partial waves, and affect elastic and inelastic cross-sections, with extensive implications for new physics scenarios, such as dark-matter freeze-out, indirect detection and self-interactions.


Introduction
Unitarity has long stood as a fundamental principle in particle physics, crucially invoked in the past to argue for the existence of then-unknown physics, including the electroweak interactions and the Higgs mechanism [1][2][3].More recently, it has been extensively employed to circumscribe the range of validity of effective theories aiming to parametrize the still-unknown physics beyond the Standard Model.In this respect, the implications of unitarity have been widely considered in the highenergy limit, relevant for new physics searches at colliders.
In this letter, we consider unitarity in the antipodal limit, the non-relativistic regime [4].In the forefront of our motivation lies dark matter (DM), which makes up 85% of the non-relativistic matter in our Universe.Besides the DM phenomenology in today's Universe, many important cosmological events, such as the DM production in the early universe, may also fall within the non-relativistic realm.Our results pertain to the cosmology and astrophysics of DM in a variety of theories at the frontiers of research, but are also more generally applicable to any non-relativistic system.
Unitarity sets upper bounds on partial-wave elastic and inelastic cross-sections that indicate the saturation of the corresponding scattering probability to one.In a theory that is by construction unitary, the ramifications of these bounds are twofold: (a) They may limit unknown physical parameters if phenomenological requirements are imposed that amount to certain processes being sufficiently efficient.(b) They indicate the point where approximations employed to compute scattering amplitudes fail.While this failure may be simply due to the break-down of a perturbative expansion in a parameter that has Email addresses: marcos.flores@phys.ens.fr(Marcos M. Flores), kalliopi.petraki@phys.ens.fr(Kalliopi Petraki) become large, in certain cases it hints towards the onset of new effects.
In the present work, we reframe these points, posing the questions: (a ′ ) What theories give rise to physical processes that can approach or attain the unitarity limits on scattering cross-sections, thus allowing to approach or attain the bounds on physical parameters set by phenomenological requirements?(b ′ ) How can we improve perturbative calculations in weaklycoupled unitary theories to ensure consistency with unitarity?We begin with reviewing the unitarity bounds, before addressing these points and discussing examples in the DM context.

Off-shell partial-wave optical theorem
To minimize technicalities that are tangential to the main goal of this work, we neglect spin in the following.We consider 2to-2 transition amplitudes that may in general be off-shell, an essential ingredient in our analysis.We expand them into partial waves as follows, where a and b denote the initial and final state, s is the first Mandelstam variable, and k a stands for the 3-momentum of each particle in the center-of-momentum (CM) frame in the state a.
Note that the two interacting particles may in general be different.Next, we define the rescaled partial-wave amplitude where we have introduced the symmetry parameter δ a = 0 or 1 if the two particles in the state a are distinguishable or identical where b runs over all on-shell states into which a can scatter, and 2  dΠ b (2π) 4 δ 4 ({k} Equation ( 4) holds for any m ∈ [−ℓ, ℓ], and can also be averaged over m.Curly brackets denote collectively the momenta of multi-particle states, and dΠ b is the phase-space integration measure for the n-particle state b, where f b is the symmetry factor of the state b that ensures the phase space is not multiply counted.For 2-particle states b, f b = 1/2 δ b , and Eq. ( 4) reduces to with |k b | determined by s due to b being on-shell.Equations ( 3) to (6) generalize the optical theorem, and we shall employ them in different ways below.Notably, X ab ℓ (s, |k|, |k|) ⩾ 0, which leads to the important observation that the sign of the imaginary part of the elastic amplitude is not arbitrary, Im[M aa ℓ (s, |k|, |k|)] ⩾ 0. This fact will prove critical in the following.

Unitarity bounds
Next, we focus on on-shell in and out states.We define the cross-section with the symmetry parameter δ and the CM momentum k referring to the incoming state.Then, all cross-sections can be expressed in the convenient form where for 2-particle states, M ab ℓ is given by Eq. ( 2  For on-shell states, Eq. (3) implies This can be recast for M elas ℓ = M aa ℓ as which is the familiar unitarity circle centered at i/2 on the complex plane.We may obtain a more restrictive combined bound on elastic and inelastic cross-sections as follows.The inequality Eq. (10) implies Im(M elas ℓ ) ⩽ |M elas ℓ |.Then, recalling Eq. ( 8), which can be recast as The above encompasses the individual bounds σ inel ℓ ⩽ σ U ℓ /4 and σ elas ℓ ⩽ σ U ℓ .It agrees with [6] except for the symmetry factor 2 δ in σ U ℓ that has been missed in the past.The unitarity circle (10), the combined bound on elastic and inelastic crosssections (11), and the mapping between them are depicted in Fig. 1.

Attaining the unitarity limit
The dependence of the cross-section (7) on the kinematic variables allows us to address question (a ′ ).Let us first motivate this phenomenologically with an example.
Reference [7] argued that the inelastic unitarity bounds imply an upper bound on the mass of thermal-relic DM.The calculation amounted to bounding the so-called canonical DM annihilation cross-section by the inelastic s-wave unitarity limit, , where we set |k| = m DM v rel /2 in the non-relativistic regime, with m DM and v rel being the mass and relative velocity of the annihilating DM particles, and we used the current measurement of the DM density to determine σ can .In this comparison, there is clearly a discrepancy between the velocity scalings of σ can and σ U 0 .Reference [7] claimed that the momentum scaling of σ U ℓ cannot be obtained in the non-relativistic regime (and simply set v freeze−out rel ∼ 0.5 in order to evaluate σ U 0 ), and that s-wave annihilation dominates.This however, is not so, as unitarity implies [4,8] and we shall now discuss.This, in turn, has implications for the upper limit m DM < 116 (82) TeV, deduced according to the considerations of Ref. [7] 1 for (non-)self-conjugate thermal-relic DM, but, more importantly, for the underlying physics of viable thermal-relic DM scenarios at the multi-TeV mass scale [4].
In the relativistic regime, |k| → √ s/2, thus σ U ℓ ≃ 2 δ 16π(2ℓ + 1)/s, leading to the well-known conclusion that cross-sections should decrease at least as fast as σ ∝ 1/s at high energies, in order to be consistent with unitarity up to arbitrary high scales.However, for processes computed within an effective theory where a large scale Λ > m has been integrated out, with m being the mass of the interacting particles, the cross-sections scale as σ ∝ (1/s)(s/Λ 2 ) p , where p > 0 depends on the effective theory.Unitarity then implies that the effective theory breaks down at s ≳ Λ 2 , where the high-energy physics must be resolved and incorporated in the computation.If doing so realizes the scaling σ ∝ 1/s, then approaching or attaining the unitarity limit in a continuum of energies s > Λ 2 amounts essentially to the couplings involved being sufficiently large.
If Λ ≲ m, then the breakdown of the effective theory and the transition described above occur in the quasi-relativistic or non-relativistic regime.As previously, the scaling σ ∝ 1/s can retain consistency with unitarity at increasing energies since s > |k| 2 , thus σ U ℓ ∝ 1/|k| 2 > 1/s.However, now, this scaling does not allow to approach or attain the unitarity limit in a continuum of low momenta, |k| ≲ m.To identify what can realize the unitarity limit in this case, we repeat the previous considerations in the non-relativistic regime.
For concreteness, we consider a toy model, and follow in part the discussion of Refs.[4,9].Similar considerations to what follows can be employed in a variety of theories, including gauge theories.Let χ be a scalar of mass m, coupled to a scalar force mediator φ of mass M via the trilinear coupling L ⊃ −(y/2)m χ 2 φ, with y being a dimensionless coupling.We define α ≡ y 2 /(16π) for convenience.The mediator φ may also couple in a similar fashion to other lighter species.χχ pairs can scatter elastically or annihilate via φ exchanges.If φ is sufficiently heavy, then it can be effectively integrated out, giving rise to the χχ elastic and annihilation cross-sections σ elas ∼ σ ann v rel ∼ α 2 s/M 4 , with v rel being the χχ relative velocity.Setting s ≃ 4m 2 and |k| ≃ mv rel /2, unitarity implies, as previously, that this effective description breaks down at sufficiently high energies -with the energy now more conveniently parametrized by |k| or v rel -in particular at |k| > M 2 /(mα).This can be recast in the more indicative form, M < m √ αv rel ≪ m.The comparison between scales that arises from this consideration hints to the reason for the apparent violation of unitarity: φ mediates a χχ interaction of range λ inter ∼ M −1 that is comparable to or longer than the geometric mean of the de Broglie wavelength of the interacting particles, λ dB ∼ |k| −1 = (mv rel /2) −1 and their Bohr radius, λ B ∼ (mα/2) −1 .Consequently, virtual φ exchanges cause long-range elastic χχ interactions that distort the wavefunction of the incoming χχ pair, and must thus be resummed in order to properly compute both the elastic and inelastic cross-sections.(More details on resummation will be provided in the next section.)In the context of inelastic scattering, this is called the Sommerfeld effect [10,11].As we discuss next, the resummation of an (attractive) long-range interaction reproduces the same energy dependence as that of the unitarity cross-section (7).
The resummation of the t/u-channel φ exchanges generates an attractive Yukawa potential between the χ particles, V(r) = −αe −Mr /r (see e.g.[12]).In the Coulomb limit (which is valid for M < mv rel /2, see e.g.[13]), this results in the cross-sections where ℓ ] with a ζ-dependent phase.The annihilation cross-sections thus scale with momentum in the same way as the unitarity limit,2 which they may approach or attain in a continuum of momenta provided that the coupling is sufficiently large.Conversely, we conclude that the inelastic unitarity limit in the non-relativistic regime can be realized in a continuum of energies only if the incoming particles interact via a longrange force [4,8].Indeed, an interaction of range comparable to or smaller than |k| −1 introduces a new scale -the range of the interaction -that can only spoil the 1/|k| 2 scaling of the cross-sections. 3eturning to the phenomenological example of DM freezeout, we may conclude that thermal-relic DM can be as heavy as unitarity allows only if it possesses attractive long-range interactions [4].Such interactions imply the existence of bound states, and it has been shown that the formation and decay of metastable bound states in the early Universe can deplete DM very efficiently [8], also via higher partial waves [4].The significant contribution of higher partial waves is in contrast to the direct annihilation processes, for which higher ℓ modes are suppressed by v 2ℓ rel or α 2ℓ at tree level and in the Sommerfeld-enhanced regime respectively, as seen above.In the case of bound-state formation, the bound-state wavefunction introduces negative powers of the coupling, compensating for this suppression [15,16].How many and which partial waves are important depends on the model, there is thus no model-independent upper bound on the mass of thermal-relic DM.Nevertheless, considering carefully the unitarity limit has allowed us to understand and predict, in a model-independent fashion, the underlying physics of viable heavy thermal-relic DM models.

Unitarity violation
Our discussion of the dynamics that can approach or attain the (inelastic) unitarity limit in a continuum of energies in the non-relativistic regime, reframes the issue of (apparent) unitarity violation.In the toy model of Section 4.1, the elastic crosssection (12a) always respects the global limit σ elas ℓ ⩽ σ U ℓ .However, the inelastic cross-sections (12b) violate σ inel ℓ ⩽ σ U ℓ /4 at large α.Moreover the cross-sections (12) violate the combined unitarity bound (11).This is made apparent, for example, in the limit σ ann ℓ → σ U ℓ /4, where the elastic cross-section oscillates in the interval [0, σ U ℓ ], in contradiction with the constraint (11) which mandates that σ elas ℓ must approach σ U ℓ /4 as well (cf.Fig. 1).This example represents one of several qualitatively distinct cases of unitarity violation observed in state-of-the-art computations in theories involving long-range interactions: (i) Attractive Coulomb(-like) potentials violate inelastic unitarity at large couplings.In a QED-like theory, this occurs at α ≳ 0.85, considering the dominant s-and p-wave inelastic processes for a fermion-antifermion pair [13].(ii) Radiative transitions between states of different Coulomb(-like) potentials, e.g.bound-state formation with emission of a charged scalar [15,[17][18][19] or a non-Abelian gauge boson [16,20], exhibit nonmonotonic dependence on the incoming momentum, with peaks that can violate unitarity for a finite velocity range, even for fairly small values of the couplings.The problem is exacerbated if the potential in the scattering state is repulsive [15,16].(iii) Potentials generated by the exchange of light but massive mediators exhibit parametric resonances at the thresholds for the existence of bound levels, violating inelastic unitarity at low v rel , even for small couplings [13,21].
These issues highlight a limitation in our understanding of this regime; they also hamper phenomenological studies of, for example, heavy thermal-relic DM or self-interacting DM, both of which feature long-range interactions.This brings us to question (b ′ ) that we now address.
The inequalities (9) cannot be satisfied if the amplitudes are calculated at any finite order in perturbation theory, as the two sides would then be of different order in the couplings of the theory.This suggests that in order to ensure consistency with unitarity, we must resum all interactions of interest.Indeed, returning to the cross-sections (12), the resummation of the elastic χχ interaction (the one-boson exchange), resulted in the elastic cross-section (12a) being consistent with the elastic unitarity bound for all momenta and values of the coupling.On the other hand, the inelastic cross-sections (12b), despite including the resummed one-boson exchanges, also involve a purely inelastic interaction that has not been resummed.However, the optical theorem implies that squaring inelastic interactions generates a contribution to elastic scattering; this must be self-consistently resummed in any computation that invokes these inelastic interactions, in order to ensure consistency with unitarity.
In the following we shall show that the proper resummation of the squared inelastic interactions indeed renders elastic and inelastic cross-sections consistent with the combined unitarity limit (11).We begin with reviewing the resummation process, establishing the appropriate framework.

Resummation
The resummation of 2-particle interactions amounts to a Dyson-Schwinger equation for the 4-point function of the two interacting particles, with the kernel being the sum of the 2particle-irreducible (2PI) diagrams.The latter are defined as those that cannot be cut to give two individual contributions to the same 4-point function.The Dyson-Schwinger equation yields the Bethe-Salpeter equation for the two-particle wavefunction, which reduces to the Schrödinger equation under the instantaneous and non-relativistic approximations (see e.g.[12] for precise definitions and derivations).In momentum space, it reads (13) with E k = k 2 /(2µ) being the kinetic energy of the system in the CM frame.Here, m t and µ stand for the total and reduced masses of the interacting particles.Equation ( 13) can be Fourier transformed to position space using where we left the dependence of V(r, r ′ ) on s implicit.
Two important remarks are in order.First, no on-shell condition should be imposed on the 2PI kernel, K 2PI (s, p, p ′ ), since resummation involves its off-shell iterations of the 2PI diagrams.Note that when on-shell, |p| = |p ′ | are determined by s.In fact, in the instantaneous approximation, the dependence of K 2PI on the energy transfer is neglected (cf.ref. [12]).In position space this yields, in general, non-local potentials V(r, r ′ ).Second, guided by the generalized optical theorem (3), we do not restrict our analysis to the case of K 2PI depending on p − p ′ only, rather than p and p ′ separately.2PI kernels that depend only on p − p ′ lead to central potentials, where all partial waves are related.We treat instead all partial waves independently, and do not specify the momentum dependence of K 2PI .
In a weakly coupled system, the potential vanishes at r → ∞ as V(r, r ′ ) ∼ r −(1+ϵ) with ϵ ⩾ 0.Then, scattering is described by wavefunctions that asymptote to an incoming plane wave plus an outgoing spherical wave at spatial infinity; in position space where k = |k|, and f k (Ω r ) encodes the scattering dynamics.(For V(r, r ′ ) ∝ 1/r, f k has a mild r dependence at r → ∞, which does not spoil any of our conclusions.) As is customary, ψ k and f k may be expanded into partial waves and re-expressed in terms of phase-shifts, where the phase shifts ∆ ℓ (k) are in general complex.The partial-wave cross-sections are (see e.g.[22]) The total cross-section is obtained from the standard form of the optical theorem, and Eq.(20c) includes all inelastic channels.The factors 2 δ incorporated in σ U ℓ arise from the (anti)symmetrisation of the wavefunction in the case of identical particles that doubles the strength of the ℓ modes that survive; for all other ℓ modes, σ ℓ vanish.
To ensure that the combined unitarity bound (11) is respected by the above solution, we must still show that Im∆ ℓ ⩾ 0. This is related to the imaginary part of K 2PI that we discuss next.

Imaginary potential from generalized optical theorem
The generalized optical theorem (3) suggests that inelastic processes generate an imaginary K 2PI component.Resumming this contribution renders inelastic cross-sections consistent with unitarity, as we show next.
An imaginary potential due to inelastic processes has been previously considered in Ref. [23], for the purpose of unitarizing s-wave annihilation.However, the connection with the generalized optical theorem was not identified, resulting in ambiguity on the sign of the potential, which is an important issue.Moreover, a central potential was adopted, which does not allow to treat different partial waves independently [24], and the chosen potential, a δ-function in position space, amounted to a specific choice for the underlying inelastic interaction.(See also [25] for a study of similar scope, and [26,27] for works that do not resum inelastic processes but propose related ansätze to improve perturbative amplitudes.)Our treatment relies on resummation, handles all partial waves independently, and does not make model assumptions for the momentum dependence of the resummed process.The sole constraint is the generalized optical theorem (3), as it should be.
To obtain Im(K 2PI ), we must remove from the right-hand side of Eq. (3) all elastic contributions (which yield 2-particlereducible diagrams), and consider only inelastic contributions with any 2PI factors amputated.We denote the latter as A inel, j .The one-boson exchange diagram is a purely elastic interaction that generates the leading order contribution to Re(K 2PI ).Bottom: 2PI diagrams arising from squaring the 2-to-2 inelastic processes, χχ → φφ and χχ → B(χχ) + φ, where B(χχ) is bound state, described at leading order by an infinite ladder of one-boson exchange diagrams.These diagrams generate subleading contributions to Re(K 2PI ) from off-shell intermediate states, and leading contributions to Im(K 2PI ) from on-shell intermediate states.The latter are described by Eq. ( 21) (see text for explanations).
Note that A inel, j is not the full amplitude for the inelastic channel j; the latter requires resummation of the 2PI diagrams of the incoming particles (cf.e.g.[12]).For simplicity, we will neglect contributions from 2 → n ⩾ 3 scattering.
We partial-wave analyze K 2PI and A inel, j according to Eq. ( 1), and define their corresponding rescaled versions, that we shall call K 2PI ℓ and a j ℓ , according to Eq. ( 2).Then, the generalized optical theorem, Eqs. ( 3) and ( 6), imply where we do not specify the final-state momenta in the inelastic factors a j ℓ since they are on-shell and determined by s.Equation ( 21) is a key element in our analysis.We depict it pictorially in Fig. 2. Examples of contributions to K 2PI from purely elastic and squared inelastic diagrams are contrasted in Fig. 3.
To proceed, we transcribe Eq. ( 21) in position space, according to the Fourier transform (15).We analyze V(r, r ′ ) in partial waves, and obtain with where we set √ s ≃ m t in the prefactor, consistently with the non-relativistic approximation leading to Eq. ( 15) [12].
We shall employ Eq. ( 23) in Section 7, solve Schrödinger's equation and obtain the corrected wavefunction that, as we show, renders elastic and inelastic cross-sections consistent with unitarity.Before doing so, we offer here a simpler proof of this unitarization procedure, valid at leading order.
We consider the current j k (r) ≡ Im[ψ * k (r)∇ψ k (r)]/µ.By virtue of Stokes' theorem and the continuity equation, Using the first Born approximation for the right-hand side, and expanding correspondingly the left-hand side in Im∆ ℓ (cf.Eq. ( 20c)), keeping lowest order terms, we find which implies Im∆ ℓ ⩾ 0 as per Eq. ( 21).(This remains so if we include 2 The above analysis, distilled in Eqs. ( 20), ( 21) and ( 26), shows that the resummation of the imaginary part of the 2PI diagrams -which arises from squaring inelastic processes, an essential point -ensures that the elastic and inelastic crosssections respect the combined unitarity bound (11).While this proof is at leading order in the couplings that give rise to inelastic scattering, and concerns the sum of all inelastic processes, we show in the following how the resummation of Eq. ( 21) ensures that unitarity is respected at all orders, and obtain unitarized solutions for the exclusive inelastic cross-sections.

Unitarization
Next, we compute how the consistent resummation of the inelastic contributions to elastic scattering, and in particular of Eq. ( 21), unitarizes the partial-wave elastic and exclusive inelastic processes.
To begin, we note that the elastic scattering cross-section will be computed according to Eq. (20b), using the phase shifts that take into account the imaginary part of the potential given by Eq. ( 23).The cross-section for an inelastic channel j is , where M inel, j ℓ stands for the corresponding amplitude rescaled according to Eq. ( 2), and incorporates both the hard-scattering inelastic amplitude, a j ℓ , and the effect of the non-relativistic potential on the incoming particles, We must now determine the wavefunctions ψ |k|,ℓ .

Schrödinger's equation with complex potential
We shall consider a complex potential whose real part we take to be central for simplicity, Re[V(r, r ′ )] = V(r)δ 3 (r − r ′ ), while its imaginary part is given by Eq. ( 23).Then, the partialwave Schrödinger equation for u |k|,ℓ (r) ≡ kr ψ |k|,ℓ (r) reads where S ℓ (r) is the differential operator that includes the real part of the potential only, and E k = k 2 /(2µ), as before.We define two versions of the amplitudes and cross-sections: the unregulated ones, for which we consider the real but neglect the imaginary part of the potential, and the regulated ones, for which we take into account the entire complex potential.We denote them by the subscripts 'unreg' and 'reg' respectively, and aim to express the latter in terms of the former.
To do so, we consider the solution of Schrödinger's equation in the former case and the Green's function, where F k,ℓ and G k,ℓ are the regular and irregular families of solutions that feature the following asymptotic behaviors at r → ∞, e ikr e i2θ ℓ (k) + e −i(kr−ℓπ) .(33b) with θ ℓ (k) ∈ R. For some well-known real potentials, such as the Coulomb potential, F k,ℓ and G k,ℓ are known functions [28,29].The functions F k,ℓ (r) describe states of an incident plane wave plus an outgoing spherical wave, and are the solutions of interest.G k,ℓ (r) will be useful for computational purposes.We shall require that G k,ℓ (r, r ′ ) acts as an outgoing spherical wave at r → ∞.
We may now write an implicit solution of Eq. ( 28) for u k,ℓ (r) for the full potential, where we used to express the r-independent factor of the second term in terms of the rescaled regulated inelastic amplitudes, which depend on the wavefunction u k,ℓ .The unregulated amplitudes will become useful shortly.
Inserting Eq. ( 34) into (35a), and introducing the matrix we obtain a matrix equation between regulated and unregulated inelastic amplitudes.Upon inversion, it gives Equation ( 37) is the key result of this work: it expresses the regulated inelastic amplitudes in terms of the unregulated ones.
Similarly, using Eq.(37), the solution (34) can be re-expressed in terms of input parameters only: the imaginary potential and unregulated wavefunctions.

Green's function
To further simplify the results (37) and (38), we focus on the Green's function G k,ℓ (r, r ′ ).We seek an expression for G k,ℓ (r, r ′ ) that is regular at r → 0 and has the asymptotic behavior of an outgoing spherical wave at r → ∞.Following Refs.[24,[29][30][31][32], we expand the r dependence of G k,ℓ (r, r ′ ) in terms of F q,ℓ (r), which, being eigenfunctions of the hermitian operator S ℓ (r), constitute a complete orthonormal set of states with the desired asymptotic behavior,4 With the above decomposition, Eq. ( 32) becomes We must now invert the above to obtain g k,ℓ (q, r ′ ).For this, we need the orthonormality relation for F q,ℓ (r), Considering this, Eq. (40) yields Then, Eq. (39) becomes In order to evaluate G k,ℓ (r, r ′ ), it will be convenient to extend the range of integration over dq to (−∞, ∞).To analytically continue to q < 0, we consider the asymptotic behaviors (33), and set self-consistently With this, Eq. (43) becomes with the prescription ϵ → 0 + chosen such that G k,ℓ (r, r ′ ) behaves as an outgoing spherical wave at r → ∞, as will be shown below.Since F q,ℓ (r) and G q,ℓ (r) are real up to the same ℓand q-dependent but r-independent phase, it follows that Here, this implies G k,ℓ (r, r ′ ) = G k,ℓ (r ′ , r), as expected.
To evaluate Eq. (45), we consider the two independent eigenfunction families of S ℓ , characterized by the asymptotic behaviors of Eqs.(33).We then define The functions H (±) k,ℓ (r) have the properties We invert Eq. ( 47) to express F k,ℓ(r) in terms of H (±) k,ℓ (r), which allows to separate Eq. ( 45) into two integrals, We evaluate the above via contour integration, using the residue theorem, and choosing the contours based on the asymptotic behavior of H (±) k,ℓ (r).For r < r ′ , the first integral must be evaluated in the upper-half plane while the second in the lower half plane.The case r > r ′ follows a similar logic.Considering also the analytic continuations (44b) and (48c), as well as Eqs.(46), we find, where r < ≡ min{r, r ′ } and r > ≡ max{r, r ′ }.The necessary asymptotic behaviour of G k,ℓ (r, r ′ ) is made apparent by the appearance of H (+) k,ℓ in the above expression.If the spectrum of S ℓ includes bound states, they should be included in the Green's function spectral decomposition of Eq. (39).However, in this case, the scattering state wavefunctions F q,ℓ exhibit poles at the imaginary values of q that yield the bound-state energies.Upon contour integration in Eq. ( 49), the residues of those poles cancel exactly the bound-state contributions from the spectral decomposition (cf.e.g.[33]), thereby leaving the final result, Eq. (50), unaffected.

Regulated cross sections
The results of the previous section can be used to derive simple expressions for the regulated elastic and inelastic cross sections in terms of the corresponding unregulated quantities.To begin, using Eq.(45), we can re-express the normalization matrix, Eq. (36), as where, for q ∈ R, we define Equations ( 44b) and (46a) imply h i ℓ (−q) h j ℓ (−q) = h i ℓ (q) h j ℓ (q).With these definitions, the integrand in Eq. ( 51) can be analytically continued in the complex q plane, as it is an analytic function of q except perhaps for isolated poles (meromorphic).If the functions h i ℓ (q) are holomorphic in the upper half of the complex plane and increase at complex infinity more slowly than |q| 1/2 , then we can evaluate the contour integral to obtain, 5   [N ℓ (k)] i j = δ i j + M inel,i ℓ,unreg (k)M inel, j * ℓ,unreg (k).
From here, it is easy to demonstrate that the amplitudes M inel, j ℓ,unreg (k) form an eigenvector of N ℓ (k), (54) 5 For physical, non-relativistic momenta, 0 < q ≪ µ, the functions h i ℓ (q) are related to the unregulated inelastic amplitudes given by Eq. (35b), h i ℓ (q) = √ qM inel,i ℓ,unreg (q).However, M inel,i ℓ,unreg (q) cannot be analytically continued on the entire real q axis, due to the rescaling by the non-analytic factor √ q relative to the full amplitudes in Eq. ( 2).In Eq. (51), we use h i ℓ (q) instead of M inel,i ℓ,unreg (q) to make analyticity manifest.To ensure analyticity in the complex plane, we also use hi ℓ (q), which depends only on q, instead of using [h i ℓ (q)] * = hi ℓ (q * ), which introduces dependence on q * .Finally, h i ℓ (q) is defined via Eq.(52a) for arbitrarily large |q|, while Eq.(35b) holds only for nonrelativistic momenta (cf.Footnote 4).The large-momentum modes are needed for the spectral analysis of the Green's function, employed in Eq. (51).

Phenomenological implications
The ramifications of the regularization Eqs.(60) (or of the more general Eqs.(37), ( 38) and ( 51)) are far-reaching: • The resummation of the squared inelastic diagrams affects both the inelastic and elastic cross-sections rather significantly.It is thus expected to affect the DM density in various production mechanisms, the DM indirect detection signals, and the DM self-interaction cross sections; the latter are important for the galactic structure (see [34,35] for reviews).
• Each inelastic channel is regulated by the inclusive inelastic cross-section.It follows that even if a certain inelastic channel appears seemingly irrelevant for a particular purpose due to its nature, it may in fact be very important, provided that it is sufficiently strong.This is because such a process regulates the pertinent inelastic processes.
For example, in DM freeze-out, metastable bound-state formation processes that are comparable to or faster than direct annihilation into radiation can significantly renormalize down the strength of the latter (cf.Ref. [18,Fig. 11]), even when they themselves cannot deplete DM efficiently due to their rapid dissociation by the radiation bath.This may suppress the DM depletion and modify predictions.
Similarly, for DM indirect detection, rapid inelastic processes with (nearly) unobservable products can renormalize down inelastic processes that produce observable signals, thereby relaxing constraints.
• For y ℓ,unreg ≪ 1, the relative corrections to the scattering rates are δx ℓ /x ℓ ≃ δy j ℓ /y j ℓ ≃ −2y ℓ,unreg .For DM freeze-out, this implies that the imaginary potential affects the DM density at a level similar to or larger than the experimental uncertainty if y ℓ,unreg ≳ O(10 −2 ), i.e., even when the annihilation cross-section is far below the unitarity limit.
• At the opposite limit, y ℓ,unreg ≫ 1, we find x ℓ,reg ≃ 1 − x ℓ,unreg and y ℓ,reg ≃ 1/y ℓ,unreg , i.e., perhaps counterintuitively, the regulated cross-sections decrease as the unregulated ones increase.(See below for an interpretation.)Note that y ℓ,unreg ≫ 1 can occur even for (very) small couplings, e.g. at resonant points, as well as in bound-state formation processes where the potentials of the incoming scattering and outgoing bound state are different [15][16][17][18].The small couplings ensure the validity of the instantaneous and non-relativistic approximations, and thus of the framework for the computations of this work.
With respect to the DM density, this implies that there could be two branches of solutions (mass-coupling correlation) that yield the observed value.
• If σ inel ℓ,unreg scales with momentum in the same way as σ U ℓ , then this feature is retained by the regularization.If not, the regularization generates a non-trivial velocity dependence.The earlier conclusion thus remains: the unitarity limits can be attained in a continuum of non-relativistic momenta only by long-range interactions.
• The upper bounds implied by unitarity on the mass of frozen-out DM annihilating via a single partial wave ℓ, or all partial waves in the range 0 ⩽ ℓ ⩽ L, are, for selfconjugate and non-self-conjugate DM, m U DM /(140 TeV) solely ℓ 0 ⩽ ℓ ⩽ L non-self-conj.DM where we used the numerical results of [4,8], assuming for definiteness that DM annihilates into plasma of the same temperature as that of the Standard Model, and the number of extra relativistic degrees of freedom is negligible.It is straightforward to adjust these assumptions; this would affect the overall mass scale (here, 140 TeV), but not the scaling with ℓ or L. Several remarks are in order: -The difference between the figures of 140 TeV denoted above, and 82 TeV deduced from [7] (cf.Section 4) is due to the momentum scaling of σ U ℓ that has been here properly taken into account [4,8].
-For self-conjugate DM, the above results incorporate the factor 2 δ = 2 of Eq. (7) that has not been previously appreciated in similar estimations.They also take into account that only ℓ even or odd modes survive, depending on the statistics.
-For non-self-conjugate DM, the above takes into account only particle-antiparticle interactions, as is standard.However, DM may be depleted also via particle-particle inelastic scatterings (see e.g.[15] for a model), which increases the bounds on the DM mass shown in the table by a factor √ 3.
-We reiterate that there is no model-independent bound on the mass of thermal-relic DM due to unitarity, since which and how many partial waves contribute significantly depends on the model.
We close this discussion with an interpretation of Eqs.(60) in hydrodynamic terms.The resummation of inelastic interactions accounts for the reduction in the particle flux due to inelastic scatterings, as suggested by the continuity equation [23].This reduction suppresses both elastic and inelastic scattering, by the factor (1 + y ℓ,unreg ) −2 .
For inelastic scatterings, this effect competes against the fact that the probability for inelastic scattering increases with y ℓ,unreg .For y ℓ,unreg ⩽ 1, the increase of the inelastic scattering probability with y ℓ,unreg dominates, while for y ℓ,unreg > 1, the suppression of the particle flux takes over.
Elastic scattering arises from two contributions.Particles may scatter due to the purely elastic couplings with probability x ℓ,unreg .They evade this scattering with probability 1 − x ℓ,unreg , but may scatter inelastically, with the products of the inelastic processes scattering back into the original species via the inverse reactions; this introduces an additional probability factor y 2 ℓ,unreg .The regenerated flux dominates elastic scattering at sufficiently large y ℓ,unreg (except if x ℓ,unreg = 1), where the regeneration probability ∝ y 2 ℓ,unreg exactly balances out the suppression (1 + y ℓ,unreg ) −2 , resulting in x ℓ,reg ≃ 1 − x ℓ,unreg .

Conclusions
The violation of inelastic unitarity bounds in the nonrelativistic regime by state-of-the art computations hampers a variety of phenomenological investigations in the frontiers of DM research.Resolving this issue is important for improving our theoretical understanding, as well as interpreting and guiding experimental searches.We have derived a simple analytical regularization scheme, given by Eqs.(60), whose input are the inelastic cross-sections as affected only by the real part of the potential, and output are the unitarized elastic and inelastic cross-sections.The scheme applies to any partial wave, and is model-independent as it makes no assumptions about the momentum dependence of the unregulated inelastic amplitudes, except for analyticity and convergence requirements.(The more general Eqs.(37) and (38) along with Eq. (51) can be employed when these requirements do not hold.)These results can be, and must be, employed in investigations of new physics scenarios in the non-relativistic regime, with wide-ranging implications.

Figure 1 :
Figure 1: Upper: The unitarity circle bounding the on-shell partial-wave elastic matrix elements M elas ℓ .The border of the circle corresponds to vanishing inelasticity.Lower: Combined bound on the partial-wave elastic and inelastic cross-sections.The symbols identify points between the two plots.The filled red circle corresponds to maximum inelastic scattering.

Figure 2 :
Figure 2: The imaginary part of the 2PI kernel arises from squaring inelastic processes (cf.Eq. (21)): the intermediate particles must be on-shell and different from the incoming / outgoing particles.Here we show only 2-to-2 processes; multiparticle intermediate states require also integration over the allowed intermediate phase space.

Figure 3 :
Figure 3: Examples of contributions to K 2PI in the toy model discussed in Section 4. Top:The one-boson exchange diagram is a purely elastic interaction that generates the leading order contribution to Re(K 2PI ).Bottom: 2PI diagrams arising from squaring the 2-to-2 inelastic processes, χχ → φφ and χχ → B(χχ) + φ, where B(χχ) is bound state, described at leading order by an infinite ladder of one-boson exchange diagrams.These diagrams generate subleading contributions to Re(K 2PI ) from off-shell intermediate states, and leading contributions to Im(K 2PI ) from on-shell intermediate states.The latter are described by Eq. (21) (see text for explanations).