On the classical description of the recombination of dark matter particles with a Coulomb-like interaction

Cold dark matter (DM) scenario may be cured of several problems by involving self-interaction of dark matter. Viability of the models of long-range interacting DM crucially depends on the effectiveness of recombination of the DM particles, making thereby their interaction short-range. Usually in numeric calculations, recombination is described by cross section obtained on a feasible quantum level. However in a wide range of parameter values, a classical treatment, where the particles are bound due to dipole radiation, is applicable. The cross sections, obtained in both approaches, are very different and lead to diverse consequences. Classical cross section has a steeper dependence on relative velocity, what leads to the fact that, after decoupling of DM particles from thermal background of"dark photons"(carriers of DM long-range interaction), recombination process does not"freeze out", diminishing gradually density of unbound DM particles. Our simplified estimates show, that at the taken parameter values (the mass of DM particle is $100$ GeV, interaction constant is $100^{-1}$, and quite natural assumptions on initial conditions, from which the result is very weakly dependent) the difference in residual density reaches about $5$ orders of magnitude on pre-galactic stage. This estimate takes into account thermal effects induced by dipole radiation and recombination, which resulted in the increase of both temperature and density of DM particles by a half order of magnitude.

is very different and whether or not the model gets constraint from observations [3,12,13]. Description of recombination process plays a clue role here. Usually a quantum approach is used for it. However a classical approach, which was used, in particular, for magnetic monopoles and heavy neutrinos [22,23], seems to be valid in a broad interval of parameter values. It leads to a result very different from that obtained on quantum level, which in the commonly accepted form does not come to the classical limit.
Classical recombination cross section is obtained from condition that the scattered particles lose, due to dipole radiation, sufficient energy to get bound [24,25], and is given by where ρ max is the maximal impact parameter at which a pair is bound, α y is the constant of y-interaction, µ = mam b ma+m b is the reduced mass of the pair of the scattered particles with m a and m b being their masses (m a ≤ m b ), v is their initial relative velocity. This cross section has a steeper dependence on velocity with respect to that of usually accepted quantum recombination cross section. The latter in form of Kramers' formula [26] summed over all quantum levels is (valid for v α y ): Steeper behaviour of classical cross section leads to the fact that recombination process does not freeze out on both radiation dominated (RD) and (even faster expanding) matter dominated (MD) stages, and relative number of unbound y-charged particles falls down gradually with time.
Under this condition, binding is found to occur predominantly due to multiple soft photon emission, what allows classical treatment. However considering on quantum level, only one-, two-photon final states are usually taken into account. Eq.(3) can be formally deduced from condition that binding of two particles (i.e. when initial kinetic energy of relative motion, E rel = µv 2 /2, is lost) happens at distance (R b ) much greater than the radius of the respective ground bound state (a B ). In this case, action of the system, as will be shown, becomes much greater 1 (in units = 1) and thus the classical approach is reasonable. Below we shortly discuss classical approach implications in recombining DM cosmological evolution. It includes estimations, some of which may seem to be simple, to trace explicit dependence of result on the parameters.
Our results will basically relate to parameter region lying around fiducial values m a ∼ µ = 100 GeV and α y = 1/100 used for numerical estimations. Also for definiteness we assume (as in case of heavy neutrino model [22]) that before a direct annihilation of a and b (where either b =ā or b =ā), happening when the temperature becomes below their mass (of the lightest), T = T * ∼ m a /10, y-plasma (consisting of a, b with their antiparticles and y-photons -massless y-force carriers) has the same temperature as ordinary matter. Right after annihilation, y-photon background decouples from a and b as well as from ordinary matter (O), while the opposite is not true. That is the y-background is no longer influenced by a and b, and by O through O-a, b possible coupling, but a and b are influenced by y. Starting from this moment, the temperature of y-background (T y ), as of a closed system, changes as inverse scale factor, whereas that of ordinary photons feels also entropy re-distribution between ordinary matter components. So for y-and O-photons' temperature relation one has where g s,o (T ) is the effective number of ordinary matter species (excluding y) contributing into entropy density (s). For the chosen numeric values, contribution of y in density at nucleosynthesis (BBN) makes up κ 4/3 (T ∼ 1 MeV) ≈ 0.06 from that of O-photons, what has no effect on BBN data. The temperature of a and b is equal to T y until they decouple from y, influence of O-matter (being determined by some weak scale interaction) is negligible (see Appendix A). After decoupling it evolves as (see Appendix A) Here g , being effective number of matter species (including y) contributing into energy density, as well as κ are taken at T = 0.2 MeV. Evolution of number density of unbound a-and b-species can be approximately described by equation system (see Appendix B) Here r is the number density conventionally expressed in units of s, the brackets " " mean averaging over velocity distribution of a and b, θ shows deviation of T a from Eq.(5), E pair = E a + E b is the energy lost by a-,b-gas (thermal bath) in the result of the pair binding (that is the case of annihilation, i.e. when b =ā, or when the bound systems are thermally decoupled from a-, b-gas). The second equation takes into account thermal effects, caused by scattering of particles: (*) presumably slower particles to be bound "go out" (annihilate or decouple) of a-,b-plasma, effectively heating it; (**) scattered but unbound pairs experience dipole energy losses cooling plasma. Evolution is considered in terms of O-photon temperature (T ).
Here we do not take into account some recombination process details which are more appropriate for quantum case (such as recombination into different level bound states, red-shifting of recombination photons), as well as inverse processes (which are not important for the big parameter space of question), and also in quantum case the second equation of (6) is omitted (θ ≡ 1).
The bound states start to form when T a becomes much lower than ionization potential I = µα 2 y /2, T a < T a rec ∼ I/10. For the chosen values, it is T a rec = 0.5 MeV, at which T a = T y , what corresponds to O-matter temperature T = T rec = κ −1/3 T a rec ≈ 1 MeV. Depending on behaviour of σ rec v with T and epoch, recombination process should flow in different regimes. Basically, it either damps (freezes out) and it is initial moment (when it starts) what predetermines the residual density of free a and b, or it "burns" continually with a self-adjusted rate and final moment, at which we need to know the density, defines its value. For classical cross section at T ∼ T rec recombination process gives effect and temporarily freezes out, but after a-y decoupling it is restored and goes with a steady rate until the galactic stage 1 . Herewith, recombination rate Γ rec = n σ rec v turns out to be of the same order as Hubble rate (H), if thermal effects (second equation of system (6)) are ignored (θ ≡ 1), otherwise (θ = 1) the ratio Γ rec /H slowly goes down with time. So the value of residual density is not sensitive to the initial moment and all early history, but fully defined by the final one, and is just weakly dependent on initial abundance for θ = 1 (see below).
Classical approach is assumed to be valid when R b as well as energy loss length-scale 2 , L loss , are much greater than dark atomic size a B . From other side, they must be much less than spacing between y-charged particles, L sp (a fortiori a-a(b), y interaction lengths). R b and L loss can be found by taking mechanical Here M is the angular momentum, which can be assumed to be conserved. In the region of interest, solution of Eq.(7) can be simplified by neglecting "E" in the square root. Condition R b , L loss L sp is found to be true by a wide margin for the most of parameter values at any epoch of question [27]. The values R b and L loss are of the same order of magnitude and both depend on impact parameter ρ and on energy, which falls down with the Universe expansion. Dependences of R b /a B from ρ/ρ max and from the temperature of O-photon are shown in Figs (1(a), 1(b)). One can see that R b /a B at the most of ρ values is close to that at ρ = 0, therefore R b /a B (ρ = 0) is used in all further calculations.
Since the late period is determinative for residual density value in case of classical cross section, we single out late expansion moment, when the galaxies start to form, z ∼ 10 (T ∼ 30 K). As one can see from Fig (1(b)), R b /a B 1 there. This ratio depends on parameters as R b /a B ∝ α 2 y µ 2/5 /T  (1) at T < T q-c . As one can see, in the late period this solution comes to that obtained in pure classical approximation. A small deviation is caused by thermal effects (described by second equation of the system (6)) which make the result being a little sensitive to the initial density (see Eq. (27) of Appendix B). In calculations the following values had been used: α y = 1/100, m a = 100 GeV, m b = 1 TeV, r(T rec ) = ρ CDM /m b /s(T = 2.7 K) ≈ 4.6 × 10 −13 with ρ CDM ≈ 1.4 keV/cm 3 being the modern CDM density. Also, quantum recombination rate was taken from [1], where σ rec is a little corrected as compared to Eq. (2).
Note that the mentioned above thermal effects in classical case are found to be weak (θ(T ) ∼ T −1/23 ), but nonetheless noticeable. At T ∼ 30 K for µ = 100 GeV and α y = 1/100, θ ≈ 4. The ratio of r(T ∼ 30), obtained with and without account for thermal effects (using only classical cross section), is ∼ 3.
Condition (3) for period z ∼ 10 allows to outline the region of parameters values (µ and α y ) when the classical approximation for σ rec is applicable. It is shown in the Fig.(2(b)). There can exist a region of parameter space where calculations should be re-considered, since interaction radius ∼ α y /T a exceeds spacing between DM particles [27].
Finally, we show arguments in favour of condition Eq.(3) on the base of action in the Appendix C.
To conclude, we have shown that classical recombination cross section Eq.(1) can be applicable in a broad parameter region of self-interacting DM models, and it leads to strongly different result comparing to obtained with usually accepted quantum cross section. It may change (extend) the parameter region of the corresponding models' viability.

Acknowledgements
We would like to thank M. Yu. Khlopov, cooperative work with whom initiated this paper, S. G. Rubin and A. E. Dmitriev for interest and useful discussions, M. N. Laletin for his help in writing the text.
This work was performed within the framework of Fundamental Interactions and Particle Physics Research Center supported by MEPhI Academic Excellence Project (contract № 02.a03.21.0005, 27.08.2013), grant of RFBR (№14-22-03048) and in part by grants of Ministry of Education and Science of the Russian Federation (№ 3.472.2014/K) and the Russian Science Foundation (№ 15-12-10039).

Appendix A
Here we trace shortly the temperature evolution of y-charged particles, a (b), before recombination starts. They are assumed to experience energy exchange with y-background and, possibly, ordinary matter. In calculations we suppose that a interacts with O as Dirac heavy neutrino.
To find T a one can formally use the first law of thermodynamics, which for a-particles can be reduced to where o = e, ν, p, n, ... are available species of O-matter, ∆Eσv ai is the kinematically averaged energy transfer in ai-interaction multiplied by respective cross section and relative velocity, averaged in thermal distribution, n i is the respective number density (i = y or o).
In all scattering processes of interest a is non-relativistic, i is ultra-relativistic (p, n turn out to be too suppressed in density, so play no role in heat transfer with a). We adopt Boltzmann approximation for all species distributions. The ay-scattering is well described by Thomson cross section over the great part of parameter space. Then calculation of ∆Eσv ay with accuracy ∼ T /m a gives where ω and E are the energies of y-and a-particles respectively. To estimate ∆Eσv ao , assuming that a is heavy neutrino, one needs to take the cross sections of relevant processes (aν, aν, ae − , ae + ), which in our limit are undistinguishable for particle and antiparticle Here ω lab is the energy of incident i-particle in the reference frame where a is in the rest, G F is the Fermi constant, ξ e = 1 − 4ξ + 8ξ 2 ≈ 0.50 with ξ = sin 2 θ W being the weak mixing parameter. Unlike ay-scattering, ae-, aν-cross sections depend on energy, what accounts for higher power of the temperature in the final expression Eq.(8) can be then re-written as T ν = T e = T was put. Coefficients T ay ∼ 0.1 MeV, T ao ∼ 10 MeV. Second term in the r.h. of Eq.(12) has no effect on solution. Excluding it, solution can be expressed in the form At T T ay (i.e. after a-y decoupling) Appendix B Evolution of abundance of free a-, b-particles can be described by Boltzmann equation which is easily reduced to the first equation of the system Eq. (6) with the help of replacements: n = rs = r 2π 2 gs 45 T 3 , −dt = 1 H dT T (for g ,s ≈ const). For calculation one parametrizes σ rec = σ 0 /v β . For σ rec v then one can get Since recombination rate is strongly temperature dependent (especially in classical case), thermal effects of recombination process itself can be important, correcting temperature evolution. These effects are relevant, obviously, after a-y decoupling.
One can take again the first law of thermodynamics, dQ = δA + dU . One has the total particle number in some volume N ab = n ab V with n ab = n a + n b = 2n being their number density, the pressure p = n ab T a . Expansion of the Universe is treated as a work of gas: δA = pdV = n ab T a 3HV dt. Inner energy gain is dU = 3 2 d(N ab T a ) = 3 2 N ab dT a + 3 2 T a dN ab . Here we assume that dN ab = −2 σ rec v n 2 V dt as if the recombined pairs disappear from a-, b-gas. It would be true when b =ā, and also when b =ā if the bound systems are out of thermal equilibrium with free a and b, however it is not always the case [12]. We do not see the opposite since it is not of principle here.
One can define dQ through the energy lost during ab-scattering. If the pair is combined (impact parameter ρ < ρ max ), then their energy E pair = E a + E b = mav 2 is lost completely. Otherwise (ρ > ρ max ), their energy is lost partially due to dipole radiation, what is for given ρ and v = | v a − v b | (in large scattering angle limit, which is realized when ρ α y /(µv 2 )) [24] 3 where E rel = µv 2 /2 is the energy of relative motion. So, energy losses rate by ab-gas per unit volume can be given bẏ ε = n a n b where f a,b is the distribution in velocity (Maxwell). Upper limit D should be given by Debye length of ab-plasma, but thanks to fast convergence of ∆E(v, ρ) with ρ → ∞, we put D → ∞ in Eq. (18). Averaging over velocity distributions gives E pair σ rec v = (7 − β)Γ(2 − β 2 ) 2 for initial conditions θ(T 0 =T ay ) = 1, r(T 0 ) = r 0 . Being interested in r(T ) on MD-stage, we will put the solution of system for RD-stage to be initial conditions for that for MD-stage. Substituting r(θ) of Eq. (25) in this manner into second equation of the system (24) yields θ(T ) = 1 +γ β R D R r 0 1