How heavy can neutralino dark matter be?

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.


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.

JHEP04(2019)107
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][46]. 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 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 ∼ 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 nonstandard-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.

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.

Squark annihilation
If the mass difference between the neutralino LSP and squark NLSP is less than ∼ 1 GeV, the cosmological evolution of the dark matter shows interesting features. Once the temperature gets lower than around m 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 JHEP04(2019)107 process freezes out again. In the following, we describe these two perturbative and nonperturbative freeze-outs. 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, nq/nχ0 1 = n eq q /n eq χ 0 1 . 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, 2Ỹ ≡ (nq + nq * + nχ0 G N is the gravitational constant, mχ0 1 is the neutralino LSP mass, g * s and g * 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, σ eff v , should include contributions of all (co)annihilation channels which can changeỸ . However, for the case of highly degenerated colored NLSP and neutralino LSP we focus in this paper, the dominant contributions to σ eff v come fromqq * 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,qq * → gg andqq * → qq, in calculating σ eff v . These two processes are universal for all squark flavors and are usually the dominant ones. 3 We also include in σ eff v theqq * Sommerfeld and bound-state effects as discussed in detail 1 To 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 We assume that there is no asymmetry betweenq andq * , and we take nq = nq * in our calculation. 3 Ref. [50] discussed the case when the dominant annihilation channels arett * → hh, W + W − , ZZ, when the interactions between stop and Goldstone bosons are large.

JHEP04(2019)107
in ref. [44]. 4 With these considerations, we use where gχ0 1 = 2 is the degrees of freedom for a neutralino, assumed to be a bino. gqq * = 6 is the sum of the degrees of freedom for squark and antisquark. σv qq * is the thermallyaveraged annihilation cross section times relative velocity, 5 given as In the above expression, the first term is the sum of the Sommerfeld enhancedqq * → gg andqq * → qq cross sections, and the second term is the contribution from theqq * boundstate effect, 6 where σv bsf is the bound state formation cross section, Γ η is the bound state annihilation decay rate, and Γ 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 refs. [44,72]. We evaluate eq. (2.1) from T ∼ mχ0 1 to lower temperatures. The usual freeze-out (whenỸ /Ỹ eq − 1 becomes larger than 1) happens around T ∼ mχ0 1 /30. After that,Ỹ continues decreasing, and in particular the bound-state effect becomes effective when T is around and below the bound-state binding energy, E B ∼ mχ0 1 /100.Ỹ essentially ceases to decrease after T ∼ mχ0 1 /1000, until the new annihilation process kicks in after the QCD phase transition which we describe in the next subsection.

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]. 7 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 4 The Sommerfeld and/or bound-state effects are also recently discussed in . 5 As is noted in [44], since gqq * = 6 is the sum of the degrees of freedom for squark and antisquark, a factor of 1/2 needs to be put in σv qq * , if using the usual spin and color averaged cross section. This is explained in the appendix of [73]. 6 We 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. 7 Also, see ref. [76] for the analytic calculation, which is basically consistent with ref. [75]. becomes a SUSY hadron such as a mesino. The SUSY hadron has a large geometrical size Λ −1 QCD and large scattering cross section with another SUSY hadron. Once such a scattering occurs, a squark-squark bound stateq * 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 theq * 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 figure 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 theq * q bound state formation. Soon after the QCD phase transition at temperature T c 200 MeV, the squarks start to form mesinos, for example,qq andq * q. Since the squark is much heavier than the QCD scale, the squark is localized at the center in the mesino bound state, ∼ m −1 q in radius, and the quark/gluon partons surround it with a radius of R had = Λ −1 had , where Λ had is the hadronic mass scale. Then, these heavy hadrons interact with each other, forming a squarkoniumq * q , a squark-antisquark bound state, where the cross section is the geometrical one, 8 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 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 max ≡ pR had , where p is the squark momentum. Summing up the partial wave cross section up to L max , we obtain the geometrical cross section [78]; Ref. [76] also discusses this behavior analytically. Then let us estimate R had = Λ −1 had . In ref. [75], Λ had is taken to be 1 GeV. On the other hand, the pion charge radius is experimentally 0.67 fm ∼ (300 MeV) −1 [79], which JHEP04(2019)107 suggests to use Λ had ∼ 300 MeV. Also, ref. [77] uses a slightly larger value, R had ∼ 0.5 fm ∼ (400 MeV) −1 . Moreover, ref. [80] also claims that such larger value can be available. 9 Thus, combining these contexts, we take the Λ 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 ∼ 2mqT , and the impact parameter is r ∼ Λ −1 had , the angular momentum is typically 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 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 O(L) angular momenta. In terms of the energy, the large L bound state loses O(Λ had ).
Ref. [75] estimates the de-excitation rate for the squark by the dipole emission of photons as where Q is the electromagnetic charge of the squark, α is the electromagnetic structure constant and α QCD (µ) is the QCD structure constant at an energy scale µ. 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. 10 Also, refs. [77] and [80] point out the existence of collisional de-excitation process, where electromagnetic or hadronic collision cools the squarkonium down, Here, X is γ or hadrons. It is expected that these cross sections roughly are equal to the geometrical cross section, This cross section is estimated just like the formation cross section. 11 Note that refs. [77] and [80] discuss that the cross section should be smaller than the destruction cross section 9 They 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. 10 Although 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. 11 Unlike a model discussed in ref. [76], the quark cloud has the electromagnetic charge and the electromagnetic collision is suppressed only by α.

JHEP04(2019)107
by a factor of (α QCD (Λ had )E X /Λ had ) 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 deexcitation rate, eq. (2.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 ) T c , where E B is the typical binding energy for the lower angular momentum bound states and satisfies E B ∼ α 2 QCD (E B )mq, 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 σ form can be written as 12) and the BF((q * q ) → SM s) is the effective branching fraction of the squarkonium decays into SM particles and written as where Γ destruction is the destruction rate ofq * q bound state and · · · represents thermal average. The destruction rate is the sum of the decay rate of the squark, Γ decay , and the collisional process, Γ dest, coll , Γ destruction = Γ decay + Γ dest, coll . (2.14) 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 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, ∆m ≡ mq − mχ0 1 , with ∆m − Λ 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 JHEP04(2019)107 as ∆m − Λ 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 O(10) MeV is the inverse process of the bound state formation π + (q * q ) → (q * q) + (qq) although it is suppressed by the Boltzmann factor compared with the formation rate, since we need additional binding energy around Λ 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, σ dest, coll, break , is expected to be geometrical: O(πR 2 had ). We can also derive the result by applying the Milne relation 12 to eq. (2.6); the Milne relation is where p is the typical momentum of the mesino for the bound state formation, ω 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 max , σ form ∼ πL 2 max /p 2 . The resultant squarkonium has an angular momentum up to L max and thus gqq ∼ Lmax L (2L+1) ∼ L 2 max . The cross section is then and indeed geometrical. However, the relevant number density of the hadron to be multiplied to the cross section in eq. (2.16) to obtain the reaction rate includes the Boltzmann suppression factor for the meson to be as energetic as E Λ had . Consequently, the destruction rate can be esimated as Γ dest, coll σ dest, coll, break × n π,E Λ had ∼ πR 2 had Λ 2 had T exp(−Λ had /T ), (2.17) and this destruction rate is suppressed for the temperature much less the hadron scale. Also, hadronic or electromagnetic collisions can make the conversion ofq +q →χ 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 where ψ is the wave function for the squarks in the bound state. For the ground state, this reduces to Γ annihilate ∼ α 2 QCD (mq)α 3 QCD (a −1 Bohr )mq, where α QCD from the wave function is  Bohr ∼ α QCD (a −1 Bohr )mq. 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 figure 2, we show the effective branching of the squarkonium (2.13) as a function of the temperature. Here we set Λ had = 0.4 GeV. In figure 2a, we show the case ofũ R NLSP with the squark-DM mass difference ∆m = 0, 0.5 and 1 GeV in green, blue and red lines, respectively. For the large mass difference case, the internal conversionq →χ 0 1 by either the decay or collision dominates for lower temperature region. In the case of ∆m = 0, this internal conversion is prevented due to the squarkonium binding energy. In figure 2b, we show the case oft R with a flavor mixing with the first generation δ 13 = 10 −3 . With this mixing, the internal conversion is prevented and the effective branching fraction is greater than theũ 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] whereas for the radiative correction [84], it is ∆m ∼ 70 MeV. (2.20) Thus, if the squark is heavier than 10 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,qqq, for example. Roughly speaking, the formation cross section should be geometrical, but the Q-value JHEP04(2019)107 is different from theq * 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 theq * q state case. As in the case of theq * q bound state, if the de-excitation process continues, the final form of theqqq bound state is a hadron with rotating squark cores. Unlikeq * 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 Pwave 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 theq * q bound state.

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 (figure 3a) is the most efficient conversion process. The rate forχ 0 1 →q is where m h is the incoming hadron mass and c denotes the hypercharge of the squark. We take the pion mass m π as m h if the NLSP squark is the first generation,ũ ord. 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. (2.21), given that the mass difference is the same order as m π . For Higgsino or wino LSP, we need to use appropriate gauge/Yukawa couplings instead of c 2 α.  Charged flavor violation coupling via the weak interaction (figure 3b) 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 figure 4, we show the conversion rate of χ 0 1 →q. In figure 4a, 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 figure 4b, 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 δ 13 , the conversion rate of stop is increased as Γ(χ 0 1 →t R ) ∼ δ 2 13 Γ(χ 0 1 →ũ R ). Even if the flavor mixing is tiny, the chemical equilibrium can be realized. For instance, when the mixing parameter δ 13 is greater than O(10 −6 ), the scalar top can be in chemical equilibrium with the bino LSP for T 10 MeV and ∆m 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 figure 5.

JHEP04(2019)107
The conversion rate between the bino LSP and the gluino NLSP is then for the tree level process with zero mass difference, where c and mq 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 ∼ T c , we obtain 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 section 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][88]. The dominant contribution comes from the stop-top loop and the conversion rate is roughly given by Then, comparing with the Hubble expansion rate at T = T c , we obtain mt < O(10) TeV.
(2.25) Thus, it may be possible to exceed the perturbative annihilation for this case. We will estimate the abundance in future [89].

Relic abundance
Now, we calculate the relic abundance, combining the dynamics discussed above. First, we compute the yield for superparticles based on section 2.1.1. Then, we need to keep track of the non-perturbative annihilation process described in section 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: where BF stands for the thermal averaged effective branching fraction of theq * q bound state. 13 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 section 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ũ R is shown in figure 6, where m DM = 10 6 GeV and we take Λ had = 400 MeV. We show the cases of the mass difference ∆m = 0 and 0.2 GeV in red and blue lines, respectively. We see that, at T ∼ m 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 ∆m = 0 GeV, at T ∼ 10 MeV, the annihilation happens further again, since the effective annihilation rate is enhanced as seen in figure 2.
In figure 7, we show the contour plot of Ω DM h 2 = 0.12 [91] on m DM − ∆m plane. We assume the LSP is the bino. The blue, red and green lines show the cases ofũ R ,b R andt R NLSP, respectively. Without flavor violation, the bino andt R are not in chemical equilibrium and we assume a small mixing between the first generation δ 13 = 10 −8 , which leads to partial chemical equilibrium. We also plot in dotted purple line the current constraint on the dark matter direct detection with XENON1T [92], for the case of the bino LSP andũ R NLSP. The region below this line is disfavored.
resonance (q * q ) into SM particles also increases (figure 2) accordingly the dark matter abundance is smaller than theũ R case. For theb R case, the chemical decoupling occurs at T ∼ 100 MeV andb R cannot annihilate at T ∼ 10 MeV, where the annihilation cross section is enhanced. Thus at small mass differences, the abundance for theb R case is larger than theũ 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 theq * 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.

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 [21,93,94]. The Feynman diagram is shown in figure 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] L eff = f q m qBBq q + g whereB is the bino fermion field and O q µν ≡ 1 2q D µ γ ν + D ν γ µ − 1 2 g µν / D q. The coefficients f q and g (1,2) q are given as and g (2) q = 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 where M red is the reduced mass of the bino and nucleon system.
for a nucleon N . The form factor f N T q is given by and q(2) andq(2) are the second moments of the parton distribution function. Using the lattice calculation for f N T q [96] and the CTEQ parton distribution function [97] at the Z Boson mass scale, we obtain   for the nucleon N and squarkq i =ũ,d ands. For the heavier flavorq i =c andb, assuming ∆m ∼ m q i , We summarize the coefficient s N i in table 1. Note that, however, if the mass difference is too small, the estimation at tree-level is not accurate [21,93,94]. Moreover in such a case, the uncertainty from the non-perturbative QCD effects will be also sizable. In figure 7, we show the constraint for the case ofũ 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 thet R is the NLSP and there is a very small flavor mixing with the first generation, e.g., δ 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 (δ 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.

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 JHEP04(2019)107 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.