Coannihilation without chemical equilibrium

TL;DR

This paper introduces a new conversion-driven freeze-out mechanism for dark matter, where chemical equilibrium fails, leading to relic densities governed by conversion rates rather than self-annihilation, with distinctive collider signatures.

Contribution

It demonstrates a novel freeze-out process without chemical equilibrium, expanding the understanding of dark matter relic density determination beyond traditional assumptions.

Findings

Identifies parameter regions with small couplings matching observed relic density.

Shows freeze-out occurs out of chemical equilibrium with inefficient self-annihilation.

Predicts collider signatures like disappearing tracks or displaced vertices.

Abstract

Chemical equilibrium is a commonly made assumption in the freeze-out calculation of coannihilating dark matter. We explore the possible failure of this assumption and find a new conversion-driven freeze-out mechanism. Considering a representative simplified model inspired by supersymmetry with a neutralino- and sbottom-like particle we find regions in parameter space with very small couplings accommodating the measured relic density. In this region freeze-out takes place out of chemical equilibrium and dark matter self-annihilation is thoroughly inefficient. The relic density is governed primarily by the size of the conversion terms in the Boltzmann equations. Due to the small dark matter coupling the parameter region is immune to direct detection but predicts an interesting signature of disappearing tracks or displaced vertices at the LHC. Unlike freeze-in or superWIMP scenarios, conversion-driven freeze-out is not sensitive to the initial conditions at the end of reheating.

Figures

Tables

Table 1. Table 1: List of all included coannihilation processes. We use the abbreviations q = 𝑞 absent q= all quarks, l = 𝑙 absent l= all leptons and V = g , γ , Z , W 𝑉 𝑔 𝛾 𝑍 𝑊 V=g,\gamma,Z,W . The b ~ ~ 𝑏 \widetilde{b} annihilation into b ​ b ¯ 𝑏 ¯ 𝑏 b\bar{b} also has contributions scaling with λ χ 2 superscript subscript 𝜆 𝜒 2 \lambda_{\chi}^{2} and λ χ 4 superscript subscript 𝜆 𝜒 4 \lambda_{\chi}^{4} .

initial state		final state		scaling
$\chi$	$\chi$	$b$	$\overline{b}$	${\lambda }_{\chi }^4$
$\chi$	$\tilde{b}$	$b$	$g,\gamma ,Z,H$	${\lambda }_{\chi }^2$
		$W^{-}$	$t,u,c$
$\tilde{b}$	${\tilde{b}}^{†}$	$V$	$V$	${\lambda }_{\chi }^0$
		$q$	$\overline{q}$
		$Z$	$H$
		$l$	$\overline{l}$
$\tilde{b}$	$\tilde{b}$	$b$	$b$	${\lambda }_{\chi }^4$

Table 2. Table 2: List of all considered conversion processes. Processes marked with ∗ have soft divergences, processes with ∗∗ can have t 𝑡 t -channel divergences.

initial state		final state		symbol	scaling
$\chi$	$b$	$\tilde{b}$	$g^{\ast },{\gamma }^{\ast },Z,H$	${\Gamma }_{\chi \to \tilde{b}}$	${\lambda }_{\chi }^2$
	$g,\gamma ,Z^{\ast \ast },H^{\ast \ast }$		$\overline{b}$
	$W^{-}$		$\overline{t},{\overline{u}}^{\ast \ast },{\overline{c}}^{\ast \ast }$
	$t^{\ast \ast },u,c$		$W^{+}$
$\tilde{b}$		$\chi$	$b$	${\Gamma }_{\tilde{b}}$	${\lambda }_{\chi }^2$
$\chi$	$\chi$	$\tilde{b}$	${\tilde{b}}^{†}$	$⟨{\sigma }_{\chi \chi \to \tilde{b}{\tilde{b}}^{†}}v⟩$	${\lambda }_{\chi }^4$

Equations

$${\mathcal{L}}_{\text{int}}=|D_{\mu }\tilde{q}{|}^2-{\lambda }_{\chi }\tilde{q}\bar{q}\frac{1-{\gamma }_5}{2}\chi +\text{h.c.},$$

$${\mathcal{L}}_h={\lambda }_h{h}^{†}h{\tilde{b}}^{†}\tilde{b}$$

$$\frac{n_{\chi }}{n_{\chi }^{\mathrm{eq}}}=\frac{n_{\tilde{b}}}{n_{\tilde{b}}^{\mathrm{eq}}},$$

$$\frac{\mathrm{d}{Y}_{\chi }}{\mathrm{d}x}=\frac{1}{3H}\frac{\text{mbox}ds}{\text{mbox}dx}$$

$$\frac{\mathrm{d}{Y}_{\tilde{b}}}{\mathrm{d}x}=\frac{1}{3H}\frac{\text{mbox}ds}{\text{mbox}dx}$$

$$\Gamma_{\chi\rightarrow\widetilde{b}}=2\sum_{k,l}\big{\langle}\sigma_{\chi k\rightarrow\widetilde{b}l}v\big{\rangle}\,n_{k}^{\mathrm{eq}}\,,$$

$$\begin{array}{rl}& ⟨{\sigma }_{ij}v⟩n_i^{\mathrm{eq}}{n}_j^{\mathrm{eq}}=\\ & T\frac{g_i{g}_j}{256{\pi }^5}\int \frac{p_{ij}{p}_{ab}}{\sqrt{s}}|\overline{M}{|}^2{K}_1\left(\frac{\sqrt{s}}{T}\right)\mathrm{d}s\mathrm{d}\cos \theta ,\end{array}$$

$${\Gamma }_{\tilde{b}}\equiv \Gamma ⟨\frac{1}{\gamma }⟩=\Gamma \frac{K_1\left(m_{\tilde{b}}/T\right)}{K_2\left(m_{\tilde{b}}/T\right)},$$

$${\Gamma }_{\tilde{b}}\frac{Y_{\tilde{b}}^{\mathrm{eq}}}{Y_{\chi }^{\mathrm{eq}}}\sim H,$$

$${\mathcal{F}}_{\text{pass}}^i=e^{-\ell /\left(c\tau \beta \gamma \right)},$$

$${\overline{\mathcal{F}}}_{\text{pass}}=\frac{{\sum }_i{\mathcal{F}}_{\text{pass}}^i{\mathcal{P}}_{\text{on}}^i{\mathcal{P}}_{\text{off}}^i}{{\sum }_i{\mathcal{P}}_{\text{on}}^i{\mathcal{P}}_{\text{off}}^i},$$

$$V(r)\approx \frac{{\alpha }_s}{2r}\left[C_Q-C_R-C_{R^{\prime }}\right]$$

$$S_0=-\frac{\pi \alpha /\beta }{1-e^{\pi \alpha /\beta }},$$

$$S_{0}\rightarrow\frac{2}{7}S_{0}\Bigg{|}_{\alpha=-4/3\alpha_{s}}+\frac{5}{7}S_{0}\Bigg{|}_{\alpha=1/6\alpha_{s}}\;.$$

$$f_X(p,t)=f_X^{\mathrm{eq}}(p,T)\frac{Y_X(t)}{Y_X^{\mathrm{eq}}(T)},$$

$$\left({\partial }_t-Hp{\partial }_p\right)f_{\chi }(p,t)=\frac{1}{E_{\chi }}C\left[f_{\chi }\right],$$

$$x(t,p)=m_{\chi }/T,q(t,p)=p/T,$$

$$\begin{array}{rl}\left({\partial }_t-Hp{\partial }_p\right)& =\frac{H}{1-\frac{x}{3h_{\text{eff}}}\frac{\mathrm{d}{h}_{\text{eff}}}{\mathrm{d}x}}\left(x{\partial }_x+\frac{x}{3h_{\text{eff}}}\frac{\mathrm{d}{h}_{\text{eff}}}{\mathrm{d}x}q{\partial }_q\right)\\ & \approx Hx{\partial }_x.\end{array}$$

$$C_{12}[f_{\chi }]=\frac{1}{2}\int \mathrm{d}{\Pi }_b\mathrm{d}{\Pi }_{\tilde{b}}(2\pi {)}^4{\delta }^4\left(\sum p_i\right)|\overline{M}{|}^2\left[f_{\tilde{b}}-f_{\chi }{f}_b\right],$$

$$\begin{array}{rl}{C}_{22}[f_{\chi }]=& \frac{1}{2}\int \mathrm{d}{\Pi }_A\mathrm{d}{\Pi }_{\tilde{b}}\mathrm{d}{\Pi }_B(2\pi {)}^4{\delta }^4\left(\sum p_i\right)|\overline{M}{|}^2\\ & \times \left[f_{\tilde{b}}{f}_B-f_{\chi }{f}_A\right],\end{array}$$

$$f_{\tilde{b}}-f_{\chi }{f}_b=f_{\tilde{b}}^{\mathrm{eq}}\frac{Y_{\tilde{b}}}{Y_{\tilde{b}}^{\mathrm{eq}}}-f_{\chi }{f}_b^{\mathrm{eq}}=f_b^{\mathrm{eq}}\left(f_{\chi }^{\mathrm{eq}}\frac{Y_{\tilde{b}}}{Y_{\tilde{b}}^{\mathrm{eq}}}-f_{\chi }\right).$$

$$C[f_{\chi }]=\tilde{C}(q,x)\left(f_{\chi }^{\mathrm{eq}}\frac{Y_{\tilde{b}}}{Y_{\tilde{b}}^{\mathrm{eq}}}-f_{\chi }\right)E_{\chi }.$$

$$Hx{\partial }_x{f}_{\chi }(q,x)=\tilde{C}(q,x)\left(f_{\chi }^{\mathrm{eq}}\frac{Y_{\tilde{b}}}{Y_{\tilde{b}}^{\mathrm{eq}}}-f_{\chi }\right).$$

$$f_{\chi }(q,x_0)=f_{\chi }^{\mathrm{eq}}(q,x_0)=\exp \left(-\sqrt{q^2+x_0^2}\right)$$

$$\begin{array}{rl}{f}_{\chi }(q,x)=& f_{\chi }^{\mathrm{eq}}(q,x)\frac{Y_{\tilde{b}}}{Y_{\tilde{b}}^{\mathrm{eq}}}-{\int }_{x_0}^x\frac{\mathrm{d}\left(f_{\chi }^{\mathrm{eq}}(q,y)Y_{\tilde{b}}(y)/Y_{\tilde{b}}^{\mathrm{eq}}(y)\right)}{\mathrm{d}y}\\ & \times \exp \left(-{\int }_y^x\frac{\tilde{C}(q,z)}{zH(z)}\mathrm{d}z\right)\mathrm{d}y.\end{array}$$

$$n_{\chi }(x)=\frac{4\pi g_{\chi }{m}_{\chi }^3}{x^3}\int \frac{q^2}{(2\pi {)}^3}{f}_{\chi }(q,x)\mathrm{d}q,$$

$${\tilde{C}}_{12}=\frac{2g_{\tilde{b}}{g}_b|M{|}^{2}T}{16\pi p_{\chi }{E}_{\chi }}\left(e^{-p_{\min }/T}-e^{-p_{\max }/T}\right),$$

$$|M{|}^2$$

$$p_{\min /\max }$$

$$\begin{array}{rl}{\tilde{C}}_{22}=& \frac{2g_{q}T}{16{\pi }^2{p}_{\chi }{E}_{\chi }}{\int }_{m_{\tilde{b}}^2}^{\infty }(s-m_{\chi }^2)\sigma (s)\\ & \times \left(e^{-p_{\min }(s)/T}-e^{-p_{\max }(s)/T}\right)\mathrm{d}s\end{array}$$

$$\sigma (s)=\frac{1}{4|p_{\chi }\cdot p_q|}\int \mathrm{d}{\Pi }_{\tilde{b}}\mathrm{d}{\Pi }_g(2\pi {)}^4\delta (p_{\chi }+p_q-p_{\tilde{b}}-q_g)|\mathcal{M}{|}^2$$

Full text

††preprint: TUM-HEP 1085/17††preprint: TTK-17-18

I Introduction

The origin and the nature of dark matter (DM) in the Universe is one of the most pressing questions in particle- and astrophysics. Despite impressive efforts to uncover its interactions with the Standard Model (SM) of particle physics in (in)direct detection and accelerator-based experiments, DM remains elusive and, so far, our understanding is essentially limited to its gravitational interactions (see e.g. Bertone et al. (2005); Arcadi et al. (2017)). It is therefore of utmost interest to investigate mechanisms for the generation of DM in the early Universe that go beyond the widely studied paradigm of thermal freeze-out, and that can point towards non-standard signatures.

In this spirit we subject the well-known coannihilation scenario Griest and Seckel (1991) to further scrutiny and investigate the importance of the commonly made assumption of chemical equilibrium (CE) between DM and the coannihilation partner. This requires solving the full set of coupled Boltzmann equations which has been performed in the context of specific supersymmetric scenarios Chung et al. (1997); Ellis et al. (2015). Here we consider a simplified DM model and explore the break-down of CE in detail finding a new, conversion-driven solution for DM freeze-out which points towards a small interaction strength of the DM particle with the SM bath. While the smallness of the coupling renders most of the conventional signatures of DM unobservable, new opportunities for collider searches arise. In particular we find that searches for long-lived particles at the LHC are very powerful tools for testing conversion-driven freeze-out.

The structure of the paper is as follows: We begin by introducing a simplified model for coannihilations in Sec. II. In Sec. III we present the Boltzmann equations which govern the DM freeze-out and investigate conversion-driven solutions before we confront the regions of parameter space which allow for a successful generation of DM with LHC searches in Sec. IV. Finally, we summarize our results and conclude.

In the appendix we provide further details and results. In particular, we discuss the appearance of divergencies in the conversion rates in Appendix A, describe our treatment of Sommerfeld enhancement in Appendix B and justify the assumption of kinetic equilibrium by providing solutions of the full momentum-dependent Boltzmann equation in Appendix C.

II Simplified model for coannihilation

While the precise impact of the breakdown of CE between DM and its coannihilation partner will in general depend on the details of the considered model, the key aspects of the phenomenology can be expected to be universal. As a representative case we choose a simplified model for DM interacting with quarks. We extend the matter content of the SM minimally by a Majorana fermion $\chi$, being a singlet under the SM gauge group, and a scalar quark-partner $\widetilde{q}$, mediating the interactions with the SM and acting as the coannihilation partner. The interactions of the new particles among themselves and with the SM are given by Garny et al. (2012)

[TABLE]

where $q$ is a SM quark field, $D_{\mu}$ denotes the covariant derivative, which contains the interactions of $\widetilde{q}$ with the gauge bosons as determined by its quantum numbers, and $\lambda_{\chi}$ is a Yukawa coupling. Here we choose $q=b$ and $Y=-\frac{1}{3}$. For the coupling $\lambda_{\chi}=\frac{1}{3}\sqrt{2}\frac{e}{\cos\theta_{W}}\approx 0.17$ our simplified model makes contact with the Minimal Supersymmetric SM (MSSM) where $\widetilde{b}$ can be identified with a right-handed sbottom and $\chi$ with a bino-like neutralino. However, here we will not consider the MSSM but vary $\lambda_{\chi}$ freely.111For a realization of small $\lambda_{\chi}$ in extensions of the MSSM see Ibarra et al. (2009). Note that choosing a top-partner mediator instead yields similar results although quantitative differences arise due to the large top mass.

On top of the gauge and Yukawa interactions described above a Higgs-portal interaction given by

[TABLE]

is also allowed. This interaction does not involve $\chi$ directly and has no impact on the conversion rates $\chi\leftrightarrow\widetilde{b}$, that are responsible for establishing CE. Nevertheless, it can modify the annihilation rate of $\widetilde{b}$. Since the additional contributions involving the scalar coupling compete with QCD processes, i.e. $\widetilde{b}\widetilde{b}^{\dagger}\rightarrow gg$, they are sub-leading unless $\lambda_{h}$ is very large. Even in this case we do not expect qualitative differences, and therefore neglect this contribution in the following.

III Freeze-out without chemical equilibrium

For coannihilation to be effective the coannihilating particles – here $\chi$ and $\widetilde{b}$ – have to be in thermal contact through efficient conversion rates $\chi\leftrightarrow\widetilde{b}$. For couplings $\lambda_{\chi}$ of the order of the electroweak coupling strength, conversion rates are much larger than the Hubble rate $H$ during freeze-out, guaranteeing CE, i.e.

[TABLE]

where $n$ ($n^{\mathrm{eq}}$) is the (equilibrium) number density. While CE holds the results are not sensitive to the size of the conversion rates and the Boltzmann equations can be reduced to a single equation that does not contain conversion terms Edsjö and Gondolo (1997). This approach is solved in standard tools Belanger et al. (2006); Gondolo et al. (2004); Backović et al. (2014).

For smaller couplings, however, CE can break down and the full coupled set of Boltzmann equations has to be solved Chung et al. (1997); Ellis et al. (2015). In our case:

[TABLE]

where $Y=n/s$ is the comoving number density, $s$ the entropy density and $x=m_{\chi}/T$. We take the internal degrees of freedom $g_{\chi}=2$ and $g_{\widetilde{b}}=3$. $Y_{\widetilde{b}}$ represents the summed contribution of the mediator and its anti-particle. Since the cross sections are averaged over initial state degrees of freedom, this leads to the factor $1/2$ in the respective Boltzmann equation, Eq. (5). Equally, $\Gamma_{\chi\rightarrow\widetilde{b}}$ is understood to contain the conversion into both.

Apart from the familiar annihilation and coannihilation terms displayed in the first lines of Eqs. (4) and (5) three additional rates for the conversion processes enter in the second lines. The first term includes all the scattering processes which convert DM in its coannihilation partner. The scattering rate is given by

[TABLE]

where $k,l$ denote SM particles. The factor of two arises from annihilation into the mediator and its anti-particle. Neglecting quantum statistical factors and assuming Boltzmann distributions the thermal average for scattering reads Edsjö and Gondolo (1997)

[TABLE]

where $K_{i}$ denotes a modified Bessel function of the second kind. Here $g_{i}$ are the internal degrees of freedom of species $i$, $p_{ij}$ and $p_{ab}$ denote the absolute value of the three momentum of the initial and final state particles in the center-of-mass frame, respectively. In certain scattering processes soft- or $t$-channel divergences can appear in the thermal averages. We discuss these issues and how we resolve them in detail in Appendix A. The next term in the second lines of Eqs. (4) and (5) captures the conversion induced by the decay and inverse decay of $\widetilde{b}$. This rate is controlled by the thermally averaged decay width $\Gamma_{\widetilde{b}}$. The thermally averaged decay rate is given by

[TABLE]

where $\gamma$ is the Lorentz factor. Finally, the last term takes the scattering processes in the odd-sector into account. We include all diagrams that are allowed at tree-level and use FeynRules Alloul et al. (2014) and CalcHEP Belyaev et al. (2013) to generate the squared matrix elements $|\overline{M}|^{2}$, see Tables 1 and 2 for a summary of the included processes. We take into account the Sommerfeld enhancement of the $\widetilde{b}\widetilde{b}^{\dagger}$ annihilation rates as detailed in Appendix B.

A first naive estimate that allows us to determine the ballpark of the parameters where conversion processes become relevant follows from demanding

[TABLE]

for temperatures relevant to the freeze-out dynamics ($m_{\chi}/T\sim 30$). Using as a representative benchmark $m_{\chi}=500\,$GeV and a value of the mediator mass that allows for coannihilations ($m_{\widetilde{b}}=510$ GeV) this relation indicates $\lambda_{\chi}\sim\mathcal{O}(10^{-7})$. The order of magnitude is largely insensitive to the precise choice of masses, as long as coannihilations can occur.

For such a small coupling a clear hierarchy emerges between the different rates, see left panel of Fig. 1. The annihilation $\chi\chi\rightarrow$ SM SM and $\chi\chi\rightarrow\widetilde{b}\widetilde{b}^{\dagger}$ that are proportional to $\lambda_{\chi}^{4}$ and thermally suppressed by $n_{\chi}^{\mathrm{eq}}$ are exceedingly small and cannot compete with the Hubble expansion. Even though the coannihilation rate $\chi\widetilde{b}\rightarrow$ SM SM, which scales as $\sigma_{\chi\widetilde{b}}v\propto\lambda_{\chi}^{2}g^{2}$ (where $g$ is a SM gauge coupling) is enhanced relative to this by many orders of magnitude it is also negligible compared to $H$. In contrast, the leading contribution to $\widetilde{b}\widetilde{b}^{\dagger}\rightarrow$ SM SM is set by the gauge interactions of $\widetilde{b}$ and, therefore, the rate remains comfortably larger than $H$ until $T\approx m_{\chi}/30$. The most important annihilations, especially for very small $\lambda_{\chi}$ are the $\widetilde{b}$ annihilations into gluons. Since the interaction rates are suppressed exponentially by the masses of external particles, it is clear that the conversion processes containing external gluons dominate over the rates containing weak scale particles. The conversion rates are close to the Hubble rate and, for this choice of couplings, just about sufficient to make conversion processes relevant for the freeze-out.

Taking these rates and solving the Boltzmann equations we find the results presented in the right-hand side of Fig. 1. We solve the system of coupled equations from $x=1$ up to $x=1000$.222Due to efficient annihilations, the $\widetilde{b}$ abundance is very close to equilibrium at early times. For numerical convenience, it is sufficient to track its deviation from equilibrium starting from $x\sim 15$. The $\chi$ abundance leaves its equilibrium value already at rather high temperatures, well before the freeze-out of a typical thermal relic or the $\widetilde{b}$ freeze-out. The slow decline of the $\chi$ abundance after this point is due to the close-to-inefficient conversion terms which remove over-abundant $\chi$s.

In Fig. 2 we show the dependence of the final freeze-out density on the coupling $\lambda_{\chi}$ (red solid line). For large enough coupling, the solution coincides with the result that would be obtained when assuming CE (blue dotted line). The relic density is in this case largely set by the strength of $\widetilde{b}$ self-annihilation into gluons. When lowering the value of $\lambda_{\chi}$, conversions $\chi\leftrightarrow\widetilde{b}$ become less efficient and one obtains a relic density that lies above the value expected for CE. For the benchmark scenario shown in Fig. 2, the freeze-out density matches the value determined by Planck Ade et al. (2016) for a coupling of $\lambda_{\chi}\approx 2.6\times 10^{-7}$.

Above we assumed that both $\chi$ and $\widetilde{b}$ have thermal abundances for $T\gg m_{\chi}$. While this assumption is certainly well justified for $\widetilde{b}$, one may question the dependence on the initial condition for $\chi$ due to its small coupling to the thermal bath. We check the dependence on this assumption by varying the initial abundance at $T=m_{\chi}$ within the range $(0\!-\!100)\times Y_{\chi}^{\mathrm{eq}}$. The evolution of the abundances for our benchmark point are shown in Fig. 3, for early times ($x<20$). We find that all trajectories converge before $x\lesssim 5$, thereby effectively removing any dependence of the final density on the initial condition at $x=1$. The dependence of the final freeze-out density on the initial condition is also indicated in Fig. 2 by the area shaded in red, and is remarkably small. Therefore, conversion-driven freeze-out is largely insensitive to details of the thermal history prior to freeze-out and in particular to a potential production during the reheating process. Note that this feature distinguishes conversion-driven freeze-out from scenarios for which DM has an even weaker coupling such that it was never in thermal contact (e.g. freeze-in production Hall et al. (2010)). Thus, while requiring a rather weak coupling, the robustness of the conventional freeze-out paradigm is preserved in the scenario considered here.

As discussed before, conversions $\chi\leftrightarrow\widetilde{b}$ are driven by two types of processes, decay and scattering. It turns out that quantitatively both are important for determining the freeze-out density. To illustrate the importance of scattering processes, we show the freeze-out density that would be obtained when only taking decays into account by the gray dashed line in Fig. 2. Furthermore, the gray shaded area indicates the dependence on initial conditions that would result neglecting scatterings. We find that scattering processes, which are active at small $x$, are responsible for wiping out the dependence on the initial abundance in the full solution of the coupled Boltzmann equations.

While Eqs. (4) and (5) do not require chemical equilibrium to hold, they rely on the implicit assumption of kinetic equilibrium, i.e. the assumption that the momentum dependence of the distribution functions is proportional to a thermal distribution. For the freeze-out of weakly interacting massive particles, this assumption is typically well justified because scattering processes are enhanced compared to annihilations by a Boltzmann factor of order $e^{m_{\chi}/T}$. However, in the present scenario, the small coupling $\lambda_{\chi}$ renders elastic scatterings of the form $\chi b\leftrightarrow\chi b$ inefficient. Instead, the leading processes to establish kinetic equilibrium are the inelastic conversion processes discussed above. Since their rate is, by definition, comparable to the Hubble rate in the interesting regime of parameters, one may wonder whether the treatment based on integrated Boltzmann equations is justified for $\chi$. In order to check this point, we solved the full, momentum-dependent Boltzmann equation for $\chi$, taking the leading decay and scattering processes into account. We find that, while the distribution function can indeed deviate from the thermal distribution at intermediate times, the final relic abundance differs only mildly from the integrated treatment (below the 10% level). The main reason is that the collision operator does not depend strongly on the momentum mode, such that all modes behave in a similar way. For a detailed discussion we refer to appendix C.

Let us briefly comment on possible refinements. Apart from quantum statistics, also thermal effects could play a role at small $x$. In particular, the thermal mass for the $b$-quark can lead to a thermal blocking of the decay at high temperatures and for very small mass splitting. Since a consistent inclusion of this effect would require to take also further thermal processes into account, and since (hard) scatterings dominate for small $x$, we do not expect these corrections to significantly affect our conclusion. Additionally, bound state effects could play a role for the $\widetilde{b}$ annihilation Kim and Laine (2016, 2017); Liew and Luo (2017); Mitridate et al. (2017); Keung et al. (2017).

IV Viable parameter space

We will now explore the parameter space consistent with a relic density that matches the DM density measured by Planck, $\Omega h^{2}=0.1198\pm 0.0015$ Ade et al. (2016). In the considered scenario, for small couplings, $\widetilde{b}\widetilde{b}^{\dagger}$ annihilation is the only efficient annihilation channel. Hence the minimal relic density that can be obtained for a certain point in the $m_{\chi}$-$m_{\widetilde{b}}$ plane is the one for a coupling $\lambda_{\chi}$ that just provides CE (but is still small enough so that $\chi\chi$- and $\chi\widetilde{b}$-annihilation is negligible). The curve for which this choice provides the right relic density defines the boundary of the valid parameter space and is shown as a black, solid curve in Fig. 4. Below this curve a choice of $\lambda_{\chi}$ sufficiently large to support CE would undershoot the relic density. In this region a solution with small $\lambda_{\chi}$ exists that renders the involved conversion rates just large enough to allow for the right portion of thermal contact between $\widetilde{b}$ and $\chi$ to provide the right relic density. The value of $\lambda_{\chi}$ ranges from $10^{-7}$ to $10^{-6}$ (from small to large $m_{\chi}$). These values lie far beyond the sensitivity of direct or indirect detection experiments.

For the solutions providing the right relic density, during typical freeze-out (i.e. when $T\sim m_{\chi}/30$) the conversion rates have to be on the edge of being efficient, cf. Eq. (9). From this simple relation (and assuming that the decay width, $\Gamma_{\widetilde{b}}$, is similar in size as the other conversion rates) we can already infer that the decay length of $\widetilde{b}$ is of the order of 1–100 cm for a DM particle with a mass of a few hundred GeV predicting the signature of disappearing tracks or displaced vertices at the LHC.

The decay length in our model is shown as the gray dotted lines in Fig. 4. It ranges from $25\,\text{cm}$ to below $2.5\,\text{cm}$ for increasing mass difference (the dependence on the absolute mass scale is more moderate).

In proton collisions at the LHC pairs of $\widetilde{b}$s could be copiously produced. They will hadronize to form $R$-hadrons Farrar and Fayet (1978) which will, for the relevant decay length, either decay inside or traverse the sensitive parts of the detector. Accordingly, the signatures of displaced vertices and (disappearing) highly ionizing tracks provide promising discovery channels. Similar searches have, e.g., been performed for a gluino $R$-hadron (decaying into energetic jets) Aad et al. (2015) or a purely electrically charged heavy stable particle Aad et al. (2013); Khachatryan et al. (2015) but have not been performed for the model under consideration (see also Mahbubani et al. (2017); Buchmueller et al. (2017) for a recent account on simplified DM models providing displaced vertices). However, some constraints on the model can already be derived from existing searches.

Searches for detector-stable $R$-hadrons Chatrchyan et al. (2013); Khachatryan et al. (2016); Aad et al. (2015); Aaboud et al. (2016) can be reinterpreted for finite decay lengths by convoluting the signal efficiency with the fraction of $R$-hadrons that decay after traversing the relevant parts of the detector. This reinterpretation provides limits down to a decay length of $c\tau\simeq 0.1\,$m for an $R$-hadron mass around 100 GeV and can be used to estimate excluded parameter regions in our model.

To this end we compute the weighted fraction of $R$-hadrons that decay after traversing the relevant parts of the detector in a Monte Carlo simulation as follows. For a given $R$-hadron in an event $i$ this fraction is

[TABLE]

where $\ell=\ell(\eta)$ is the travel distance to pass the respective part of the detector which depends on the pseudo-rapidity $\eta$ while $\gamma$ is the Lorentz factor according to the velocity $\beta$. We use a simple cylindrical approximation for the CMS tracker333We considered the tracker-only and tracker+muon-system analysis of Chatrchyan et al. (2013) finding the higher sensitivity for the former one. with a radius and length of 1.1 m and 5.6 m, respectively. For the weighting we compute444For simplicity we display the formula for one $R$-hadron candidate per event, for events with two candidates we follow the prescription in Khachatryan et al. (2015) (with the replacement $\mathcal{P}^{i}_{\text{off}}\to\mathcal{F}^{i}_{\text{pass}}\mathcal{P}^{i}_{\text{off}}$ in the respective sum in the numerator of Eq. (11)).

[TABLE]

where $\mathcal{P}^{i}_{\text{on}}$ and $\mathcal{P}^{i}_{\text{off}}$ are the probabilities of the respective event to be triggered and pass the selection cuts, respectively, and the sum runs over all generated events. We use the tabulated probabilities  $\mathcal{P}^{i}_{\text{on}},\mathcal{P}^{i}_{\text{off}}$ for lepton-like heavy stable charged particles following the prescription in Khachatryan et al. (2015) (see also Heisig et al. (2015) for details of the implementation of isolation criteria and validation). We expect this to be a good approximation as the selection criteria for lepton-like heavy stable charged particles and $R$-hadrons are identical and differences in the overall detector efficiency cancel out in Eq. (11). We simulate events with MadGraph5_aMC@NLO Alwall et al. (2014), performing showering and hadronization with Pythia 6 Sjostrand et al. (2006).

We use the cross section predictions from NLLFast Beenakker et al. (2010) and rescale the signal by $\overline{\mathcal{F}}_{\text{pass}}$. The 95% CL exclusion limits are then obtained from a comparison to the respective cross section limits from searches for (top-squark) $R$-hadrons presented in Chatrchyan et al. (2013). The results are shown in Fig. 5. We show limits for two models regarding the hadronization and interaction of the $R$-hadron with the detector material, the generic model Kraan (2004); Mackeprang and Rizzi (2007) and Regge (charge-suppressed) model de Boer et al. (2008); Mackeprang and Milstead (2010) as the red solid and blue dashed line, respectively.

In addition to the results for the 8 TeV LHC we show results from an analogous reinterpretation of the preliminary results from $12.9\,\text{fb}^{-1}$ of data from the 13 TeV LHC run Khachatryan et al. (2016). Since the tabulated probabilities in Khachatryan et al. (2015) are only provided for 8 TeV we use these also for the analysis of the 13 TeV simulation assuming a similar detector efficiency for $R$-hadrons in both runs.

The fraction of $R$-hadrons passing the tracker is exponentially suppressed for small life-times significantly weakening the respective sensitivity. However, there are two competing factors that nevertheless result in meaningful limits for $c\tau$ smaller than the detector size. On the one hand, for small masses the production cross section rises quickly. On the other hand, for smaller masses a larger fraction of $R$-hadrons is significantly boosted enhancing the travel distance in the detector. However, this (latter) effect does not significantly enhance the sensitivity as the signal efficiency for largely boosted $R$-hadrons decreases rapidly (as tracks become indistinguishable from minimal ionizing tracks for $\beta\to 1$).

Note that the above CMS analysis has been interpreted for $R$-hadrons containing top-squarks. As discussed in Mackeprang and Milstead (2010) the expected energy loss for an $R$-hadron containing sbottoms is smaller. This results in an efficiency around 30–40% smaller relative to the case of the stop and therefore in slightly weaker limits on the sbottom mass, see e.g. Aad et al. (2015). However, we use the above limit taking the result for the Regge model (that provides the weaker limits) as a realistic estimate of the LHC limits for the (sbottom-like) mediator in our model considering the fact that the uncertainties in the hadronization are of similar size as the difference between the sbottom and stop case.

The resulting limits from the 8 TeV Chatrchyan et al. (2013) and 13 TeV Khachatryan et al. (2016) LHC data are superimposed in Fig. 4. For mass splittings below $m_{b}$ (below gray dashed curve) the 2-body decay is not allowed and the $R$-hadrons can be considered detector-stable. Towards large mass splittings (smaller life-times) the limits fall off significantly providing no constraint above $\Delta m_{\chi\widetilde{b}}\simeq 13\,\text{GeV}$.

In addition, a large number of experimental results for a sbottom-neutralino simplified model assuming a prompt sbottom decay exist, see e.g. Sirunyan et al. (2017); Khachatryan et al. (2017); Sirunyan et al. (2017); Aaboud et al. (2016). While most of these searches are not applicable to non-prompt decays, monojet searches, targeting small mass splittings, have been performed that do not rely on the prompt decay of the mediator Khachatryan et al. (2015); Aaboud et al. (2016). We superimpose the (stronger) limit from Aaboud et al. (2016) that uses $3.2\,\text{fb}^{-1}$ of 13 TeV data.

V Conclusion

In this work we have considered the possibility that the common assumption of chemical equilibrium during DM freeze-out does not hold. For definiteness, we have focused on a simplified model with particle content inspired by supersymmetry, comprising a neutral Majorana fermion as the DM candidate and a colored scalar particle that mediates a coupling to bottom quarks. For small mass splitting between the mediator and the DM particle, the freeze-out is dominated by self-annihilation of the mediator. This process can be efficient enough to deplete the DM density below the observed value, thus giving rise to a portion of parameter space in which thermal freeze-out cannot account for all of the DM abundance. In this work we have demonstrated that this conclusion hinges on the assumption of chemical equilibrium, and that the freeze-out process can account for the DM density determined by Planck when relaxing this assumption. This occurs when the DM particle interacts very weakly with both the SM and the mediator, such that conversion processes have to be taken into account explicitly. We find that this opens up new regions in parameter space which lead to characteristic signatures of long-lived particles at collider experiments. $R$-hadron searches performed at the $8$ and $13$ TeV LHC runs already constrain part of the parameter space providing conversion-driven freeze-out. A dedicated search for disappearing $R$-hadron tracks and displaced vertices targeting decay lengths in the range 1–100 cm is expected to probe an even larger portion of the allowed parameter space.

The mechanism discussed here is distinct from the freeze-in scenario Hall et al. (2010), for which the DM particle was never in thermal equilibrium, and which would require a much smaller coupling strength than considered here. On the other hand, it shares some similarities with the superWIMP scenario (see e.g. Feng et al. (2003)), but also differs in various respects. In particular, the relic density is set by the interplay of conversion and annihilation processes during freeze-out, unlike for superWIMPs, where DM is produced from the late decay of a heavier state that undergoes a standard thermal freeze-out. In addition, for the mechanism considered here, the life-time of the coannihilation partner is short enough such that constraints from Big Bang nucleosynthesis are generally avoided, provided that the decay rate gives a sizable contribution to conversion. Unlike both freeze-in and superWIMP scenarios in the considered mechanism the final relic density is insensitive to the initial condition of the abundance at the end of the reheating process.

We expect that conversion-driven freeze-out can be realized generically in DM models featuring strong coannihilations. If the coannihilation partner is not colored but only electrically charged, one may expect signatures related to lepton-like highly ionizing tracks. Finally, it is possible that the efficient self-annihilation of the coannihilation partner is itself driven by a new interaction beyond the SM Baker et al. (2015). In this case the mechanism described here can be relevant even if the coannihilating state is a SM singlet with macroscopic decay length, potentially leading to displaced vertex signatures.

Note added: Reference D’Agnolo et al. (2017), which appeared recently, discusses a related mechanism.

Acknowledgements

We thank Michael Krämer, Björn Sarrazin, Pasquale Serpico and Wolfgang Waltenberger for helpful discussions. We acknowledge support by the German Research Foundation DFG through the research unit “New physics at the LHC”.

Appendix A Divergences in conversion rates

Due to the inclusion of scattering processes, two issues arise in the thermal averages. Since we do not consider loop corrections to the two-body decay or $1\to 3$ processes $\widetilde{b}\to\chi bg$, we cannot use them to cancel soft contributions from $g\widetilde{b}\leftrightarrow b\chi$ scatterings (of course, the $\gamma$ scattering also has this problem, as noted by ∗ in Table 2). Instead, we regularize these processes by imposing a cut on the minimal process energy of $s_{\text{min}}=(m_{\widetilde{b}}+m_{\text{cut}})^{2}$, with fiducial value $m_{\text{cut}}=0.5\,\text{GeV}$. We checked that our numerical results are stable when varying $m_{\text{cut}}$ over a wide range, see Fig. 6, indicating that the bulk of the scattering processes occurs at energies above the $b$ mass. On top of this, we find that in processes of the type $\widetilde{b}\bar{b}\leftrightarrow\chi H$ the $b$-quark in the $t$-channel is allowed to go on-shell for some center-of-mass energies (the affected processes are marked with ∗∗ in Tab. 2). This corresponds to a double counting of the on-shell two-body decay. We choose to suppress the on-shell part by introducing a large Breit-Wigner width for the $b$-quark, taking $\Gamma_{b}=m_{b}$ in our numerical calculations. Since this issue occurs only in processes involving weak scale particles ($H,W,Z$) and the scattering rate is dominated by gluons the precise value for the width does not have an appreciable impact on the results.

Appendix B Sommerfeld enhancement

In the presence of light degrees of freedom non-perturbative corrections to the annihilation rates are known to become relevant in the non-relativistic limit Hisano et al. (2004, 2005). Between pairs of color charged particles the exchange of gluons generates a potential which modifies the wave function of the initial state particles and leads to a non-negligible correction of the tree-level cross section De Simone et al. (2014); Ibarra et al. (2015); Hryczuk (2011); El Hedri et al. (2017).

To leading order the effect of the QCD potential can be described by a Coulomb-like potential Fischler (1977)

[TABLE]

where $C_{R}$ and $C_{R^{\prime}}$ denote the Casimir coefficients of the incoming particles while $C_{Q}$ is the Casimir coefficient of the final state. For a general Coulomb-potential with $V(r)=\alpha/r$ the $s$-wave Sommerfeld correction factor $S_{0}$ is given by De Simone et al. (2014)

[TABLE]

where $\beta=v/2$ and the total annihilation cross section of particles moving in this potential is given by $\sigma_{Somm}=S_{0}\cdot\sigma_{tree}$.555In principle the Sommerfeld factors have to be determined separately for each partial wave. For the model considered here the total Sommerfeld effect can be approximated to good accuracy by applying the $s$-wave correction to the full cross section. For final states which are exclusively in a singlet, i.e. $ZZ,W^{+}W^{-},\gamma\gamma$, or an octet representation, i.e. $\gamma g,Zg$, the enhancement is given by Eq. (13) with $\alpha=-4/3\alpha_{s}$ or $\alpha=1/6\alpha_{s}$, respectively. The $gg$ final state is slightly more complicated since it can be in a singlet or octet representation. After summing over the different contributions the total Sommerfeld correction factor for this case reads De Simone et al. (2014)

[TABLE]

Since this channel dominates the annihilation rates by orders of magnitude, we only take the correction for annihilation to gluons into account.

Appendix C Kinetic equilibrium

In this section we compare solutions of the differential, momentum-dependent Boltzmann equation for $\chi$ with the integrated Boltzmann equation used in the main text. The latter relies on the assumption of kinetic equilibrium. This assumption may be questionable for the range of parameters we are interested in, because the leading interactions of $\chi$ are the conversion processes that become inefficient around the time of freeze-out.

If a particle $X$ is in kinetic equilibrium with the thermal bath of SM particles at temperature $T$, its distribution function is given by

[TABLE]

where $f_{X}(p,t)$ is the phase-space density. We assume that due to efficient coupling to the SM, the mediator $\widetilde{b}$ is in kinetic equilibrium for all relevant times. When dropping the assumption of kinetic equilibrium for $\chi$, the unintegrated Boltzmann equation for $f_{\chi}$,

[TABLE]

has to be solved, see e.g. D’Agnolo et al. (2017); Binder et al. (2017). The collision operator $C$ determines the differential interaction rate of each momentum mode of $\chi$. Introducing the notation

[TABLE]

the Liouville operator on the left-hand side of Eq. (16) can be brought into the form

[TABLE]

In the last line we assume $h_{\text{eff}}$ to be constant during freeze-out which is approximately satisfied above temperatures of $\mathcal{O}(\text{GeV})$. Here we adapt the notation of Edsjö and Gondolo (1997) for the effective degrees of freedom $g_{\text{eff}}$ and $h_{\text{eff}}$ associated with the energy and entropy densities, respectively.

C.1 Collisional operator for conversion

For the purpose of comparing the integrated with the differential Boltzmann equation, we focus on the dominant interaction processes that drive the distribution function towards kinetic equilibrium. For the parameter range we are interested in, these are conversion processes, in particular (inverse) decays $\chi b\leftrightarrow\widetilde{b}$, and, for earlier times, also inelastic scatterings $\chi A\leftrightarrow\widetilde{b}B$. We neglect the small contributions from $\chi\chi$ annihilation, $\chi\widetilde{b}$ coannihilation and elastic scattering $\chi b\leftrightarrow\chi b$. Note that if we would take these additional processes into account, one expects them to drive the distribution function closer to its equilibrium shape. Therefore the following analysis may be regarded as a conservative estimate of the deviations from kinetic equilibrium.

With these assumptions Eq. (16) becomes a linear differential equation and the remaining contributions to $C\left[f_{\chi}\right]$ are the (inverse) decay $\chi b\leftrightarrow\widetilde{b}$,

[TABLE]

and the inelastic scatterings $\chi A\leftrightarrow\widetilde{b}B$,

[TABLE]

where $\mathrm{d}\Pi_{X}=g_{X}\mathrm{d}^{3}p_{X}/\left((2\pi)^{3}E_{X}\right)$. Apart from the different integrations, the terms in brackets differ from those in the case of the integrated Boltzmann equation. For the decay term we find

[TABLE]

using the relation of detailed balance for the equilibrium distributions. A similar simplification can be performed for inelastic scattering such that $C$ factorizes into

[TABLE]

As before, we neglect quantum statistical factors in the calculation. The unintegrated Boltzmann equation (16) can hence be written in the form

[TABLE]

This ordinary differential equation together with the boundary condition

[TABLE]

can be solved with separation of variables and variation of constants. The result reads (cf. D’Agnolo et al. (2017))

[TABLE]

C.2 Simplifications for the numerical solution

In order to trace the evolution of the number density

[TABLE]

the phase space density $f_{\chi}(q,x)$ has to be computed for a large number of momentum modes $q$ and temperature parameters $x$ providing a sufficiently accurate numerical approximation. This is a computationally expensive task. In the following we therefore introduce analytic simplifications of the collision operators.

For the two-body decay we can find the analytic result for the collision operator, neglecting the $b$-quark mass,

[TABLE]

where

[TABLE]

and $g_{\chi}=2$, $g_{b}=6$, $g_{\tilde{b}}=3$. The factor of $2$ in the numerator of Eq. (27) accounts for the two different processes when considering the mediator and its anti-particle, which is not included in $g_{\tilde{b}}$.

For the inelastic scattering we consider the process $\chi b\leftrightarrow\widetilde{b}g$. Neglecting, again, the bottom mass we can express the collision operator as

[TABLE]

where

[TABLE]

and

[TABLE]

Note that while $\widetilde{C}_{12}$ is an analytic function, $\widetilde{C}_{22}$ has to be numerically evaluated. We choose to precompute $\widetilde{C}_{22}$ on a two-dimensional grid in $x$ and $q=p_{\chi}/T$ and use an interpolation for the numerical evaluation of Eq. (25). We provide most of the discussion for the case when taking into account the decay term only, which is expected to capture the main effects. We comment on the effect of scatterings at the end of this section.

C.3 Iterative solution of the coupled system

The solution Eq. (25) of the Boltzmann equation for $f_{\chi}(q,x)$ requires as an input the evolution of the mediator abundance, $Y_{\widetilde{b}}(x)$. The latter can be obtained by solving the corresponding integrated Boltzmann equation, which in turn involves $Y_{\chi}(x)$, that is determined by integrating $f_{\chi}(q,x)$ over all momentum modes. Therefore the equations for $f_{\chi}(q,x)$ and $Y_{\widetilde{b}}(x)$ form a coupled set of equations.

Here we solve this coupled set of differential equations in an iterative process. We start with an initial “guess” for $Y_{\widetilde{b}}(x)$, which we take to be the solution when assuming kinetic equilibrium (see below for a discussion of different choices). We then solve for $f_{\chi}(q,x)$ on a momentum-grid, and numerically compute $Y_{\chi}(x)$ using Eq. (25) as described in the last subsection. With this solution for $Y_{\chi}(x)$ we recalculate $Y_{\widetilde{b}}(x)$ using the integrated Boltzmann equation. We subsequently iterate between solving for $f_{\chi}(q,x)$ and $Y_{\widetilde{b}}(x)$, until we encounter sufficient convergence. In order to solve the differential Boltzmann equation in an acceptable time, we neglect the bottom mass and choose $h_{\text{eff}}$ and $g_{\text{eff}}$ to be evaluated at $x=50$ and constant for all times. We do not expect a strong dependence on these simplifications.

The resulting evolution of the abundance $Y_{\chi}(x)$ for the benchmark point $m_{\chi}=500\,$GeV, $m_{\widetilde{b}}=510\,$GeV is shown in Fig. 7 (upper panel) as a red solid curve. We compare the result to the solution of the coupled integrated Boltzmann equation (red dotted curve) obtained under the same approximations. We adjust the coupling $\lambda_{\chi}=4.03\times 10^{-7}$ so as to obtain the measured DM relic density for the solution of the coupled integrated Boltzmann equation. The lower panel of Fig. 7 shows the ratio of the differential and integrated solutions for $Y_{\chi}(x)$. While the dark matter abundance differs by up to a factor of two at intermediate times, the final relic abundance agrees well with the corresponding result when assuming kinetic equilibrium, with deviations below the 10% level.

The main reason is that, for the process and the kinematical situation that is relevant here, the collision term does not depend strongly on the momentum mode, see Fig. 8 (dot-dashed lines). In the same figure, we also show the result for $f_{\chi}(q,x)$ at various times $x$ (blue lines), which indeed differs from an equilibrium distribution (indicated by the red dashed lines) at intermediate times (upper and middle panel in Fig. 8). Nevertheless, around the time when the dark matter abundance freezes out, the remaining decays of thermalized $\widetilde{b}$ tend to restore an equilibrium distribution (lowest panel).

It is interesting to observe that the total abundance obtained from the unintegrated Boltzmann equation is slightly below the result when assuming kinetic equilibrium. This can also be understood from Fig. 8. For high temperatures, the momentum modes obtained from the differential solution essentially change only due to redshift. In contrast, the kinetic equilibrium distribution populates somewhat higher modes. By the time the conversion gets efficient, the collision term is larger for smaller momentum modes. Therefore, the conversion into $\widetilde{b}$s is stronger for the differential solution, rendering a slightly smaller abundance.

Let us now discuss the impact of the initial “guess” for $Y_{\widetilde{b}}(x)$ used for the iterative solution. We check that the converged result is independent of the starting point of the iteration by using two rather different initial abundances. First, we use the equilibrium abundance $Y_{\widetilde{b}}^{\mathrm{eq}}(x)$ as a starting point. The results are shown in the left panel of Fig. 9. The evolution of $Y_{\chi}(x)$ obtained in the first iteration step is shown by the thick orange line, and the successive iterations are indicated by the thin orange lines. The final, converged result (thick red line) agrees well with the result obtained in Fig. 7. The same is true for $Y_{\widetilde{b}}(x)$ (solid blue line). On the other hand, we would like to point out that the first iteration and the converged result are rather far apart. This means that it is crucial to solve for the coupled set of equations, allowing for deviations $Y_{\widetilde{b}}(x)\not=Y_{\widetilde{b}}^{\mathrm{eq}}(x)$. For curiosity, we note that if one would compare the differential with the integrated result for $Y_{\chi}(x)$ while fixing $Y_{\widetilde{b}}(x)=Y_{\widetilde{b}}^{\mathrm{eq}}(x)$, one would find an ${\cal O}(10)$ difference in the final abundance (see orange dotted versus solid line in the left panel of Fig. 9), while the corresponding difference for the converged results is below $\sim 10\%$. Hence, the partial freeze-out of the mediator $\widetilde{b}$ and its subsequent decay into $\chi$ are crucial for the conclusion that the impact of deviations from kinetic equilibrium on the relic density is small.

Second, we consider an extreme possibility and initially set $Y_{\chi}(x)$ to be constant and equal to the relativistic equilibrium density. In this case we start the iteration with the computation of $Y_{\widetilde{b}}(x)$. The resulting iterative solutions for $Y_{\chi}(x)$ are shown in right panel of Fig. 9 (orange lines). Again, the converged result for $Y_{\chi}(x)$ (red solid line) and $Y_{\widetilde{b}}(x)$ (solid blue line) agree well with those shown in Fig. 7.

The convergence of the final relic density for the three different choices of initial abundances is shown in Fig. 10. Indeed, the converged results agree, indicating that the iterative scheme is stable and leads to a unique result.

Next we want to check if the situation changes drastically when including also $2\rightarrow 2$ scattering processes. Due to the increase in numerical complexity described above, we consider the leading process $\chi b\leftrightarrow\widetilde{b}g$ expected to capture the main effects. In order to estimate the physical contributions from hard scatterings, we perform regularizations on the level of the scattering cross section by introducing a cut-off $s_{\min}=(m_{\widetilde{b}}+1\,\text{GeV})^{2}$ and additionally a regulator at the matrix element level of $1/t^{2}\rightarrow 1/(t^{2}+(1\,\text{GeV})^{4})$. In addition, we restrict the integration over the angle $\theta_{t}$ between $b$ and $g$ in the center-of-mass frame to $\cos\theta_{t}\in[-0.9,0.9]$.

Again, we solve the coupled system in an iterative approach as described above, but taking scatterings into account. As before, we then compare the converged result for the final relic density with the corresponding result obtained when assuming kinetic equilibrium. We find that the relative deviations in the resulting relic density stay below $10\%$. Furthermore, the deviation for $Y_{\chi}(x)$ for intermediate times becomes smaller. This is expected, because scatterings increase the conversion rates at smaller $x$.

Altogether, we find that the impact of deviations from kinetic equilibrium on the final relic abundance is rather mild, below the $10\%$ level. This justifies using integrated Boltzmann equations for $Y_{\chi}(x)$ and $Y_{\widetilde{b}}(x)$.