How Heavy can Neutralino Dark Matter be?

TL;DR

This paper investigates the maximum possible mass of neutralino dark matter consistent with observed relic abundance, focusing on scenarios with near-degenerate masses leading to a second freeze-out after the QCD phase transition.

Contribution

It introduces the concept of a second freeze-out due to near-degenerate neutralino and colored sparticles, allowing for neutralino masses well above 100 TeV.

Findings

Neutralino mass can exceed 100 TeV with near-degenerate sparticles.

Second freeze-out significantly enhances dark matter annihilation.

Implications for dark matter detection strategies.

Abstract

What is the upper limit of the mass of the neutralino dark matter whose thermal relic is consistent with the observation? If the neutralino dark matter and colored sparticles are extremely degenerated in mass, with a mass difference less than the QCD scale, the dark matter annihilation is significantly increased and enjoys the "second freeze-out" after the QCD phase transition. In this case, the neutralino dark matter with a mass much greater than 100 TeV can realize the correct dark matter abundance. We study the dark matter abundance and its detection in the case of such highly degenerated mass spectrum of the neutralino dark matter and colored supersymmetric particles.

Figures

Tables

Table 1. Table 1 : The coefficients for the direct detection in Eqs. ( 34 ) and ( 35 ).

nucleon	$\tilde{u}$	$\tilde{d}$	$\tilde{s}$	$\tilde{c}$	$\tilde{b}$
proton	$7.2\times {10}^{-46}$	$1.6\times {10}^{-47}$	$3.0\times {10}^{-48}$	$1.8\times {10}^{-47}$	$9.5\times {10}^{-52}$
neutron	$2.4\times {10}^{-46}$	$4.7\times {10}^{-47}$	$3.0\times {10}^{-48}$	$1.8\times {10}^{-47}$	$9.5\times {10}^{-52}$

Equations

$$\frac{d\tilde{Y}}{dx}=-\frac{xs}{H(m_{{\tilde{\chi }}_1^0})}\left(1+\frac{T}{3g_{\ast s}}\frac{d{g}_{\ast s}}{dT}\right)⟨{\sigma }_{eff}v⟩\left({\tilde{Y}}^2-{\tilde{Y}}_{eq}^2\right),$$

$$x\equiv \frac{m_{{\tilde{\chi }}_1^0}}{T},s=\frac{2{\pi }^2}{45}{g}_{\ast s}{T}^3,H(m_{{\tilde{\chi }}_1^0})={\left(\frac{4{\pi }^3{G}_N{g}_{\ast }}{45}\right)}^{\frac{1}{2}}{m}_{{\tilde{\chi }}_1^0}^2.$$

$$⟨{\sigma }_{eff}v⟩=⟨\sigma v{⟩}_{\tilde{q}{\tilde{q}}^{\ast }}\frac{g_{\tilde{q}{\tilde{q}}^{\ast }}^2(1+\Delta {)}^3{e}^{-2\Delta x}}{g_{eff}^2},$$

$$\Delta \equiv (m_{\tilde{q}}-m_{{\tilde{\chi }}_1^0})/m_{{\tilde{\chi }}_1^0},g_{eff}\equiv g_{{\tilde{\chi }}_1^0}+g_{\tilde{q}{\tilde{q}}^{\ast }}(1+\Delta {)}^{3/2}{e}^{-\Delta x},$$

$$⟨\sigma v{⟩}_{\tilde{q}{\tilde{q}}^{\ast }}=⟨\sigma v(\tilde{q}{\tilde{q}}^{\ast }\to gg,q\bar{q})⟩+⟨\sigma v{⟩}_{bsf}\frac{⟨\Gamma {⟩}_{\eta }}{⟨\Gamma {⟩}_{\eta }+⟨\Gamma {⟩}_{dis}}.$$

$${\sigma }_{\text{form}}((\tilde{q}\bar{q})+({\tilde{q}}^{\ast }q)\to ({\tilde{q}}^{\ast }\tilde{q})+X)\sim \pi {\Lambda }_{\text{had}}^{-2}.$$

$${\sigma }_{\text{form}}\sim \sum _L^{L_{\text{max}}}\frac{4\pi (2L+1)}{p^2}\simeq \pi \frac{L_{\text{max}}^2}{p^2}=\pi R_{\text{had}}^2.$$

$$L\sim 10{\left(\frac{m_{\tilde{q}}}{\text{TeV}}\right)}^{\frac{1}{2}}{\left(\frac{T}{T_c}\right)}^{\frac{1}{2}}{\left(\frac{{\Lambda }_{\mathrm{had}}}{1\mathrm{GeV}}\right)}^{-1}.$$

$${\Gamma }_{\text{de-excitation}}\sim \frac{Q^2\alpha {\Lambda }_{\text{had}}^3}{{\alpha }_{\text{QCD}}({\Lambda }_{\text{had}})m_{\tilde{q}}^2},$$

$$({\tilde{q}}^{\ast }\tilde{q})+X\to ({\tilde{q}}^{\ast }\tilde{q}{)}^{\prime }+X.$$

$${\sigma }_{\text{de-excitation, col, π}}$$

$${\sigma }_{\text{de-excitation, col, γ}}$$

$${\sigma }_{\text{form}}^{\prime }(T)=\mathrm{BF}(({\tilde{q}}^{\ast }\tilde{q})\to \mathrm{S}{\mathrm{M}}^{\prime }\mathrm{s})\times {\sigma }_{\text{form}},$$

$$\mathrm{BF}(({\tilde{q}}^{\ast }\tilde{q})\to \mathrm{S}{\mathrm{M}}^{\prime }\mathrm{s})=\frac{⟨{\Gamma }_{\text{de-excitation}}⟩}{⟨{\Gamma }_{\text{de-excitation}}⟩+⟨{\Gamma }_{\text{destruction}}⟩},$$

$${\Gamma }_{\text{destruction}}={\Gamma }_{\text{decay}}+{\Gamma }_{\text{dest, coll}}.$$

$${\Gamma }_{\mathrm{decay}}=2\times \Gamma (\tilde{q}\to {\tilde{\chi }}_1^0+\mathrm{SM}).$$

$${\sigma }_{\text{dest, coll, break}}=\frac{p^2}{{\omega }^2}\frac{g_{\tilde{q}\bar{q}}^2}{g_{\tilde{q}{\tilde{q}}^{\star }}{g}_{\pi }}{\sigma }_{\text{form}},$$

$${\sigma }_{\text{dest, coll, break}}\sim \pi R_{\mathrm{had}}^2,$$

$${\Gamma }_{\text{dest, coll}}\simeq {\sigma }_{\text{dest, coll, break}}\times n_{\pi ,E\gtrsim {\Lambda }_{\mathrm{had}}}\sim \pi R_{\mathrm{had}}^2{\Lambda }_{\mathrm{had}}^{2}T\exp (-{\Lambda }_{\mathrm{had}}/T),$$

$${\Gamma }_{\text{annihilate}}\sim \frac{{\alpha }_{\text{QCD}}^2(m_{\tilde{q}})}{m_{\tilde{q}}^2}{\left|\psi (0)\right|}^2,$$

$$m_{{\tilde{u}}_L}-m_{{\tilde{d}}_L}\simeq \frac{(m_u^2-m_d^2)+\cos (2\beta )m_W^2}{2m_{\tilde{q}}}\sim 100\mathrm{MeV}{\left(\frac{m_{\tilde{q}}}{10\mathrm{TeV}}\right)}^{-1},$$

$$\Delta m\sim 70\text{MeV}.$$

$${\Gamma }_{\text{conv}}\sim c^2\alpha T\times \{\begin{array}{ll}\exp (-m_h/T)& \mathrm{for}{m}_h\gg \Delta m\\ \exp (-\Delta m/T)& \mathrm{for}{m}_h\ll \Delta m\end{array}$$

$${\Gamma }_{\text{conv, tree}}\sim c^2\alpha \frac{T^5}{m_{\tilde{q}}^4}\exp (-m_h/T)$$

$$m_{\tilde{q}}\lesssim O(1)\text{TeV}$$

$${\Gamma }_{\text{conv, loop}}\sim \frac{T^3{m}_t^2}{(16{\pi }^2{)}^2{m}_{\tilde{t}}^4}\exp (-m_{\pi }/T).$$

$$m_{\tilde{t}}<O(10)\text{TeV}.$$

$$\frac{d{n}_{\tilde{q}}}{dt}=$$

$$\frac{d{n}_{{\tilde{\chi }}_1^0}}{dt}=$$

$${\mathcal{L}}_{\text{eff}}=f_q{m}_q\overline{\stackrel{~}{B}}\tilde{B}\bar{q}q+\frac{g_q^{(1)}}{m_{{\tilde{\chi }}_1^0}}\overline{\stackrel{~}{B}}i{\partial }^{\mu }{\gamma }^{\nu }\tilde{B}{\mathcal{O}}_{\mu \nu }^q+\frac{g_q^{(2)}}{m_{{\tilde{\chi }}_1^0}^2}\overline{\stackrel{~}{B}}i{\partial }^{\mu }i{\partial }^{\nu }\tilde{B}{\mathcal{O}}_{\mu \nu }^q,$$

$$f_q$$

Full text

IPMU18-0201

**How Heavy can Neutralino Dark Matter be? **

Hajime Fukuda1, Feng Luo2 and Satoshi Shirai1

1*Kavli Institute for the Physics and Mathematics of the Universe (WPI),

The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa 277-8583, Japan

2School of Physics and Astronomy, Sun Yat-sen University, Zhuhai 519082, China *

1 Introduction

The supersymmetric (SUSY) standard model (SSM) is one of the most attractive candidate of the physics beyond the standard model. In the SSM, the lightest SUSY particle (LSP) is stable if the $R$-parity is conserved. The most important feature of the SSM is that the LSP neutralino can be a weakly interacting massive particle (WIMP) dark matter (DM).

The WIMP dark matter abundance is mainly determined by the freeze-out mechanism, if the reheating temperature of the Universe is high enough and there is no additional entropy production such as a moduli field decay. An interesting feature of the freeze-out mechanism is that the DM abundance can be determined by the DM annihilation rate at the early Universe and predicted solely from the low-energy property of the DM sector, regardless of the initial condition of the Universe. Generally speaking, as the DM mass gets larger, the DM abundance also gets larger. In order to keep the thermal relic DM density being consistent with the observation, we can obtain an upper limit of the DM mass. This upper limit is a very important guideline for the DM search experiments.

In the case of the SSM, the upper limit of the LSP DM mass depends on the details of the mass spectrum. If the lightest neutralino is the wino, the upper limit of the LSP mass is around 3 TeV [1]. In fact, this wino LSP is a generic prediction of the anomaly mediation model [2, 3]. This model is the simplest realization of the (mini-)split SUSY [4, *Wells:2004di, 6, *Giudice:2004tc, *ArkaniHamed:2004yi, *ArkaniHamed:2005yv], which is getting more and more attention [10, *Hall:2012zp, *Nomura:2014asa, 13, *Ibe:2012hu, 15, 16]. The wino has rich signatures of (in)direct dark matter [17, *Hisano:2004ds, *Hisano:2005ec, 20, *Hisano:2010ct, *Hisano:2012wm, *Hisano:2015rsa] and collider searches [24, *Buckley:2009kv, 26, *Asai:2008sk, *Asai:2008im]. The detection of the 3 TeV wino will be within reach within the next decade.

In the case of a Higgsino dark matter, the upper limit is around 1 TeV. The Higgsino DM is also intensively studied. Future hadron collider and direct detection experiments will enable almost complementary searches up to the 1 TeV Higgsino [29, 30, *Fukuda:2017jmk, 32, *Matsumoto:2018ioi, 34, 35].

If the SSM mass spectrum is fine-tuned, the relic abundance of the LSP can be exceptionally smaller, accordingly, the upper limit of the DM mass is also increased [36]. An example is the so-called Higgs-funnel region, where the neutralino/chargino mass is close to half of the heavier Higgs mass. The annihilation process of the neutralino is significantly enhanced as the $S$-channel process hits the pole of the heavy Higgs propagator. In this case the upper limit of the DM mass can be around 10 TeV [37, 38, 39, 40].

Another exceptional case is the coannihilation. If the lightest neutralino and next-to-LSP (NLSP) is almost degenerated in mass, the effective annihilation rate of the neutralino LSP can be enhanced through the NLSP annihilation. A prime example is coannihilation with colored sparticles. In particular the coannihilation with gluino drastically increases the effective annihilation rate and the upper limit of the DM is around 10 TeV [41, 42, 43, 44, 45, *Nagata:2017gci]. This 10 TeV gluino can be also probed at a future hadron collider [47].

In this paper, we further investigate the coannihilation with colored sparticles. If the mass difference between the colored NLSP and DM is less than around the QCD scale, the annihilation rate of the DM is affected by non-perturbative QCD effects. In this case, in addition to the conventional freeze-out which occurs at temperature around $m_{\rm DM}/30$, there is another annihilation stage in the QCD phase transition era. Due to the non-perturbative strong interaction, the effective annihilation cross section drastically increases, allowing even PeV scale dark matter. This large mass is above the upper-bound on the WIMP mass $\sim 100$ TeV, which is a generic upper-limit of the WIMP DM mass based on $S$-wave unitarity bound [48]. This DM reduction mechanism is studied in the context of the non-standard-model strong interaction [49] and suggested in the context of the SSM [38]. We revisit this possibility, paying attention to the non-perturbative annihilation and chemical equilibrium between the LSP and colored sparticles. We also study the detection of such a heavy DM. Due to the small mass difference with the colored NLSP, the direct detection rate is significantly increased.

2 Coannihilation

In this section, we discuss the coannihilation in the case that the colored NLSP is highly degenerated with the neutralino LSP dark matter in mass. First, we discuss the annihilation processes of the colored NLSP in the high temperature and QCD phase transition eras. Then, we study the chemical equilibrium condition of the neutralino LSP and the colored NLSP. In this paper, we mainly focus on the squark-neutralino coannihilation.

2.1 Squark Annihilation

If the mass difference between the neutralino LSP and squark NLSP is less than $\sim 1\,\text{GeV}$, the cosmological evolution of the dark matter shows interesting features. Once the temperature gets lower than around $m_{\rm DM}/30$, the dark matter coannihilation process freezes out, as in the case of conventional WIMP freeze-out scenario. However, if the colored sparticle is still abundant in the QCD phase transition era, non-perturbative QCD effect restarts the coannihilation. When the temperature gets low enough, the coannihilation process freezes out again. In the following, we describe these two perturbative and non-perturbative freeze-outs 111To be precise, the word “perturbative” is not very appropriate here, since we include the Sommerfeld and the bound-state effects, which are inherently non-perturbative and important before the QCD phase transition, in the “perturbative annihilation”. However, in order to be concise when referring to the two phases before and after the QCD phase transition, as well as to highlight the key differences of them, we choose to still use the word “perturbative annihilation” to refer to the annihilation process through the QCD in the Coulomb phase and “non-perturbative annihilation” to the QCD in the confined phase..

2.1.1 Perturbative Annihilation

We first review the squark-neutralino coannihilation process in the era long before the QCD phase transition. First of all, coannihilation can happen only if the interconversion rate between a squark and a neutralino (including decays/inverse decays, and conversions by scattering with the standard model particles in the thermal bath) is sufficiently large compared to the Hubble expansion rate, otherwise the two particle species will freeze out independently. As we will discuss in section 2.2, this condition of coannihilation can be satisfied in this era, so that to a very good approximation the number density ratio of squarks and neutralinos equals to their thermal equilibrium number density ratio at the same temperature, that is, $n_{\tilde{q}}/n_{\tilde{\chi}^{0}_{1}}=n_{\tilde{q}}^{eq}/n_{\tilde{\chi}^{0}_{1}}^{eq}$. Therefore we can use a single Boltzmann equation to track the evolution of the total yield, which is the number density sum of squarks, antisquarks and neutralinos divided by the entropy density 222We assume that there is no asymmetry between $\tilde{q}$ and $\tilde{q}^{*}$, and we take $n_{\tilde{q}}=n_{\tilde{q}^{*}}$ in our calculation., $\tilde{Y}\equiv(n_{\tilde{q}}+n_{\tilde{q}^{*}}+n_{\tilde{\chi}^{0}_{1}})/s$,

[TABLE]

where

[TABLE]

$G_{N}$ is the gravitational constant, $m_{\tilde{\chi}^{0}_{1}}$ is the neutralino LSP mass, $g_{\ast s}$ and $g_{\ast}$ are the numbers of effectively massless degrees of freedom associated with the entropy density and the energy density, respectively.

In principle, the thermally-averaged effective annihilation cross section, $\langle\sigma_{eff}v\rangle$, should include contributions of all (co)annihilation channels which can change $\tilde{Y}$. However, for the case of highly degenerated colored NLSP and neutralino LSP we focus in this paper, the dominant contributions to $\langle\sigma_{eff}v\rangle$ come from $\tilde{q}\tilde{q}^{*}$ annihilations. Also, since we are not committed to the details of the SSM mass spectrum which affects the size of electroweak (co)annihilation channels through for example squark left-right mixing angle, we include only the pure strong interaction channels, $\tilde{q}\tilde{q}^{*}\to gg$ and $\tilde{q}\tilde{q}^{*}\to q\bar{q}$, in calculating $\langle\sigma_{eff}v\rangle$. These two processes are universal for all squark flavors and are usually the dominant ones. 333Ref. [50] discussed the case when the dominant annihilation channels are $\tilde{t}\tilde{t}^{*}\to hh,W^{+}W^{-},ZZ$, when the interactions between stop and Goldstone bosons are large. We also include in $\langle\sigma_{eff}v\rangle$ the $\tilde{q}\tilde{q}^{*}$ Sommerfeld and bound-state effects as discussed in detail in Ref. [44]444The Sommerfeld and/or bound-state effects are also recently discussed in [51, *Ellis:2014ipa, *vonHarling:2014kha, *Ibarra:2015nca, *Cirelli:2016rnw, *Kim:2016zyy, *Kim:2016kxt, *ElHedri:2016onc, *Petraki:2016cnz, *ElHedri:2017nny, *Pierce:2017suq, *Keung:2017kot, *Garny:2017rxs, *Harz:2017dlj, *Biondini:2018pwp, *ElHedri:2018atj, *Binder:2018znk, *Biondini:2018xor, *Biondini:2018ovz, *Harz:2019rro, *Mitridate:2017izz, 72].. With these considerations, we use

[TABLE]

where

[TABLE]

where $g_{\tilde{\chi}^{0}_{1}}=2$ is the degrees of freedom for a neutralino, assumed to be a bino. $g_{\tilde{q}\tilde{q}^{*}}=6$ is the sum of the degrees of freedom for squark and antisquark. $\langle\sigma v\rangle_{\tilde{q}\tilde{q}^{*}}$ is the thermally-averaged annihilation cross section times relative velocity 555As is noted in [44], since $g_{\tilde{q}\tilde{q}^{*}}=6$ is the sum of the degrees of freedom for squark and antisquark, a factor of $1/2$ needs to be put in $\langle\sigma v\rangle_{\tilde{q}\tilde{q}^{*}}$, if using the usual spin and color averaged cross section. This is explained in the Appendix of [73]., given as

[TABLE]

In the above expression, the first term is the sum of the Sommerfeld enhanced $\tilde{q}\tilde{q}^{*}\to gg$ and $\tilde{q}\tilde{q}^{*}\to q\bar{q}$ cross sections, and the second term is the contribution from the $\tilde{q}\tilde{q}^{*}$ bound-state effect 666We update in this work the bound-state formation cross section and dissociation rate calculations compared to [44] by including the contribution from the diagram of gluon emission from gluons. This type of diagram is only present when the gauge interaction is non-abelian, and was first pointed out in [74]. We follow [72] for the calculation., where $\langle\sigma v\rangle_{bsf}$ is the bound state formation cross section, $\langle\Gamma\rangle_{\eta}$ is the bound state annihilation decay rate, and $\langle\Gamma\rangle_{dis}$ is the bound state dissociation rate. All quantities are thermally averaged. We note that the bound-state effect we are discussing in this subsection is in the perturbative regime of QCD at high temperature, in contrast to the QCD bound state which we will discuss in the next subsection, when the QCD enters the non-perturbative regime after the QCD phase transition. For the detail calculations in the perturbative regime, we refer the readers to Ref. [44, 72].

We evaluate Eq. (1) from $T\sim m_{\tilde{\chi}^{0}_{1}}$ to lower temperatures. The usual freeze-out (when $\tilde{Y}/\tilde{Y}_{eq}-1$ becomes larger than 1) happens around $T\sim m_{\tilde{\chi}^{0}_{1}}/30$. After that, $\tilde{Y}$ continues decreasing, and in particular the bound-state effect becomes effective when $T$ is around and below the bound-state binding energy, $E_{B}\sim m_{\tilde{\chi}^{0}_{1}}/100$. $\tilde{Y}$ essentially ceases to decrease after $T\sim m_{\tilde{\chi}^{0}_{1}}/1000$, until the new annihilation process kicks in after the QCD phase transition which we describe in the next subsection.

2.1.2 Non-perturbative Annihilation

After the first freeze-out, as the Universe gets cooler, the QCD phase transition occurs. With non-perturbative QCD effect, the squark annihilation rate drastically increases. In this subsection, we review the annihilation process using the non-perturbative effect of the QCD based on Ref. [75] 777Also, see Ref. [76] for the analytic calculation, which is basically consistent with Ref. [75].. The annihilation process is in the following way. First, a squark forms a QCD bound state with a quark/gluon, which are abundant in the environment, and becomes a SUSY hadron such as a mesino. The SUSY hadron has a large geometrical size $\Lambda_{\rm QCD}^{-1}$ and large scattering cross section with another SUSY hadron. Once such a scattering occurs, a squark-squark bound state ${\tilde{q}}^{*}\tilde{q}$ of a large angular momentum forms. Then, the bound state de-excites into the ground state. Finally, the $S$-wave state annihilates into quarks and gluons. If the rates of de-excitation and annihilation of the ${\tilde{q}}^{*}\tilde{q}$ bound state are fast enough, the annihilation cross section of the squarks is virtually equivalent to the bound state formation cross section, which is much larger than the perturbative annihilation. In Fig. 1, we show a schematic picture of this annihilation process. In the following, we examine the rates of these non-perturbative processes.

First, we discuss the ${\tilde{q}}^{*}\tilde{q}$ bound state formation. Soon after the QCD phase transition at temperature $T_{c}\simeq 200\,\text{MeV}$, the squarks start to form mesinos, for example, $\tilde{q}\bar{q}$ and ${\tilde{q}}^{*}{q}$. Since the squark is much heavier than the QCD scale, the squark is localized at the center in the mesino bound state, $\sim m_{\tilde{q}}^{-1}$ in radius, and the quark/gluon partons surround it with a radius of $R_{\text{had}}=\Lambda_{\text{had}}^{-1}$, where $\Lambda_{\text{had}}$ is the hadronic mass scale. Then, these heavy hadrons interact with each other, forming a squarkonium ${\tilde{q}}^{*}\tilde{q}$, a squark-antisquark bound state, where the cross section is the geometrical one, 888Note that Ref. [77] has claimed that $\tilde{q}qq$ bound states are involved in the squarkonium formation process instead of an anti-mesino, $\tilde{q}\bar{q}$. However, they concluded that even including $\tilde{q}qq$ state, the squarkonium formation rate is not so different. Thus, we simply use Eq. (6) here.

[TABLE]

This geometrical cross section can be explained in another way. As long as the impact parameter of the collision is less than the size of the mesino, $R_{\text{had}}$, we may assume the interaction is strong enough. Thus, the cross section for a partial wave $L$ saturates the unitarity bound if $L<L_{\text{max}}\equiv pR_{\text{had}}$, where $p$ is the squark momentum. Summing up the partial wave cross section up to $L_{\text{max}}$, we obtain the geometrical cross section [78];

[TABLE]

Ref. [76] also discusses this behavior analytically.

Then let us estimate $R_{\text{had}}=\Lambda_{\text{had}}^{-1}$. In Ref. [75], $\Lambda_{\text{had}}$ is taken to be $1\,\text{GeV}$. On the other hand, the pion charge radius is experimentally $0.67\,\text{fm}\sim(300\,\text{MeV})^{-1}$ [79], which suggests to use $\Lambda_{\text{had}}\sim 300\,\text{MeV}$. Also, Ref. [77] uses a slightly larger value, $R_{\text{had}}\sim 0.5\,\text{fm}\sim(400\,\text{MeV})^{-1}$. Moreover, Ref. [80] also claims that such larger value can be available 999They also discuss the energy dependence of the cross section. However, in our region of interest, the annihilation process is almost instantaneous in the temperature and we do not discuss it.. Thus, combining these contexts, we take the $\Lambda_{\rm had}$ as a free parameter which varies from $400$ MeV to $1$ GeV in the following analysis.

Next, we discuss the de-excitation process. Initially, the bound state has a very large angular momentum: since the momentum of the squark is $p\sim\sqrt{2m_{\tilde{q}}T}$, and the impact parameter is $r\sim\Lambda_{\text{had}}^{-1}$, the angular momentum is typically

[TABLE]

The formed bound state is a rather excited state with a large angular momentum $L$. In order that the bound state annihilates efficiently, it must de-excite into the lower $L$ state. Typically, after a bound state loses $\mathcal{O}(L)$ angular momenta and falls into a deep bound state, the energy gaps between the excited states are larger than the hadronic scale and the bound state can immediately de-excite through spontaneous pion emission. Thus, we below discuss the de-excitation process for the large $L$ bound state to lose $\mathcal{O}(L)$ angular momenta. In terms of the energy, the large $L$ bound state loses $\mathcal{O}(\Lambda_{\text{had}})$.

Ref. [75] estimates the de-excitation rate for the squark by the dipole emission of photons as

[TABLE]

where $Q$ is the electromagnetic charge of the squark, $\alpha$ is the electromagnetic structure constant and $\alpha_{\text{QCD}}(\mu)$ is the QCD structure constant at an energy scale $\mu$. This formula comes from the classical electrodynamics. Since the angular momentum and the quantum number of the bound state is very large, it is reasonable to regard them as a classical bound state. Ref. [80] also estimates a similar formula for the radiative de-excitation 101010Although they do not use higher angular momentum region, as Ref. [75] has discussed, it takes longer for higher angular momentum state to de-excite and that process is more important here..

Also, Refs. [77] and [80] point out the existence of collisional de-excitation process, where electromagnetic or hadronic collision cools the squarkonium down,

[TABLE]

Here, $X$ is $\gamma$ or hadrons. It is expected that these cross sections roughly are equal to the geometrical cross section,

[TABLE]

This cross section is estimated just like the formation cross section 111111Unlike a model discussed in Ref. [76], the quark cloud has the electromagnetic charge and the electromagnetic collision is suppressed only by $\alpha$.. Note that Refs. [77] and [80] discuss that the cross section should be smaller than the destruction cross section by a factor of $\left({\alpha_{\rm QCD}(\Lambda_{\text{had}})E_{X}/\Lambda_{\rm had}}\right)^{2}$, where $E_{X}$ is the energy of $X$. On the other hand, if we apply the naive dimensional analysis [81], these rates may be comparable to each other. This collisional cooling needs $O(L)$ collisions to loose angular momentum enough. In the following we adopt the estimation of Ref. [77].

Just as the de-excitation process, particles in the thermal bath may excite the bound state or even break the bound state. For the transition between the bound states, if the de-excitation rate, Eq. (8), is faster than the cosmic expansion, the bound states are at least in the chemical equilibrium. Since the energy gap between bound states with lower angular momenta is of order $O(E_{B})\gg T_{c}$, where $E_{B}$ is the typical binding energy for the lower angular momentum bound states and satisfies $E_{B}\sim\alpha_{\text{QCD}}^{2}(E_{B})m_{\tilde{q}}$, and this corresponds to the gap of chemical potentials for the bound states, we conclude that the number of the $S$-wave bound state is exponentially larger than the others. Thus, as in Ref. [75], we do not need to include the effect as long as the de-excitation is effective enough.

On the other hand, the destruction of the bound state into SUSY particles should be also considered, as these processes cannot reduce the number of the SUSY particles. If the decay/annihilation rate of the bound state is much greater than the Hubble rate, we can treat this resonance as an intermediate state of the squark annihilation. Then the effective annihilation cross section $\sigma_{\text{form}}^{\prime}$ can be written as

[TABLE]

and the ${\rm BF}((\tilde{q}^{*}\tilde{q})\to{\rm SM^{\prime}s})$ is the effective branching fraction of the squarkonium decays into SM particles and written as

[TABLE]

where $\Gamma_{\text{destruction}}$ is the destruction rate of $\tilde{q}^{*}\tilde{q}$ bound state and $\langle\cdots\rangle$ represents thermal average. The destruction rate is the sum of the decay rate of the squark, $\Gamma_{\text{decay}}$, and the collisional process, $\Gamma_{\text{dest, coll}}$,

[TABLE]

The constituent squark in the bound state can decay into the LSP, if the binding energy is smaller than the mass difference between the squark and LSP. In this estimation, we adopt

[TABLE]

The factor $2$ is for the decay of both squark and anti-squark. In our estimation, in order to take the binding energy into account, we estimate the decay width of the squark mesino by replacing the mass difference, $\Delta m\equiv m_{\tilde{q}}-m_{\tilde{\chi}^{0}_{1}}$, with $\Delta m-\Lambda_{\rm had}$. This replacement can be understood as follows; first, consider the non-relativistic effective theory of the squarkonium bound state, just like Ref. [17]. Then, the width of the bound state comes from the imaginary part of the squark propagator. Only the difference from the squark decay width is the momentum squared, which now corresponds to replacing the mass diffrence as $\Delta m-\Lambda_{\rm had}$. Also, for the beta decay of various atomic nuclei, a similar calculation of the width is known to be consistent with the observation [82].

For the collisional process, the dominant one for $T\gtrsim O(10)\,\text{MeV}$ is the inverse process of the bound state formation $\pi+(\tilde{q}^{*}\tilde{q})\to(\tilde{q}^{*}q)+(\tilde{q}\bar{q})$ although it is suppressed by the Boltzmann factor compared with the formation rate, since we need additional binding energy around $\Lambda_{\text{had}}$ to break the bound state apart into mesinos. The process is regarded as the collision between finite-size objects, just like the bound state formation, and the cross section, $\sigma_{\text{dest, coll, break}}$, is expected to be geometrical: $O(\pi R_{\rm had}^{2})$.

We can also derive the result by applying the Milne relation 121212See Eqs. (19)-(21) in Ref. [42], for example. to Eq. (6); the Milne relation is

[TABLE]

where $p$ is the typical momentum of the mesino for the bound state formation, $\omega$ is the typical energy of the pion for the bound state destruction, being around the hadronic scale, and $g_{i}$ is the degree of freedom for particle $i$. As we have discussed above, the bound state formation process can be decomposed into the partial wave and with the maximal angular momentum for the bound state formation $L_{\text{max}}$, $\sigma_{\text{form}}\sim\pi L_{\text{max}}^{2}/p^{2}$. The resultant squarkonium has an angular momentum up to $L_{\text{max}}$ and thus $g_{\tilde{q}\tilde{q}^{\star}}\sim\sum_{L}^{L_{\text{max}}}(2L+1)\sim L_{\text{max}}^{2}$. The cross section is then

[TABLE]

and indeed geometrical.

However, the relevant number density of the hadron to be multiplied to the cross section in Eq. (16) to obtain the reaction rate includes the Boltzmann suppression factor for the meson to be as energetic as $E\gtrsim\Lambda_{\text{had}}$. Consequently, the destruction rate can be esimated as

[TABLE]

and this destruction rate is suppressed for the temperature much less the hadron scale.

Also, hadronic or electromagnetic collisions can make the conversion of $\tilde{q}+\bar{q}\to\tilde{\chi}^{0}_{1}$ through which a squarkonium can change into a mesino and a LSP. This collision rate is roughly proportional to the temperature $T$. This process may be suppressed if the binding energy is much larger than the squark-DM mass difference.

After the de-excitation, squarks annihilate efficiently. The decay rate is similar to the case of a positronium and given as

[TABLE]

where $\psi$ is the wave function for the squarks in the bound state. For the ground state, this reduces to $\Gamma_{\text{annihilate}}\sim\alpha_{\text{QCD}}^{2}(m_{\tilde{q}})\alpha_{\text{QCD}}^{3}(a_{\rm Bohr}^{-1})m_{\tilde{q}}$, where $\alpha_{\text{QCD}}$ from the wave function is evaluated at the scale of inverse Bohr radius of the bound state, $a_{\rm Bohr}^{-1}\sim\alpha_{\text{QCD}}(a_{\rm Bohr}^{-1})m_{\tilde{q}}$. Since this is much larger than the rate to break such deeply bounded state, which is highly suppressed by the Boltzmann factor, $\exp(-E_{B}/T)$, and the Hubble expansion rate, we can assume that the annihilation is instantaneous once bound states fall into the ground state.

In Fig. 2, we show the effective branching of the squarkonium (13) as a function of the temperature. Here we set $\Lambda_{\rm had}=0.4$ GeV. In Fig. 2(a), we show the case of $\tilde{u}_{R}$ NLSP with the squark-DM mass difference $\Delta m=0,0.5$ and 1 GeV in green, blue and red lines, respectively. For the large mass difference case, the internal conversion $\tilde{q}\to\tilde{\chi}^{0}_{1}$ by either the decay or collision dominates for lower temperature region. In the case of $\Delta m=0$, this internal conversion is prevented due to the squarkonium binding energy. In Fig. 2(b), we show the case of $\tilde{t}_{R}$ with a flavor mixing with the first generation $\delta_{13}=10^{-3}$. With this mixing, the internal conversion is prevented and the effective branching fraction is greater than the $\tilde{u}_{R}$ case.

If the NLSP squark is left-handed, we need to consider the mass difference between the up-type and down-type left-handed squarks. The mass difference comes from the Higgs coupling and the radiative correction. For the Higgs coupling, the mass difference is roughly [83]

[TABLE]

whereas for the radiative correction [84], it is

[TABLE]

Thus, if the squark is heavier than $10\,\text{TeV}$, then both up-type and down-type squarks contribute to the annihilation in the QCD era.

Finally, we briefly comment on the squark-squark bound state, $\tilde{q}\tilde{q}q$, for example. Roughly speaking, the formation cross section should be geometrical, but the $Q$-value is different from the $\tilde{q}^{*}\tilde{q}$ bound state formation process, and thus whether they are formed or not is not obvious [80]. However, once they are formed, the de-excitation and destruction processes will be similar to the $\tilde{q}^{*}\tilde{q}$ state case. As in the case of the $\tilde{q}^{*}\tilde{q}$ bound state, if the de-excitation process continues, the final form of the $\tilde{q}\tilde{q}q$ bound state is a hadron with rotating squark cores. Unlike $\tilde{q}^{*}\tilde{q}$, because of the wavefunction symmetrization, the angular momentum between squarks must change by $2$ during the de-excitation process, but at least by the collisional process can efficiently de-excite them eventually to the $P$-wave ground state. Finally, squarks annihilate exchanging virtual neutralino, whose rate is much larger than the Hubble. Hence, we expect that if they are formed, the annihilation process is not much different from the case of the $\tilde{q}^{*}\tilde{q}$ bound state.

2.2 Chemical Equilibrium

In order that the coannihilation processes reduce the LSP abundance significantly, the LSP and squark NLSP must maintain the chemical equilibrium when the squark effective annihilation is efficient. For higher temperature, the chemical equilibrium is easily maintained. However, in the QCD phase transition era, the temperature of the Universe is rather low. We examine the condition of the chemical equilibrium.

The conversion via hadron collisions (Fig. 3(a)) is the most efficient conversion process. The rate for $\tilde{\chi}^{0}_{1}\to\tilde{q}$ is

[TABLE]

where $m_{h}$ is the incoming hadron mass and $c$ denotes the hypercharge of the squark. We take the pion mass $m_{\pi}$ as $m_{h}$ if the NLSP squark is the first generation, $\tilde{u}$ or $\tilde{d}$. For the other generations, there is another Boltzmann factor such as $\exp(-m_{B}/T)$, as the initial of final state needs to include heavy flavor mesons, if there is no flavor violation in the squark sector. If the mass difference is larger than the meson mass, the rate is Boltzmann suppressed. We also have the weak process or electromagnetic process, but they are generally smaller than Eq. (21), given that the mass difference is the same order as $m_{\pi}$. For Higgsino or wino LSP, we need to use appropriate gauge/Yukawa couplings instead of $c^{2}\alpha$.

Charged flavor violation coupling via the weak interaction (Fig. 3(b)) can be relevant for the heavy flavor scalar quark. This process is suppressed by the CKM mixing and the Fermi constant. However there is no significant Boltzmann suppression due to the heavy flavor meson mass, this process can dominate in low temperature region.

In Fig. 4, we show the conversion rate of $\chi^{0}_{1}\to\tilde{q}$. In Fig. 4(a), we show the cases of right-handed scalar up quark, for the LSP-NLSP mass difference 0, 0.5 and 1 GeV. For larger mass differences, the conversion rate is suppressed by the Boltzmann factor. Still, the chemical equilibrium can be maintained even if the temperature is around 10 MeV. In Fig. 4(b), we show the case of heavy flavor scalar quarks, with minimal flavor violation assumed. Generally speaking, in the SSM, there can be non-minimal flavor violations. If there is a flavor mixing between the first generation and third generation with a mixing parameter $\delta_{13}$, the conversion rate of stop is increased as $\Gamma(\tilde{\chi}^{0}_{1}\to\tilde{t}_{R})\sim\delta_{13}^{2}\Gamma(\tilde{\chi}^{0}_{1}\to\tilde{u}_{R})$. Even if the flavor mixing is tiny, the chemical equilibrium can be realized. For instance, when the mixing parameter $\delta_{13}$ is greater than $O(10^{-6})$, the scalar top can be in chemical equilibrium with the bino LSP for $T\gtrsim 10$ MeV and $\Delta m\lesssim 1$ GeV.

Next, let us comment on the conversion process for the gluino NLSP scenario. If gluino is the NLSP, the virtual squark is involved in the conversion process, as is shown in Fig. 5. The conversion rate between the bino LSP and the gluino NLSP is then

[TABLE]

for the tree level process with zero mass difference, where $c$ and $m_{\tilde{q}}$ are the hypercharge and mass of the lightest squark, respectively, and $m_{h}$ is the mass of hadron involved in the squark-quark-gluino vertex. Requiring that the conversion rate is greater than the Hubble rate at least around $T\sim T_{c}$, we obtain

[TABLE]

for the first generation squark mass. Since this must be heaver than the LSP and the NLSP, we conclude that the gluino must be lighter than TeV for the bino LSP and the gluino NLSP scenario. Also, even if the LSP is wino, the mass bound is small as well. For such cases, the perturbative annihilation is effective enough to reduce the abundance, as discussed in Sec. 2.1.1.

On the other hand, if the Higgsino is the LSP and the gluino is the NLSP, the loop-level conversion rate is rather large, as the gluino-Higgsino-gluon dipole operator is relatively enhanced [85, 86, 87, *Sato:2013bta]. The dominant contribution comes from the stop-top loop and the conversion rate is roughly given by

[TABLE]

Then, comparing with the Hubble expansion rate at $T=T_{c}$, we obtain

[TABLE]

Thus, it may be possible to exceed the perturbative annihilation for this case. We will estimate the abundance in future [89].

2.3 Relic Abundance

Now, we calculate the relic abundance, combining the dynamics discussed above. First, we compute the yield for superparticles based on Sec. 2.1.1. Then, we need to keep track of the non-perturbative annihilation process described in Sec. 2.1.2.

For the non-perturbative annihilation process, we need to solve the Boltzmann equation for coannihilation between the LSP and super-hadrons. The Boltzmann equations are given as:

[TABLE]

where BF stands for the thermal averaged effective branching fraction of the $\tilde{q}^{*}\tilde{q}$ bound state 131313Note that we assume that the annihilation between super-hadrons also occurs in the same rate as the annihilation between super-hadron and anti-super-hadron, as discussed in Sec. 2.1.2.. We assume the non-perturbative annihilation starts at the QCD phase transition temperature $T=T_{c}$. In the following, we take this temperature 200 MeV. In the QCD confined phase, the degrees of freedom of squark hadrons is non-trivial. In our estimation, we simply assume that this value is 6 as in the case of a free squark. In this calculation, we have used values of $g_{*s}$ and $g_{*}$ estimated in Ref. [90]. The initial yield is the perturbative value derived in Sec. 2.1.1. As we have discussed in the previous subsection, the conversion process and thus the end of the integration highly depends on the squark flavor structure and mass difference and in some case fails to keep the chemical equilibrium.

Finally, we show our result. First, the temperature dependence of the yield for the bino and $\tilde{u}_{R}$ is shown in Fig. 6, where $m_{\rm DM}=10^{6}$ GeV and we take $\Lambda_{\text{had}}=400$ MeV. We show the cases of the mass difference $\Delta m=0$ and 0.2 GeV in red and blue lines, respectively. We see that, at $T\sim m_{\rm DM}/30$, the first freeze-out occurs. When $T=200$ MeV, the non-perturbative annihilation starts and the abundance rapidly reduces. In the case that $\Delta m=0$ GeV, at $T\sim 10$ MeV, the annihilation happens further again, since the effective annihilation rate is enhanced as seen in Fig. 2.

In Fig. 7, we show the contour plot of $\Omega_{\rm{DM}}h^{2}=0.12$ [91] on $m_{\rm DM}-\Delta m$ plane. We assume the LSP is the bino. The blue, red and green lines show the cases of $\tilde{u}_{R}$, $\tilde{b}_{R}$ and $\tilde{t}_{R}$ NLSP, respectively. Without flavor violation, the bino and $\tilde{t}_{R}$ are not in chemical equilibrium and we assume a small mixing between the first generation $\delta_{13}=10^{-8}$, which leads to partial chemical equilibrium. A small flavor violation $\delta_{13}$ larger than $O(10^{-6})$ can keep almost complete chemical equilibrium. Moreover the branching fraction of the resonance $(\tilde{q}^{*}\tilde{q})$ into SM particles also increases (Fig. 2) accordingly the dark matter abundance is smaller than the $\tilde{u}_{R}$ case. For the $\tilde{b}_{R}$ case, the chemical decoupling occurs at $T\sim 100$ MeV and $\tilde{b}_{R}$ cannot annihilate at $T\sim 10$ MeV, where the annihilation cross section is enhanced. Thus at small mass differences, the abundance for the $\tilde{b}_{R}$ case is larger than the $\tilde{u}_{R}$ case.

In the discussions above, we mainly consider the bino LSP. For the wino and Higgsino LSP, essentially identical estimation is possible. The DM abundance of these cases is slightly larger than the bino case, due to the large number of degrees of freedom.

Now let us discuss the uncertainties of the present estimation. As we focus on the non-perturbative QCD phase, there are so many $O(1)$ uncertainties in every stage of the estimation. One of the largest uncertainties of our calculation is the de-excitation rate of the $\tilde{q}^{*}\tilde{q}$ bound state. In this analysis, we adopt the estimation based on Refs. [77] and [80], which indicates $O(10^{-2})$ suppression, compared to the estimation by the naive dimensional analysis. If there is no such suppression, the DM abundance can be further reduced.

3 Direct Detection

In this section, we discuss the direct detection of the dark matter. Since the squark and neutralino LSP are degenerated in mass, the dark matter scattering cross section to the nucleon is enhanced [93, 21, 94]. The Feynman diagram is shown in Fig. (8). First, we consider bino LSP and right-handed squark NLSP. Then, we briefly discuss the other cases.

The effective Lagrangian for spin-independent scattering between the bino and the quark by the squark exchange is [95]

[TABLE]

where $\tilde{B}$ is the bino fermion field and $\mathcal{O}^{q}_{\mu\nu}\equiv\frac{1}{2}\bar{q}\left(D_{\mu}\gamma_{\nu}+D_{\nu}\gamma_{\mu}-\frac{1}{2}g_{\mu\nu}\not{D}\right)q$. The coefficients $f_{q}$ and $g_{q}^{(1,2)}$ are given as

[TABLE]

and $g_{q}^{(2)}=0$. In the present setup, the mixing between the right-handed and left-handed squarks is negligible.

The spin-independent scattering cross section of the bino with a nucleon (atomic number $Z$ and mass number $A$) is given by

[TABLE]

where $M_{\text{red}}$ is the reduced mass of the bino and nucleon system.

[TABLE]

for a nucleon $N$. The form factor $f_{Tq}^{N}$ is given by

[TABLE]

and $q(2)$ and $\bar{q}(2)$ are the second moments of the parton distribution function.

Using the lattice calculation for $f_{Tq}^{N}$ [96] and the CTEQ parton distribution function [97] at the $Z$ Boson mass scale, we obtain

[TABLE]

for the nucleon $N$ and squark $\tilde{q}_{i}=\tilde{u},\tilde{d}$ and $\tilde{s}$. For the heavier flavor $\tilde{q}_{i}=\tilde{c}$ and $\tilde{b}$, assuming $\Delta m\sim m_{q_{i}}$,

[TABLE]

We summarize the coefficient $s_{Ni}$ in Tab. 1. Note that, however, if the mass difference is too small, the estimation at tree-level is not accurate [93, 21, 94]. Moreover in such a case, the uncertainty from the non-perturbative QCD effects will be also sizable.

In Fig. 7, we show the constraint for the case of $\tilde{u}_{R}$ NLSP. We see that even 100 TeV DM can be excluded in this case. In fact, the light flavor squark coannihilation cases are also disfavored by the direct detection constraint. This constraint can be, however, evaded by introducing non-minimal flavor violation. For instance if the $\tilde{t}_{R}$ is the NLSP and there is a very small flavor mixing with the first generation, e.g., $\delta_{13}=10^{-4}$, we can still keep the successful coannihilation in the QCD era, and avoid the DM direct detection experiments as the cross section is suppressed by a factor of $(\delta_{13})^{4}$. A similar estimation also applies to the wino and Higgsino DM. In the case of the Higgsino DM, the direct detection rate is suppressed by the Yukawa couplings of the quarks and the constraint of the direct detection is also relaxed.

Although the DM detection cross section with very degenerated squark can have large uncertainty, we see that anyway if the mass difference is smaller, the cross section is significantly increased. The direct detection constraints will disfavor several combinations of the LSP species and squark flavors.

4 Conclusion and Discussion

In this paper, we investigate the possibility that superheavy neutralino dark matter realizes a viable thermal relic density. If the mass difference between the colored NLSP and the neutralino LSP is smaller than $O(1)$ GeV, non-perturbative QCD effects drastically enhance the NLSP annihilation cross section, and further reduce the DM density. Although there are large uncertainties in the estimation of the non-perturbative effects, we found that a sub-PeV scale neutralino DM can be consistent with the current observation. Compared to the previous study for the relic abundance estimation of the colored heavy particles, the destruction of the resonance of two-body bound state of the colored heavy particles is increased and then the effective annihilation cross section of the colored particles is reduced. Moreover, we study the chemical equilibrium condition between the LSP and NLSP and estimate the effect of out-of-chemical equilibrium on the DM abundance.

Such a highly mass degeneration between the LSP and colored NLSP increases the DM direct detection cross section. We see that several cases of DM nature and flavor of the squark NLSP are already excluded by the current DM direct detection experiments. In future direct detection experiments, much larger parameter space, typically ten times larger mass difference than the current constraint can be tested.

Now let us comment on the gluino coannihilation case. In principle, the gluino annihilation rate can be increased after the QCD phase transition era. However, as discussed before, chemical equilibrium of the gluino and neutralino is prevented if the squark is too heavy. This constraint implies that the upper-bound of the squark is $O(1-10)$ TeV, depending on the type of the LSP neutralino. In the gluino coannihilation case, the gluino mass should be less than the squark mass and accordingly the upper-bound of the gluino is also $O(1-10)$ TeV. We leave the detailed study in a future work [89].

In this paper, we consider the case that the LSP and NLSP mass is extremely close to each other. Such a fine-tuned mass spectrum may not be natural. However, it is very important to clarify what is the maximal DM mass in the SSM, for comprehensive searches of SSMs.

Acknowledgments

S.S. thanks to A. Kamada and S. Matsumoto for useful discussions. This work is supported by Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan, No. 17H02878 and 18K13535 (S.S.) and by World Premier International Research Center Initiative (WPI), MEXT, Japan. The work of H.F. is supported in part by a Research Fellowship for Young Scientists from the Japan Society for the Promotion of Science (JSPS). F.L. is supported by the One Hundred Talent Program of Sun Yat-sen University, China.

Note Added: During the completion of this paper, a related study [98] was submitted to arXiv which discusses the non-perturbative effect of the colored coannihilation.