Final state Sommerfeld effect on dark matter relic abundance

If the annihilation products of dark matter (DM) are non-relativistic and if there is some long-range force between them, there can be Sommerfeld effect for the final state particles. We study this effect on DM relic abundance in the thermal freeze-out scenario. As a proof of concept, we consider the case of a DM pair annihilation into a final state pair, assuming that the mutual interactions between the two final state particles give rise to a Coulomb-like potential, and that the masses of the initial and final state particles are similar, so that both the initial and final state particles are non-relativistic. The size of the final state Sommerfeld (FSS) effect depends on the strength of the potential, as well as on the mass ratio of the final and initial state particles. We find that the impact of the FSS effect on DM relic abundance can be significant, and an electroweak sized long-range interaction is large enough to make a correction well beyond the observational accuracy. Another feature of the FSS effect is that it could be suppressed when its time scale is longer than the lifetime of the final state particles. As a corollary, we also study in the DM coannihilation scenario where the initial state Sommerfeld effect between two coannihilators could be reduced due to their instability, which may need to be taken into account for an accurate calculation of the DM relic abundance.


Introduction
The relic abundance of Weakly Interacting Massive Particle (WIMP) dark matter (DM) is usually given by the thermal freeze-out mechanism [1,2]. During freeze-out, the annihilating DM particles are moving with non-relativistic velocities. If there is some long-range interaction between two annihilating particles, such that the two-body wave function is modified from the plane wave, the annihilation cross sections and the DM relic abundance are affected. This non-relativistic quantum mechanics effect -the so-called Sommerfeld effect [3] -on DM relic abundance has been wellstudied. Also, in coannihilation scenarios [4] the long-range interaction strengths between coannihilators can often be larger than the ones between DM particles; therefore, the Sommerfeld effect can be significant in coannihilation scenarios (see e.g. [5][6][7][8][9][10][11][12][13][14][15][16]).
Previous works about the Sommerfeld effect have been focusing on the initial DM or coannihilator pairs. However, the two conditions of the Sommerfeld effect, namely, non-relativistic motion and long-range interaction, may also apply to the annihilation products. In this work, we study the Sommerfeld effect of the final state particles on DM relic abundance.
First, it is possible that the annihilation products can move with non-relativistic velocities. An example is that the annihilation products have masses close to the DM or coannihilators in a 2-to-2 reaction. The annihilation products can be either the Standard Model (SM) particles or new particles in beyond the Standard Model (BSM) theories. The mass degeneracy between the initial and final state particles can be either a coincidence or from a symmetry. For example, if supersymmetry were not broken or only slightly broken, the R-parity odd sparticles would have been nearly mass-degenerate with their corresponding R-parity even partners which can be the annihilation products of the former. In particular, at finite temperature during thermal freeze-out, the initial annihilating particles in the tail of the Maxwell-Boltzmann velocity distribution can have enough energy to annihilate into heavier particles. As one of the "three exceptions" [4], annihilation into forbidden channels is important when these channels dominate the total annihilation cross section and determine the DM relic abundance. This scenario has been further developed recently under the names of "Forbidden DM" [17] and "Impeded DM" [18].
Second, it is natural that there is some long-range interaction between the annihilation products. A typical situation is that the final state is a particle-antiparticle pair, and gauge bosons can be exchanged between them. The gauge groups could be either the ones in the SM or the ones in BSM theories. Exchanges of the Higgs or some other new scalars are also possible.
We will show that for non-relativistically moving annihilation products with long-range interaction strengths comparable to the SM ones, the final state Sommerfeld (FSS) factors are large enough to affect significantly the DM relic abundance, well beyond the percent level accuracy of the observationally determined value [19].
In fact, the FSS factor has been discussed routinely in collider productions of charged particles near the production thresholds 1 , including W + W − , tt, charged leptons, baryons and heavy mesons (see e.g. [20][21][22][23][24][25]), where the long-range interactions are induced by the exchanges of the SM gauge bosons.
There is another feature that needs to be considered in the calculations of the FSS effect. While the DM in the initial state is stable, the final state particles may decay. The FSS effect is ineffective if its time scale is longer than the lifetimes of the final state particles. This time scale is given by the inverse of the kinetic energy of the final state particles [20,[26][27][28]. Consequently, for the case that the final state is a particle-antiparticle pair or are two same particles, with mass m 2 and decay width Γ 2 , one may introduce a cut-off velocity ∼ Γ 2 /m 2 , below which the FSS effect is suppressed.
In fact, a similar velocity cut-off also applies to the initial state Sommerfeld (ISS) effect between two coannihilators, which can convert to DM particles through scatterings and decays. Suppose that there is a direct coupling among a DM, a coannihilator and a massless particle, at the leading order the scattering rate takes the form of ∼ c scatt T 3 /m 2 1 , while the decay rate takes the form of ∼ c decay ∆m, where m 1 is the DM particle mass, T is the temperature, ∆m is the mass difference between the coannihilator and the DM particle, and the massless parameters c scatt and c decay depend on the couplings and mixings in a specific BSM model. While the time scale of scatterings is usually larger than that of the ISS effect, we find that for not very small ∆m the coannihilator decay time can be smaller than the latter and the ISS effect can be suppressed. Since ∆m determines the signatures of collider search for WIMP, the modification of DM relic abundance due to the suppression of the coannihilators' Sommerfeld effect may need to be taken into account, when one is looking for viable and interesting parameter regions in BSM models.
The rest of our paper is organized as follows. In Section 2, we calculate the thermally averaged FSS factor in WIMP DM annihilation during thermal freeze-out, estimate its impact on DM relic abundance, and discuss the modification induced by decays of the final state particles. In Section 3, we estimate in DM coannihilation scenarios the effect of the coannihilators' instability on the ISS factor between two coannihilators. We summarize our conclusions in Section 4.

Final state Sommerfeld effect
The origin of the FSS effect is the same as the more familiar initial state one, namely, the plane wave functions are modified due to some long-range force between non-relativistically moving final state particles. To illustrate the physics, in this paper we consider a simple Coulomb-like potential between two final state particles in a 2-to-2 annihilation, and we focus on the s-wave Sommerfeld effect. We assume that the two incoming particles have the same mass, m 1 , and the two outgoing particles have the same mass, m 2 .

Thermally averaged FSS factor
Following the notations in [4], the product of the annihilation cross section and the relative velocity of the incoming particles in the center-of-mass frame is written as (σv) w/o FSS = a v 2 , where a is a constant and v 2 is the velocity of one of the final state particles. The subscript 'w/o FSS' indicates that the FSS factor is not included. The dependence of v 2 comes from doing the phase space integration of the two outgoing particles. v 2 relates with v, m 1 and m 2 through where z ≡ m 2 /m 1 . For z ≤ 1, the reaction is kinematically allowed for any value of v, and the minimum value of v 2 is √ 1 − z 2 . For z > 1, the reaction is kinematically allowed only for v ≥ 2 1 − 1/z 2 , and the minimum value of v 2 is 0. Therefore, when m 2 is not too much smaller than m 1 , say, z 0.8, the final state particles can move with non-relativistic velocities, considering that the incoming particles move non-relativistically during the freeze-out and thereafter.
A Coulomb-like potential between the two non-relativistically moving final state particles, V (r) = −α f /r, modifies the otherwise two-body plane wave function and leads to the FSS effect. Depending on the charges or quantum states 2 of the final state particles under this long-range interaction, α f can be either positive or negative. The Sommerfeld-corrected s-wave cross section takes the form of (σv) with FSS = a S f v 2 , where the FSS factor is The thermally averaged cross section is where v min = 0 for z ≤ 1, and v min = 2 1 − 1/z 2 for z > 1. f (v) is the Maxwell-Boltzmann distribution for the relative velocity v, given as , and the S f v 2 in Eq. (3) can be written as Since S f is applicable to a non-relativistic v 2 , we choose to turn off the FSS effect for v 2 > 0.6. That is, when using Eq. (5) in the following calculations, we substitute where H(0.6 − v 2 ) is the Heaviside step function. This choice means that we only consider the FSS effect for z ≥ 0.8. If not considering the FSS effect for all v 2 , i.e., in the α f → 0 limit, the thermally averaged cross section becomes for which the integration can be done analytically, as shown in Eqs. (24) and (25) of [4]. The left panel of Fig. 1 plots S f v 2 as a function of z at a typical freeze-out value of x = 25, for several choices of α f . The solid blue, brown and purple lines are for α f = + 0.02, + 0.1 and + 0.5, respectively. The dashed blue, brown and purple lines are for α f = − 0.02, − 0.1 and − 0.5, respectively. A positive (negative) α f results in a Sommerfeld enhancement (suppression) for the thermally averaged cross section. The orange line in the middle is for α f = 0, that is, v 2 ; it is the same as the solid line in the upper panel of Fig. 2 in [4]. All lines merge at z = 0.8, since below this value the minimum v 2 is larger than 0.6 and we choose to turn off the FSS effect, as mentioned above. During the freeze-out and thereafter, typical incoming particles are non-relativistic, so that   To see the Sommerfeld effect more clearly, we show in the right panel of Fig. 1 the ratio of S f v 2 to v 2 , as a function of z. That is, at each z, the values on the solid and dashed blue, brown and purple lines are the corresponding ones in the left panel divided by the value on the orange line. Starting from z = 0.8, for a given α f the FSS effect increases with the increase of z until the latter is slightly bigger than 1, and after that it mildly decreases. This behavior can be understood from the fact that for z > 1 the fraction of the incoming particles that have enough kinetic energy to activate the reaction becomes smaller for larger z, and that the deviation of the FSS factor from 1 is larger for smaller v 2 according to Eq. (2), together with the above-mentioned fact that v 2 decreases with the increase of z. At x = 25 and for z = 1, there is a ∼ 15% Sommerfeld enhancement for the thermally averaged cross section for an electroweak interaction sized α f = + 0.02, a factor of ∼ 2 enhancement for a strong interaction sized α f = + 0.1, and a factor of ∼ 7 enhancement for an even stronger interaction α f = + 0.5; while for α f = − 0.02, − 0.1 and − 0.5, the suppression factor is about 0.87, 0.5 and 0.03, respectively. Fig. 1 gives us a hint that the FSS effect may bring a noticeable correction to the DM relic abundance. After the freeze-out, the DM density can continue decreasing by a factor of a few or even an order of magnitude and more. Therefore, we need to study the Boltzmann equation which describes the evolution of the DM density. By introducing the yield, which is the ratio of the DM number density to the entropy density, Y 1 ≡ n 1 /s, the Boltzmann equation can be written as

FSS effect on DM relic abundance
where G N is the gravitational constant, H(T ) is the Hubble parameter, g * s and g * are the numbers of effectively massless degrees of freedom associated with the entropy density and the energy density, respectively. Y 1,eq is the equilibrium value of the yield, Y 1,eq = n 1,eq /s, where g 1 is the DM degrees of freedom, and K 2 (x) is the modified Bessel function of the second kind. At x 1, n 1,eq ≈ g 1 3/2 e −m 1 /T . By integrating Eq. (7) from a small value of x before the freeze-out when Y 1 = Y 1,eq , to its current value which is essentially x → ∞, we get the yield of today, Y 1,0 . The DM relic abundance is given as [29][30][31] 3 where h is the present-day dimensionless Hubble parameter. We note that by writing the term Y 2 1,eq in Eq. (7), we have assumed the usual condition that the annihilation products quickly thermalize and that their number densities equal the thermal equilibrium values. This assumption ensures that the information about the final state particles' number density n 2 and its equilibrium value n 2,eq does not appear in the Boltzmann equation, due to the principle of detailed balance, σv → n 2 1,eq = σv ← n 2 2,eq , where σv → ≡ σv is the forward reaction and σv ← is the backward reaction. While this assumption is true for SM final products, one might need to check its validity when the annihilation products are new particles in BSM models. Without committing to specific BSM models, in this work we make this assumption and study the FSS effect in the simplest situation.
We are interested in the ratio of DM relic densities with and without including the FSS effect, Ω/Ω w/o S f . A good estimation of the ratio can be obtained by using an approximate form of Eq. (7), which is valid after the freeze-out, since Y 1,eq quickly becomes negligible compared to Y 1 . Eq. (11) has the solution where x fo is the freeze-out value of x, defined when Y 1 (x fo ) − Y 1,eq (x fo ) ≡ κY 1,eq (x fo ) for a numerical constant κ of order unity. Since the final yield Y 1,0 is usually much smaller than Y 1 (x fo ), we can make a further approximation to get Therefore, the ratio we want is in which we have dropped the factor 1 + T 3g * s dg * s dT in the integrals, since | T 3g * s dg * s dT | 1 except in the narrow range around the quark-hadron phase transition temperature ∼ 200 MeV [29,32]. Unless m 1 /x fo falls into this temperature range, neglecting this factor is of no harm to our estimation of the ratio.
In Eq. (14) the lower limits of the two integrals, x fo,y and x fo,n , are x fo with and without the inclusion of the FSS effect. Since (x fo + 1 2 ln x fo − ln σv x fo ) is roughly a constant (see e.g. Eq. (2) of [4]), where σv x fo means the value of σv evaluated at x fo , then Since |x fo,y − x fo,n | x fo,n , up to the first order of (x fo,y − x fo,n )/x fo,n we get where S f v 2 x fo,n means the derivative of S f v 2 with respect to x, evaluated at x fo,n . Using this relation between x fo,y and x fo,n , we see from Eq. (14) that an estimation of Ω/Ω w/o S f only depends on the values of x fo,n , z and α f . From the solid lines, we can see that the FSS effect results in a more than ∼ 5% (23%, 65%) reduction to the relic abundance for α f = 0.02 (0.1, 0.5) for z ≥ 0.85, and the reduction reaches ∼ 18% (58%, 90%) at z = 1; for α f = − 0.02 (− 0.1, − 0.5), the FSS effect increases the relic abundance by more than ∼ 6% (32%, a factor of ∼ 4.7) for z ≥ 0.85, and the value reaches ∼ 22% (2.6, 51) at z = 1. Therefore, indeed the FSS effect can be important in DM relic abundance calculations, and an electroweak interaction sized α f is already large enough to make a correction well beyond the percent level accuracy of the observational value.
In order to check the estimation, we start from Eq. (7) and choose three sets of values of (a , m 1 ) to calculate DM relic densities with and without the FSS effect, and then we get the ratios. The x fo,n =25 x fo,n =30 x fo,n =25 x fo,n =30  7), using three sets of (a /GeV −2 , m 1 /GeV) = (10 −9 , 10 3 ), (10 −9 , 10 4 ), (10 −7 , 10 4 ), and g 1 = 2. According to whether the freeze-out, without including the FSS effect, happens before m 1 /T = 22.5, between 22.5 and 27.5, or after 27.5, these ratios are marked with blue dots, orange triangles, or red squares, respectively. sets are (10 −9 GeV −2 , 10 3 GeV), (10 −9 GeV −2 , 10 4 GeV) and (10 −7 GeV −2 , 10 4 GeV). For all sets, we use g 1 = 2 in the calculations. We compute the ratios for α f = ± (0.02, 0.1, 0.5), and from z = 0.8 to 1.2 with an interval of 0.04. In Fig. 2, we mark each ratio according to the freeze-out value of x without including the FSS effect, given by when Y 1 /Y 1,eq = √ 2, as used in [4]. A ratio is marked with a blue dot if this freeze-out value of x is less than 22.5, an orange triangle if it is between 22.5 and 27.5, and a red square if it is larger than 27.5. We expect that the marks should fall closely to the lines computed by using the same α f . We also expect that for each α f , the blue dots, orange triangles and red squares should fall closer to the dashed, solid and dotted line, respectively. We can see that it is indeed the case, and this confirms our estimation formula Eq. (14).

Correction to the FSS effect induced by final state particle decays
Unlike the initial state DM particles which are stable or very long lived, the final state particles may decay. Consequently, we need to consider whether there is enough time for the FSS effect to happen before final state particles decay. The former time scale is determined by the typical long-range interaction time between the two final state particles, ∼ (m 2 v 2 2 ) −1 , while the latter is the inverse decay rate Γ   We show in the left panel of Fig. 3 the final state particle decay effect on the evolution of S f v 2 / v 2 with x. Since the behaviors of this quantity are different on the two sides of z = 1, as shown in Fig. 1, we study its evolution for three z values, 0.9, 1 and 1.1, represented by the green, orange and red lines, respectively. The solid (dashed) lines above (below) S f v 2 / v 2 = 1 are for α f = 0.1 (− 0.1); these lines are all for v 2 cut = 0, that is, the final state particle decay effect is not included. The two dotted lines and the two dash-dotted lines are for (α f = 0.1, v 2 cut = 0.1) and (α f = − 0.1, v 2 cut = 0.1), respectively. We can see that the deviations of the dotted (dash-dotted) lines from the solid (dashed) lines grow with the increase of x. This can be understood from the fact that typical velocities of final state particles decrease with the increase of x. For a given α f , the FSS factor S f deviates from 1 more significantly for smaller v 2 , but particles with v 2 < v 2 cut do not contribute. These two opposite v 2 effects compete, and the deviations of dotted and dash-dotted lines from S f v 2 / v 2 = 1 cease to grow and then reduce around x ∼ 60. We note that only the the latter one. 5 By using the Heaviside step function H(v 2 − v 2 cut ), the FSS factor becomes 1 for v 2 < v 2 cut . However, we note that the suppression of the FSS factor from v 2 > v 2 cut to v 2 < v 2 cut is a smooth transition. Nevertheless, a sharp cut-off is sufficient for an estimation of the effect of final state particle decays on the FSS factor. Similarly, in Sec. 3 we also use a sharp velocity cut-off for an estimation of the effect of the coannihilator decays on the ISS factor. cases z = 1 and 1.1 have the v 2 cut = 0.1 lines. For the case z = 0.9 the minimum v 2 is larger than 0.1 so that H(v 2 − v 2 cut ) always equals 1 for v 2 cut = 0.1.
In the right panel of Fig. 3 we show the impact of final state particle decays on the ratio of DM relic densities with and without the FSS effect. The lines above and below Ω/Ω w/o S f = 1 are for α f = − 0.1 and 0.1, respectively. x fo,n = 25 is used for all lines. The solid lines are for v 2 cut = 0, and they are the same as the middle solid lines in the two panels of Fig. 2. The dashed, dotted and dash-dotted lines are for v 2 cut = 0.03, 0.1 and 0.3, respectively. For z < 0.995, the dotted and solid lines completely overlap, since the minimum value of v 2 is larger than 0.1. For the same reason, the dash-dotted and solid lines completely overlap for z < 0.954. We can see that the deviations of the dotted and dash-dotted lines from the solid lines quickly grow until z = 1, and then gradually reduce for larger z. At z = 1, for α f = 0.1, compared to the ∼ 58% reduction of the relic abundance when not considering decays of the final state particles, the values are ∼ 46% and only 6% for v 2 cut = 0.1 and 0.3, respectively; for α f = − 0.1, the FSS effect increases the relic abundance by a factor of ∼ 1.9 and only 5% for v 2 cut = 0.1 and 0.3, respectively, compared to a factor of ∼ 2.6 increase for v 2 cut = 0. On the other hand, the deviations of the dashed lines from the solid lines are quite small.
Similar to what we did in Fig. 2, we check the results by calculating the FSS-corrected DM densities from Eq. (7), imposing v 2 cut = 0.1 or 0.3. We use (a = 10 −7 GeV −2 , m 1 = 10 4 GeV) for the orange triangles in the range 1.04 ≤ z ≤ 1.2, and (a = 10 −9 GeV −2 , m 1 = 10 4 GeV) for the rest, because for each of these points the freeze-out value of x is between 22.5 and 27.5 when the FSS factor is not included. All the orange triangles match well the corresponding lines, and this confirms the results.
We conclude that the modification of the FSS effect induced by final state particle decays is significant when Γ 2 /m 2 O(10 −1 ), and it is already not negligible when Γ 2 /m 2 ∼ O(10 −2 ).

Modification of the initial state Sommerfeld effect in DM coannihilations induced by the coannihilators' instability
As a corollary, we investigate in DM coannihilation scenarios whether there are also noticeable corrections to the Sommerfeld effect between coannihilators induced by their instability. We consider the usual situation that the coannihilators and the DM share the same discrete symmetry which makes the DM stable (e.g., the R-parity in supersymmetric models [33] and the KK-parity in Universal Extra Dimension models [34]), so that a coannihilator can convert to the DM particle or other species of coannihilators through decays or scatterings. In order to have the two-body wave function of a pair of coannihilators modified from the plane wave, the two particles need to come together from an initial separation large enough compared to the characteristic distance of the long-range interaction, ∼ (m c v c ) −1 , where m c is the coannihilator's mass and v c is one of the coannihilators' velocity in their center-of-mass frame. The initial separation is give as ∼ v c /(Γ decay + Γ scatt ), where Γ decay and Γ scatt are the decay and scattering rates. If the initial separation is larger, a coannihilator would have decayed or converted away before it meets the other one 6 . Both the decay and scattering rates are model dependent 7 . In this work, we consider the simplest situation by assuming that there is only one species of coannihilator, and that there is a direct coupling among a coannihilator, a DM and a massless particle, so that at the lowest order the rates can be written as Γ decay ∼ c decay ∆m and Γ scatt ∼ c scatt T 3 /m 2 1 , where ∆m is the mass difference between the coannihilator and the DM, ∆m ≡ m c − m 1 , c scatt and c decay are massless parameters depending on the couplings and mixings in a specific BSM model. Since c scatt and c decay are controlled by the same interaction vertex, they usually do not differ by orders of magnitude. Then Γ scatt is typically smaller than Γ decay during and after the freeze-out, unless the DM and the coannihilator are very degenerate in mass. We will only focus on Γ decay in this work.
Similar to the FSS effect case, we estimate the impact of the instability of the coannihilators on the ISS effect by applying a velocity cut-off on v c , given as v c cut = c decay ∆m/(m 1 + ∆m), below which we switch off the Sommerfeld factor between two coannihilators. We consider the s-wave annihilation cross section between a pair of coannihilators, and assume that there is a Coulomb-like potential between them, given as V (r) = −α i /r. Similar to Eqs. (2) and (3), the ISS factor is and the thermally averaged cross section is where a cc is the s-wave cross section without considering the ISS factor. Note that we do not need to write the term in the square bracket as [ , as we did for the FSS case, because during the freeze-out and later the incoming coannihilators with relativistic velocities are in the tail of the Maxwell-Boltzmann velocity distribution 8 . For simplicity, we assume that for the parameter space we will consider, the effective annihilation cross section σv ef f is always dominated by the above cross section, and we neglect the contributions from the DM-DM and DM-coannihilator (co)annihilation cross sections. That is, we use where g c is the degrees of freedom of coannihilators. g ef f is the effective degrees of freedom, given as Also, if the coannihilator is not its own antiparticle, then a cc should be understood as a cc , and we neglect the particle-particle and antiparticle-antiparticle cross sections in that case. The DM relic abundance can be solved by substituting in Eq. (7) σv by σv ef f , Y 1 by Y ≡ Y 1 + Y c and Y 1,eq by Y eq ≡ Y 1,eq + Y c,eq , where Y c,eq = n c,eq /s = T 2π 2 g c m 2 c K 2 (m c /T )/s; then, in Eq. (10) substitute Y 1,0 by Y 0 . Fig. 4 shows the impact of the coannihilators' instability on DM relic abundance. As an example, we choose g 1 = g c = 2, m 1 = 2 × 10 3 GeV and a cc = 10 −8 GeV −2 .
coannihilation may not happen, if the squark mass is very heavy [13,35]. 8 We have explicitly checked that the curves in Fig. 4 do not change if we impose the H(0.6 − v c ) factor.  In the left panel, we fix α i = 0.1, and plot the ratio of DM relic densities with and without the ISS effect, as a function of ∆m, for several choices of c decay . Ω w/o S i h 2 is calculated by letting S i = 1 in Eq. (19), and it is about 0.096 and 1.2 at ∆m = 0 and 100 GeV, respectively. The solid line is the ratio when not considering the instability of the coannihilator, namely, c decay = 0 so that v c cut = 0. The dashed, dotted and dash-dotted lines are for c decay = 0.01, 0.1 and 1, respectively. The deviations of these lines from the solid line are getting larger with the increase of ∆m, which controls the decay rate of the coannihilator for a given c decay . The impact of the instability on the ISS effect is small for c decay = 0.01, while it becomes visible for c decay = 0.1 or larger. For instance, without considering the instability of the coannihilator, at ∆m = 40 GeV the ISS effect can reduce 59% of the DM relic abundance, while the reduction becomes 56% for c decay = 0.1 and 36% for c decay = 1. At ∆m = 100 GeV, compared to the 54% reduction when taking c decay = 0, the number becomes only 18% if c decay = 1.
In the right panel, we fix ∆m = 100 GeV, and plot the ratio of DM relic densities with and without considering the instability of the coannihilator, as a function of c decay , for several choices of α i . The blue, red and purple lines are for α i = 0.02, 0.1 and 0.5, respectively. All lines tilt up because the strength of the ISS effect is more suppressed when the coannihilators became increasingly unstable. The effect of the instability is more significant for larger α i . For example, compared to a stable coannihilator, an unstable coannihilator with c decay = 0.4 can make the DM relic abundance larger by 9%, 39% and 73% for α i = 0.02, 0.1 and 0.5, respectively.
We conclude that when the DM and the coannihilator are not very degenerate in mass, for Γ decay /∆m O(10 −1 ), the modification of the ISS effect induced by the decay of the coannihilator may need to be considered for an accurate calculation of the DM relic abundance. Such modification can be quite significant when the long-range force between the coannihilators is of the strong interaction or stronger size for ∆m/m 1 1/20.

Summary
We have studied in this paper the final state Sommerfeld effect on DM relic abundance. This effect occurs when the DM annihilation products move non-relativistically and there is some long-range force between them, so that the wave function of the final state particles is modified from the plane wave. As a proof of concept, we consider the case that two WIMP DM particles annihilate into two equal mass particles, and calculate the thermally averaged s-wave FSS factor arising from a Coulomb-like potential between the two final state particles. We show the dependence of the FSS effect on the strength of the long-range interaction, as well as on the mass ratio of the final and initial state particles. We find that the impact of the FSS effect on DM relic abundance can be significant, and an electroweak sized long-range interaction is already large enough to make a correction well beyond the current percent level observational accuracy.
While the physical origin of the FSS effect is similar to the well-studied initial state Sommerfeld effect between two stable annihilating DM particles, an additional point of the former is that the final state particles are unstable, so that the FSS effect may not have enough time to happen before the particles decay. We find that when the mass ratio of the final and initial state particles is close to 1 and larger, the decay can suppress non-negligibly the FSS effect if the decay rate is larger than 1% of the final state particle's mass, and if larger than 10% the suppression becomes significant.
As a corollary to the above point, we study the impact of the instability of the coannihilators on the initial state Sommerfeld effect in the calculations of the DM relic abundance in the coannihilation scenario. Here the instability comes from decays and scatterings of a coannihilator into the DM and other species of coannihilators. Since the decay rate usually dominates over the scattering rate unless the DM and the coannihilator are very degenerate in mass, we focus on the decay rate in this work. We find that the decay of the coannihilator makes a non-negligible correction to the DM relic abundance when the decay rate is more than 10% of the mass difference between the coannihilator and the DM particle. If the long-range interaction between the coannihilators is of the strong interaction size or more, the correction can be quite large when the mass difference is over ∼ 1/20 of the DM mass.
Before we close, we note that other types of long-range force, for example a Yukawa potential, can also give rise to the FSS effect when the final state particles move non-relativistically. Since we consider two equal mass final state particles in this work, the condition of non-relativistic moving is satisfied when the final state particle's mass is close to the initial state one. However, for the situation that the masses of the two final state particles are different, or for the final states in for instance 2-to-3 annihilation, the parameter space of at least two final state particles being nonrelativistic will be different. Also, besides its effect in the calculations of the DM relic abundance in the early Universe, the FSS effect may play a role in the indirect searches for DM in the late Universe. Finally, for the impact of the instability of the coannihilators on the ISS effect, while we focus on the decay of the coannihilators in this work, the scattering is expected to be the dominant way of the coannihilator ↔ DM conversion when the coannihilator and the DM are very degenerate in mass. This scenario may be worth to be explored, since it has interesting collider signals (see e.g. [35][36][37][38][39]) and it may help answer the question of how heavy the DM can be in the WIMP DM coannihilation scenarios [13,16].