Langevin simulation of dark matter kinetic equilibration

Recently it has been questioned, notably in the context of the scalar singlet dark matter model with $m_\varphi^{ }\simeq 60$ GeV, how efficiently kinetic equilibrium is maintained if freeze-out dynamics is pushed down to low temperatures by resonant effects. We outline how Langevin simulations can be employed for addressing the non-equilibrium momentum distribution of non-relativistic particles in a cosmological background. For a scalar singlet mass $m_\varphi^{ }\simeq 60$ GeV, these simulations suggest that kinetic equilibrium is a good approximation down to $T \sim 1$ GeV, with the deviation first manifesting itself as a red-tilted spectrum. This reduces the annihilation cross section, confirming findings from other methods that a somewhat larger ($<20\%$) coupling than in equilibrium is needed for obtaining the correct abundance.


Introduction
Determining the abundance of dark matter in a given model requires ingredients from general relativity, quantum field theory, and non-equilibrium statistical physics. While the first two are well-established frameworks, non-equilibrium statistical physics is an open construction, with even the specification of the state hard to achieve in full generality. Therefore it is not rare that cosmological computations call for physical intuition.
Often, a fruitful approach is to estimate the rates at which various processes take place, and then divide the variables into two classes, the fast and slow ones. A tractable problem is found if the fast variables undergo many interactions within the observation period; they can then be assumed to thermalize, constituting a heat bath. The slow variables remain out of equilibrium, but since there are fewer of them, the problem is easier to handle.
In dark matter computations, we normally assume that all Standard Model particles are fast variables. For dark matter particles, the fastest processes are soft elastic scatterings, which may be assumed to decohere the system. Subsequently we can focus on classical notions, like the momentum distribution of the dark matter particles ("kinetic non-equilibrium"), and their overall number density ("chemical non-equilibrium").
Adjusting the momentum distribution towards the thermal (Bose-Einstein or Fermi-Dirac) form only requires elastic scatterings. In the non-relativistic regime, these are much faster than inelastic ones. Therefore, it is often a good assumption to impose kinetic equilibrium from the outset, and only focus on deviations from chemical equilibrium.
Recently, however, the validity of this picture has been questioned. In particular, the example of the scalar singlet dark matter model (cf., e.g., refs. [1][2][3][4] and references therein) has been intensively discussed in the non-relativistic freeze-out regime [5][6][7][8][9] (cf., e.g., refs. [10][11][12][13][14] for similar effects in other models). To be clear, let us remark that in the so-called freeze-in scenario, dynamics takes place in the relativistic regime, and then there is in general no hierarchy between the kinetic and chemical equilibration rates, so that kinetic non-equilibrium is certainly present.
Following their use in the context of heavy ion collision experiments [15], we propose here to employ Langevin simulations for studying the efficiency of kinetic equilibration in the non-relativistic regime. The effect of the fast variables is encoded in the values of two matching coefficients, which can be defined and computed at the NLO [16] or even at the non-perturbative level [17]. Therefore Langevin simulations offer for a systematically improvable framework for studying strongly coupled systems, notably dark matter scattering off a Standard Model plasma at temperatures of a few GeV.
This paper is organized as follows. We start by reviewing how the Langevin equation can be set up in an expanding background, in sec. 2. This is followed by a description of an algorithm for its numerical solution, and a summary of the corresponding simulation results, in sec. 3. In sec. 4 we show how the non-equilibrium momentum distribution can be implemented in a freeze-out computation. Conclusions are collected in sec. 5, relegating the computation of the matching coefficients for the scalar singlet model to appendix A.

Langevin equation in an expanding background
We assume the universe to be described by a homogeneous, isotropic and spatially flat Friedmann-Lemaître-Robertson-Walker background, with the metric (2.1) The physical 4-momentum of an on-shell particle is denoted by p µ . The covariant derivative of its spatial components reads where H ≡ȧ/a is the Hubble rate andṗ ≡ dp/dt. Viewing p i as a slow variable, the Langevin equation takes the formṗ where η is a friction (or "drag") coefficient and f i is a random force, taking care of detailed balance (i.e. returning thermal energy to the heavy particle, in exchange for that lost through friction). The force obeys the autocorrelator where ζ is called the momentum diffusion coefficient. The constraint that the system should thermalize to a temperature T imposes the fluctuation-dissipation relation The average velocity can in turn be expressed as v 2 ≈ 3T /m ϕ , where we have introduced the notation m ϕ for the mass of a generic non-relativistic dark matter particle. When we implement the Langevin equation in a cosmological context, time and temperature are not independent variables. If the system does not undergo phase transitions, so that the temperature evolves smoothly, we may take as a time-like variable (we choose T max ≡ 5 GeV). The Jacobian to physical time is where c 2 s = ∂p/∂e is the speed of sound squared. Furthermore the entropy density, s, satisfieṡ s + 3Hs = 0, and consequently sa 3 = const. If we now define dimensionless momenta aŝ then Langevin dynamics can be expressed as Given the constancy of sa 3 , we note thatp i ∝ ap i ≡ k i , known as a comoving momentum.
A key element of the dynamics is that the coefficientsη andζ are not constant but evolve rapidly with x. The Hubble rate reads where e is the energy density and m pl ≈ 1.22091× 10 19 GeV is the Planck mass. Since e ∼ T 4 in the Standard Model plasma, H scales as ∼ T 2 . The entropy density scales as s ∼ T 3 . The coefficient ζ is suppressed by the mass of the dark matter particle and that of the mediator between the visible and dark sectors. For dimensional reasons, we may write it as with ξ displaying modest temperature dependence. Different contributions to ξ in the scalar singlet model, derived in appendix A, 1 are shown in fig. 1(left). The speed of sound squared  from scatterings off quarks or gluons, converted into ξ as defined by eq. (2.12), with m ϕ ≃ 60 GeV and κ ≃ 0.00064 [18]. Confining effects in the vicinity of the QCD crossover have been modelled via the substitution N c → N c,eff < 3 [19], and the difference to the tree-level value N c = 3 has been displayed as a grey band. Middle and right: the correspondingη andζ from eq. (2.9), either with ξ as displayed in the left panel (solid line), or for various fixed values of ξ. For thermodynamic potentials we have inserted estimates from ref. [20], tabulated at http://www.laine.itp.unibe.ch/eos15/.
can often be approximated as 3c 2 s ≃ 1, though it experiences corrections when mass thresholds are crossed. Putting all of these scalings together, we expectη ∝ (T /GeV) 4 and ζ ∝ (T /GeV) 3 (cf. fig. 1(middle, right)). This implies that kinetic equilibrium is likely to be lost at low T . Equation (2.10) is a linear inhomogeneous first order differential equation, and as such it can be given an explicit formal solution, (2.13) Taking an average over the noise, moments can be obtained, for instance A numerical illustration is shown in fig. 2. Compared with such moments, the advantage of a direct numerical solution of eq. (2.10) is that all moments can be obtained at once, from the momentum distribution. Finally we note that the would-be equilibrium distribution, at any given temperature, is obtained by assuming the coefficients temperature-independent, so that there is a lot of time for the system to adjust to the given situation. If we normalize the momentum distribution in analogy with cosmological power spectra, so that to Langevin simulation data (crosses). The system started from an equilibrium configuration at T = 5.0 GeV. The larger ξ (cf. eq. (2.12)), the longer the system stays close to equilibrium. then the equilibrium form reads To see how fast the system approaches this limit, we can make the approximation of temperatureindependent coefficients in eq. (2.14), obtaining Forη(x 2 − x 1 ) ≪ 1 the terms 3ζ/(2η) cancel, so that non-equilibrium manifests itself by the system staying close to the old value. Forη(x 2 − x 1 ) ≫ 1, the system loses memory of initial conditions and moves towards the equilibrium value p 2 eq = 3ζ/(2η) ∼ m ϕ /T . As the variable x is of O(1), we can say that kinetic decoupling starts whenη < 1.

Time discretization and numerical simulations
We now move on to a numerical integration of eq. (2.10). The problem is technically nonchallenging (unless one is interested in the distribution of momenta in the far UV tail), and we employ a simple-minded approach. The time-like variable x is discretized, and we denote byp i n ,η n andζ n the values of the momenta and coefficients on the corresponding grid. For eq. (2.10) we use the Ito discretization with Gaussian noise, Here, the thermodynamic functions appearing inη n andζ n are interpolated from the tabulated values given in [20] by the cubic-spline method [21].
As the coefficients in eq. (3.1) change by 4 orders of magnitude in the temperature range studied (cf. fig. 1), it is important to have a small enough time step. We have found that this requirement can be sufficiently satisfied with dx = 10 −6 . The initialp i 's are drawn from the equilibrium distribution at T = 5 GeV. The momentum distribution at each x is obtained by histograms produced from N = 10 5 independent runs. The error in each histogram bin is calculated from a jackknife analysis, with a block size of 10 3 .
If the momentum distribution obtained from the simulation is denoted by P, then a useful observable is its ratio to the would-be equilibrium value from eq. (2.16), Snapshots of r as obtained from the simulations are illustrated in fig. 3. For ξ, we consider a number of fixed values (cf. fig. 1), spanning the range that is realistic for the model considered in appendix A. It is clear from fig. 3 that the system rapidly departs from equilibrium if ξ is small, and that it does so by retaining power at small momenta, as a remnant from an earlier time (cf. the discussion around eq. (2.17)).
We find that the simulation results are well represented by the functional form of eq. (2.16), parametrized however by a different coefficient which we denote by α, The corresponding fits are illustrated in fig. 3. These results lead to an easy crosscheck of the accuracy of the solution, as illustrated by the crosses in fig. 2. This implies that the information in fig. 2, originating from eq. (2.14), is sufficient for determining the full momentum distribution. The conclusion that kinetic non-equilibrium is captured by the value p 2 is not new but appears frequently in the literature, however we have arrived at it as a result of a systematic computation, rather than adopted it as a starting point.

Chemical equilibration with non-equilibrium momentum distribution
Having determined the non-equilibrium momentum distribution in sec. 3, the next step is to implement it in the equation governing the dark matter number density. For the scalar singlet model in the resonant regime, it was demonstrated in ref. [18] that inelastic processes are to a very good approximation described by the leading-order ϕϕ ↔ h reaction, where h stands for the Standard Model Higgs boson, set on-shell. Noting furthermore that freeze-out physics takes place deep in the non-relativistic regime, where πT ≪ m ϕ , we may employ the Boltzmann form for the equilibrium distribution function. We denote the non-equilibrium phase space distribution by f ϕ and the equilibrium one byf ϕ ≡ exp(−ǫ ϕ /T ). The Boltzmann equation for f ϕ then takes the form 2π) 3 , and unspecified notation is explained in the context of eq. (A.1).
To proceed, we integrate eq. (4.1) over p 1 . The number density is denoted by n ϕ ≡ p i f ϕ i , and the left-hand side becomes (∂ t + 3H)n ϕ after partial integration. The system is closed by noting that on the right-hand side, we may parametrize the non-equilibrium distribution function with the information obtained in sec. 3, as with the constraints that where K 2 denotes a modified Bessel function. The constraint on r follows from the definition in eq. (3.2). Going over to the variables Y ϕ ≡ n ϕ /s,Ȳ ϕ ≡n ϕ /s, and again replacing time through x from eq. (2.6), we can rewrite the evolution equation as Here the first annihilation cross section is averaged over the non-equilibrium momentum distribution, where the integration bounds can be established as The second term in eq. (4.3) contains an average with respect to equilibrium distributions, For a numerical illustration, we have fixed m ϕ ≃ 60 GeV and κ ≃ 0.00064, which would yield the correct dark matter abundance in kinetic equilibrium according to ref. [18]. However we vary ξ (cf. eq. (2.12)), in order to obtain an ensemble of non-equilibrium momentum distributions (cf. fig. 3). The corresponding σv rel / σv rel is shown in fig. 4(left). Inserting into eq. (4.3), we obtain the dark matter yield, Y ϕ , as illustrated in fig. 4(right). As the system falls out of kinetic equilibrium, the annihilation cross section is reduced, and consequently freeze-out takes place earlier, leading to a larger dark matter abundance. To keep the dark matter abundance at the correct value, the coupling would need to be correspondingly increased, however for our benchmark the effect is only on the 20 percent level.

Conclusions
The goal of this paper has been to demonstrate that Langevin simulations are well suited to studying the efficiency of kinetic equilibration of non-relativistic dark matter candidates, produced through the freeze-out mechanism. The advantage of the Langevin framework, compared with the more standard Boltzmann one, is that it cleanly factorizes the slow nonequilibrium problem from the effect of fast reactions. The role of the latter is to determine For small ξ, the non-equilibrium momentum distribution is red-tilted (cf. fig. 3), having then less weight in the domain that contributes to the annihilation cross section. Right: the corresponding Y ϕ , obtained from eq. (4.3). Because of the smaller annihilation cross section, the overall yield freezes out earlier, and thus to a larger value. Between ξ = 10 −7 and ξ = 10 −9 , there is a ∼ 45% difference in the final yield, which could be compensated for by a ∼ 20% change of κ.
the values of the matching coefficients in the Langevin description. The computation of the matching coefficients can be viewed as a quantum field theoretic problem, and therefore pursued up to higher orders of perturbation theory, or in principle even non-perturbatively, as could be relevant for strong interactions at temperatures of a few GeV (cf. the discussion around eq. (A.14)). However, we have remained at the leading order in the current study, in order to conform with the accuracy of literature studies making use of Boltzmann equations.
As far as the non-equilibrium problem goes, our numerical simulations confirm that the momentum distribution retains the Gaussian form despite the rapidly evolving matching coefficients (that said, the resolution of our setup is not sufficient for studying momenta in the far UV tail of the distribution; for this, more advanced techniques would be required). Therefore, for practical purposes, it is enough to know the width of the momentum distribution, given by the quadratic expectation value p 2 (cf. eq. (2.14)).
In order to illustrate these general points, we chose the example of the scalar singlet model, in a mass regime where an efficient s-channel resonance drives the freeze-out dynamics down to low temperatures. This is among the main examples for which the viability of the kinetic equilibrium assumption has been questioned. If the same processes are included in the computation of the momentum diffusion coefficient (cf. fig. 1(left)) as in the respective literature [5][6][7][8][9], our final phenomenological conclusion turns out to be similar. In particular, if all processes are included, then ξ > 10 −9 in the domain T > 1 GeV in which the freeze-out dynamics take place (cf. fig. 1(left)). Then kinetic non-equilibrium has a < 45% influence, as shown by a comparison of the ξ = 10 −9 and ξ = 10 −7 curves in fig. 4(right), the latter of which represents practically the equilibrium solution. In terms of the coupling κ (cf. eq. (A.1)), this corresponds to a < 20% effect.

A. Momentum diffusion coefficient in the scalar singlet model
We illustrate the computation of the momentum diffusion coefficient ζ of eq. (2.4) with the help of the scalar singlet model, defined by the Lagrangian where φ is the Higgs doublet. The computation can be carried out in a local Minkowskian frame. In the Higgs phase, where φ ≃ (0 v) T / √ 2, with v ≃ 246 GeV, the scalar singlet mass is denoted by m ϕ , and has the value m 2 ϕ ≃ m 2 ϕ0 + κv 2 /2. We consider freeze-out dynamics, taking place in the regime πT ≪ m ϕ ≪ v. Then ϕ can be represented by an effective non-relativistic field ψ, defined as Inserting this in eq. (A.1), integrating over fast oscillations, and anticipating the presence of a conserved charge (cf. eq. (A.4)), whereby we introduce a chemical potential µ ≫ T guaranteeing a small overall number density, we find that the dynamics of ψ is described by where h is an off-shell mode of the Higgs field, with energy e h ≪ m ϕ . Relative to the structures shown, the terms omitted are suppressed by powers of ∂ 0 /m ϕ , h/v, or λ ϕ . Now, the low-energy description of eq. (A.3) has an emergent U(1) symmetry, corresponding to conserved particle number. The Noether current reads where the higher-order terms are suppressed by ∇ 2 /m 2 ϕ . What is important for us is the "force" acting on the dark matter particles. By making use of equations of motion and integrating over the force density, this can be expressed as where the terms omitted are of the same type as in eqs. (A.3) and (A.4).
In accordance with eq. (2.4), the momentum diffusion coefficient ζ is given by the autocorrelator of the force. In quantum field theory, it is advantageous to define the autocorrelator as the zero-frequency limit of a time-symmetrized expectation value, In addition, correlation functions should be normalized so that the overall density drops out, which can be done with the help of the conserved Noether charge, where the time can be taken to be Euclidean. Thereby the momentum diffusion coefficient can be obtained as [17] where the subscript in F i x denotes the spatial position, and the definition is from eq. (A.5). In order to evaluate eq. (A.8), the first step is to insert eq. (A.5), go over to a path integral representation, and carry out the contractions over the fields ψ and ψ * . The propagators are non-relativistic (i.e. with poles only in one half-plane). The contractions represent the thermally averaged amplitude squared for a process in which a heavy ψ interacts with an off-shell h. Let k be the momentum transfer from h, so that q = p + k, where p and q are the momenta of ψ before and after the interaction. Estimating p ∼ m ϕ T and k ∼ T (see below), and employing non-relativistic energies ǫ p = p 2 /(2m ϕ ) and ǫ q = q 2 /(2m ϕ ), the Boltzmann weight does not depend on k to leading order in T /m ϕ , Therefore the thermal average over k effectively localizes h x , Furthermore the overall density of the charge matter particles (determined by the chemical potential µ) cancels between the numerator and denominator. Thus we are left with The next step is to consider various interactions experienced by the Higgs boson, where only operators linear in h need to be included at leading order. Going to momentum space and contracting over the (off-shell) Higgs boson yields where k ≡ |k|, k = d 3 k (2π) 3 , and we went to the static limit as required by the definition of Γ. A further simplification follows by noting that the momentum integral is saturated by k ∼ πT ≪ m h , cf. eqs. (A. 16) and (A.19). Therefore the momentum diffusion coefficient can be approximated as In other words, the Higgs exchange is a contact interaction at low energies. We stress that the operator O x is gauge invariant under QCD, and that therefore eq. (A.14) is defined and computable beyond perturbation theory.
To proceed with a leading-order evaluation of eq. (A.14), consider first a fermionic operator containing the b quark, In the numerical illustrations in fig. 1, we also include the c quark. A few-page thermal field theory computation yields where n F is the Fermi distribution. Changing variables, we finally get where the upper bound is saturated in the limit m b ≪ πT . Numerically, the upper bound gives a fairly good approximation for charm quarks.
As a second example, we consider the bosonic operator obtained by integrating out the top quark [22], whose relevance for thermal considerations has been underlined in ref. [23]. In this case we obtain where n B is the Bose distribution. Changing variables and carrying out the integral yields In order to illustrate the magnitude of these corrections, we show in fig. 1(left) how they can be converted into the coupling ξ, defined in eq. (2.12).