Thermal Relic Dark Matter Beyond the Unitarity Limit

We discuss a simple model of thermal relic dark matter whose mass can be much larger than the so-called unitarity limit on the mass of point-like particle dark matter. The model consists of new strong dynamics with one flavor of fermions in the fundamental representation which is much heavier than the dynamical scale of the new strong dynamics. Dark matter is identified with the lightest baryonic hadron of the new dynamics. The baryonic hadrons annihilate into the mesonic hadrons of the new strong dynamics when they have large radii. Resultantly, thermal relic dark matter with a mass in the PeV range is possible.


I. INTRODUCTION
Despite overwhelming evidence of the existence of dark matter, its identity has remained unknown for almost eighty years since its first postulation. We are only almost certain that dark matter is not a part of the standard model of the elementary particle physics.
Therefore, it is one of the most important tasks of modern particle physics to identify the origin of dark matter (see e.g. [1][2][3]).
Among various candidates for dark matter, thermal relic dark matter is one of the most attractive candidates [4][5][6][7][8][9]. The thermal relic dark matter explains the observed dark matter density by its freeze-out from the thermal bath. For the s-wave annihilation, for example, the observed dark matter density is reproduced when the annihilation cross section satisfies σv 3 × 10 −26 cm 3 /s. The beauty of thermal relic dark matter is that the resultant density does not depend on the initial condition as long as dark matter was in the thermal equilibrium in the early universe.
As an important consequence of thermal relic dark matter, there is an upper limit on the mass of dark matter from the so-called unitarity limit on the annihilation cross section [10].
In fact, the s-wave annihilation cross section of dark matter with a mass M is limited from above by unitarity; Combined with the required cross section mentioned above, the upper limit on the dark matter mass turns out to be about a hundred TeV.
In this paper, we challenge the unitarity limit on the mass of thermal relic dark matter.
In fact, the above unitarity limit applies when the dark matter is a point-like particle. If dark matter is a bound state with a large radius compared with its Compton length, on the other hand, it may have a geometrical cross section for annihilation [10] (see also [11][12][13]).
With the larger cross section, thermal relic dark matter with a mass much larger than a few hundred TeV is possible. We construct a simple model where bound state dark matter annihilate while they have large radii and hence have a large geometrical cross section.
This mechanism should be compared with the enhancement of the dark matter annihilation cross section by the so-called Sommerfeld enhancement [14][15][16][17]. In this case, the dark matter itself is a rather point-like particle, and hence, the enhanced cross section does satisfy the unitarity limit of point-like particles (see Ref. [18] for recent discussion). 1 In a model discussed in this paper, on the other hand, dark matter itself is a bound state and has an annihilation cross section of a geometrical size with which the number density of dark matter is significantly reduced.
The organization of this paper is as follows. In section II, we introduce a model based on a simple strongly coupled gauge theory. In section III, we discuss thermal history and the relic density of dark matter. The final section is devoted to conclusions and discussions.
There, we also comment on a possible application of the present model to explain the excess of the observed flux of extraterrestrial neutrinos in the PeV range [21][22][23]. In the appendices, we also discuss two alternative models.

II. MODEL OF DARK MATTER WITH AXION PORTAL
Let us consider an SU (N c ) gauge theory with one-flavor of Weyl fermions, (U ,Ū ), in the fundamental and the anti-fundamental representations. We call (U,Ū ) the quarks in the following. The quark does not carry any gauge charges under the Standard Model gauge groups. For a while, we assume that the quark possesses a mass, M U .
As a special feature of the present model, we arrange the dynamical scale of SU (N c ), Λ dyn , to be much smaller than M U . That is, we take the gauge coupling constant at the renormalization scale around M U small; where α Nc = g 2 Nc /4π is the fine-structure constant. Below M U , the model behaves as the pure-Yang Mills theory. According to the standard understanding of QCD, this theory also exhibits confinement, which has been confirmed by lattice simulations e.g. [24][25][26] (see also 1 The same is true in the models with the so-called Breit-Wigner enhancement [19,20]. be stretched much longer than Λ −1 dyn due to the heaviness of the quarks [31]. It eventually breaks-up and creates a pair of a quark and an anti-quark when its length becomes of O(M U /F Nc ) where F Nc denotes the string tension made by the flux tube. Therefore, the SU (N c ) gauge dynamics leads to a rather long-range force even after confinement.
The quarks are stable and can be a dark matter candidate due to a vector-like global U (1) symmetry under which the quarks are charged. We call this symmetry the U (1) B symmetry.
The quarks, however, do not become dark matter as they are. As noted above, they are confined into hadrons when the temperature of the universe becomes lower than the critical temperature T c = O(Λ dyn ). Below the critical temperature, the U (1) B charges of the quarks are inherited to the baryons, and the lightest baryon, becomes dark matter eventually. 2 The mesons, on the other hand, do not carry the U (1) B charges and are not stable. In fact, the ground state meson, for example, immediately decays into a pair of the glueballs as we will see shortly.
For a successful model of thermal relic dark matter, the above dark matter sector needs to be connected to the Standard Model. As an example of such connection, we here consider a model with "axion portal". 3 For that purpose, we first replace the mass term of the quark with an interaction term to a singlet complex scalar field φ and assume that the model possesses an approximate chiral symmetry, U (1) A . Here, g denotes a coupling constant of O(1). The quark obtains a mass M U = g φ when the U (1) A chiral symmetry is spontaneously broken by a vacuum expectation value (VEV) of φ.
At around the VEV of φ, φ = f a / √ 2, φ is decomposed into a scalar boson ρ and a pseudo Nambu-Goldstone boson a, 2 The lightest baryon, B 0 , possesses a spin N c /2 due to the fermi-statistics. 3 In the appendices A and B, we discuss models with "higgs portal" and "hypercharged particle" to the Standard Model sector as alternative examples. We may also consider models with a "vector portal" in which a dark photon connects the two sectors.
The mass of the scalar boson ρ is expected to be of O(f a ). As we will see shortly, however, the mass of ρ should be somewhat suppressed for a successful model. The "axion" component a, on the other hand, obtains a mass from explicit breaking of the U (1) A symmetry. When the explicit breaking effects are dominated by the U (1) A anomaly of SU (N c ), the axion mass is estimated to be which is much smaller than the dynamical scale. Similarly to (U,Ū ), the newly introduced (d ,d ) also couples to φ via, After integrating out (d ,d ), we obtain effective interactions of the axion to the Standard Model gauge bosons, where α QCD and α QED are the fine-structure constants of QCD and QED, respectively. The Lorentz indices of the field strengths G (QCD) and F (QED) should be understood. Now, we have all the necessary components of the model of dark matter. The relevant features for the following arguments are; • SU (N c ) gauge theory with one-flavor of quarks, (U,Ū ), whose mass is much larger than the dynamical scale (M U Λ dyn ).
• The mass of (U,Ū ) is generated as a result of spontaneous breaking of an approximate U (1) A chiral symmetry, i.e. M U = g φ .
• The axion associated with spontaneous breaking of an approximate chiral symmetry couples to both the dark matter sector and the Standard Model sector. • The U (1) B charge of the quarks are inherited to the baryons after the confinement.
• The mesons decay immediately into glueballs and axions.
• The glueballs decay into the axions which eventually decay into the Standard Model gauge bosons.
The scalar boson ρ immediately decays into a pair of axions, and hence, it does not play a crucial role in the following discussion.
Before closing this section, let us give a rough sketch of thermal history which will be discussed in the next section. When the temperature of the universe is much higher than M U , the quarks are in the thermal bath. Once the temperature becomes lower than M U , the annihilation process freezes-out and the resultant relic density per the entropy density s is given by [9], Here, T F denotes the freeze-out temperature, x the temperature mass ratio, x = M U /T , g * (T ) the massless degrees of freedom at T , and M PL 2.4 × 10 18 GeV the reduced Planck scale. The freeze-out temperature is recursively determined by, where g U denotes the degree of the freedom of U , i.e. g U = 4N c . A typical freeze-out temperature is given by x F ∼ O(10).
At around the freeze-out temperature, the quarks mainly annihilate into φ's, with the spin and color averaged annihilation cross section, where α 2 g = g 2 /(4π). We neglect the annihilation into a pair of the gluons due to Eq. (2). Below the freeze-out temperature, the number density of the quarks are diluted by cosmic expansion, and a typical distance between the quarks at a temperature T is given by, When the temperature decreases to the critical T c Λ dyn , the SU (N c ) gauge interaction becomes strong and exhibits confinement. Below this temperature, the quarks do not freely fall separately anymore. In the following, we discuss the fates of the bound states assuming that phase transition is first order according to Refs. [32,33]. 5

B. Bound State Formation
In order to trace the thermal history below the dynamical scale precisely, we need to solve the strong gauge dynamics, which is impossible with the current techniques. Here, instead, we follow the picture in Ref. [12], and treat hadrons as composites of heavy quarks which are attracted with each other by a phenomenological potential (see e.g. [34]), Here, κ is an O(1) numerical factor that depends on the color exchanged between the quarks. For a color singlet configuration of a quark and an anti-quark, for example, κ = The linear term represents the effects of non-perturbative dynamics and F Nc corresponds to the tension of the flux tube. At a high temperature, F Nc (T ) is vanishing while F Nc ∼ Λ 2 dyn below the critical temperature T c = O(Λ dyn ). 6 The gauge coupling constant α Nc in Eq. (15) is, on the other hand, estimated at the renormalization 5 The following arguments are not altered significantly as long as the growth of the string tension of the strong dynamics is fast enough. 6 The lattice simulations suggest T c / F Nc 0.6 for the pure Yang-Mills SU (N c ) (N c ≥ 3) theories [32]. scale corresponding to the Bohr radius µ κα Nc (µ)M U as the leading order approximation (see Fig. 2).
When the temperature of the universe becomes lower than T c , the SU (N c ) gauge dynamics transits into the confined phase and the quarks and gulons are confined into color singlets.
In particular, the quarks at the distance D(T c ) in Eq. (12) are pulled with each other by the linear potential, and the sizes of the quark bound states become much shorter than the original distance. 7 It should be noted that the quarks are not accelerated even when they are pulled by the strong force due to frictions caused by the interactions with the glueballs in the thermal bath.
To estimate the typical size of the bound state at a temperature, T , let us consider a partition function of a quark and anti-quark bound state by the potential in Eq. (13); Here, the reduced mass of the two body system is given by M U /2. For the negative energy states where the Coulomb potential is dominant, i.e. r < (κα Nc /F Nc ) 1/2 , we approximate their energy eigenvalues by 7 In the parameter space we are interested in, D(T c ) is shorter than the length of the string breaking, M U /F Nc , on the other hand, the strings between the quarks break up immediately and the quarks are dominantly confined not into baryons but into mesons especially for large N c . In this situation, the relic abundance of the baryon dark matter can be much smaller than the present scenario, which will be discussed elsewhere. In the blue shaded band we vary n max from one to three times of the one defined by r nmax = (κα Nc /F Nc ) 1/2 . The horizontal red line corresponds to the Bohr radius. (Right) The fractional occupation numbers of the negative energy state, ξ(E < 0), and the ground state, ξ 1 . Here, we fix n max to be the one defined by r nmax = (κα Nc /F Nc ) 1/2 . In both panels, we fix F Nc Λ 2 dyn even for T > T c Λ dyn for presentation purpose.
Here, n denotes the principal quantum number and the radii of the corresponding states are given by, For the positive energy states which correspond to r > (κα Nc /F Nc ) 1/2 , on the other hand, we approximate them by continuous spectrum (see Fig. 3). We checked that the above approximation well reproduces a quantum statistical partition function with approximate energy eigenvalues in Ref. [35]. For ease of the computation, we rely on the approximation in Eq. (14) in the following arguments.
In Fig. 4, we show the typical size of the quark bound state for a given temperature estimated by We also show the fractional occupation numbers of the negative energy state, ξ(E < 0), and the ground state, ξ 1 , respectively. Here, n max is defined by r nmax = (κα Nc /F Nc ) 1/2 , although the results do not depend on the precise value of n max significantly. The figure shows that R(T c ) = O(Λ −1 dyn ). Thus, we find that the bound states are in excited states below the critical temperature.
When the temperature decreases further, the bound states are de-excited and the typical size becomes r 1 in Eq. (16).
It should be noted that quarks in the ground state are knocked out to the excited states by scatterings with the glueballs in the thermal bath. The rate of such processes is roughly Here m S denotes the glueball mass which is slightly larger than the scale of the string tension in pure Yang-Mills theories [24][25][26]. In the parameter region we are interested in, Γ ex is larger than the Hubble expansion rate at T T c . Therefore, the each bound state transits between the ground state to the excited states rather frequently (see Fig. 5). This behavior plays a crucial role for the final dark matter abundance.

C. Fate of Mesons
As we have seen above (e.g. Fig. 4), the bound states shrink and get de-excited to the ground state once the temperature of the universe becomes much lower than T T c . Once the bound states stay in the ground state, they immediately decay into the glueballs and the scalars φ (i.e. a's and ρ's) in which the heavy quarks annihilate microscopically (Fig. 6).
The decay rate is given by the annihilation rate multiplied by the radial wave function of the ground state at around the origin, 9 Since this rate is much larger than the Hubble expansion rate, the mesons decay away very quickly. It should be also noted that the bound states spend a small fraction of their time as the ground state even around T T c . Thus, the mesons start to decay without waiting for complete de-excitation, as long as Γ M 0 × ξ 1 is larger than the Hubble expansion rate. As a result, we find that the mesons decay away from the thermal bath immediately for T T c .
Excited glueball states decay into lower-lying states immediately. 10 The ground state CP -even glueball, S 0 , decays into a pair of the axions through the mixing to ρ (see Fig. 7). 9 The Bohr radius is of the order of (α Nc M U ) −1 . 10 The masses of some low-lying states may be smaller that the twice of the mass of the ground state glueball.
Those states decay by emitting off-shell glueballs and have decay rates similar to the one of the ground state.
The CP -odd glueball decays into a pair of S 0 and an axion with a much higher rate. The decay rate of the CP -even ground state glueball is roughly estimated by, Here, the mixing angle between ρ and S 0 is estimated to be, based on the Naive Dimensional Analysis [36,37]. In terms of the cosmic temperature, the decay temperature of the glueball is roughly given by, Thus, the glueballs also decay away immediately unless ρ is very much heavier than m S .
The massive glueballs decouple from the thermal bath when their annihilation into the axions freeze-out, which leaves the yield of the glueballs, where we approximate m S Λ dyn . 11 The relic glueballs would dominate the energy density of the universe at the temperature T dom m S Y S if they are stable. To avoid large entropy production by the decay of the glueballs, we require so that T S 0 > T dom . We also require that S 0 decays before the era of the Big-Bang Nucleosynthesis. 12 Let us note here that S 0 decays more efficiently without requiring m ρ O(f a ) in the Higgs portal model discussed in the appendix A.
Finally, the axion decays into the Standard Model particles via the anomalous coupling in Eq. (8) (see Fig. 8). For m a O(1) GeV, the axion mainly decays into the QCD jets. For m a O(1) GeV, the axion decays into light hadrons through the mixing to the η and η mesons in the Standard Model [38]. It should be noted that the axion lighter than O(10 − 11 Excited glueballs have much smaller yields. 12 Even if T S0 < T dom , the present model provides a consistent dark matter model as long as this condition is satisfied. In this case, the resultant dark matter density is further reduced than the one in the following estimation.
The cross section of this process is expected to be about a geometrical one, where A = O(1). In fact, as discussed in Ref. [12], the heavy quarks inside the bound states are moving very slowly, v ∼ Λ dyn /M U when the baryons are colliding. Hence, the 13 Our assumption corresponds to the so-called the ∆-law, where the long-range potential is simply the sum of two-body potentials. See Refs. [34,40] for more on phenomenological potentials for baryons. quarks stay in overlap regions of the bound states for a long time, ∆t ∼ M U /Λ  dyn in the collisions. As a result, the quarks and anti-quarks are largely disturbed during the collision and they are well stirred. Eventually, the quarks and the anti-quarks are reconnected so that the baryons are broken into the mesons with O(1) probability in each collision (see Fig. 9).
Once the annihilation into the mesons happens, the mesons in the final state immediately decay into glueballs as discussed in the previous section.
With the above annihilation cross section, the Boltzmann equation of the total number density of the baryon, n B , is roughly given by, 14 By solving the Boltzmann equation, the number density of the baryons are reduced down leading to the relic abundance, 14 Here, σ B denotes the annihilation cross section of each baryonic bound state, which is roughly independent of the spins or any other internal degrees of freedom. Thus, if there are N B species of the baryonic bound states, the Boltzmann equation of the number density of each species, n = n B /N B , is given by, Here, the factor N c comes from the fact that the dark matter mass is M B N c × M U .
Therefore, the observed dark matter density, Ωh 2 0.1198 ± 0.0015 [41], can be explained by the dark matter mass in the PeV range.
In Fig. 10 Let us emphasize here that the number density of the quarks is conserved when the baryons annihilate into the mesons. The annihilation of the baryons just reconnects the quarks and anti-quarks inside the bound states. The actual reduction of the number of quarks happens when the meson decays. In this way, we can achieve a model of thermal relic dark matter with a mass lager than the unitarity limit although no interaction violates the unitarity limit microscopically.
The consistency with the unitarity limit can also be understood in the following way [10].
When the dark matter particle has a radius of R = O(Λ −1 dyn ), the highest partial wave that contributes to the collision is In this case, the annihilation cross section is bounded by the unitarity limit , This shows that the geometrical cross section in Eq. (26) is consistent with the unitarity limit.

IV. CONCLUSIONS AND DISCUSSIONS
In this paper, we discussed a model with thermal relic dark matter where the dark matter mass exceeds the so-called unitarity limit on the mass of point-like particle dark matter. In this model, the baryonic bound states are identified with dark matter, which possesses large radii when they are formed at the critical temperature around the the dynamical scale.
With the large radii, they annihilate into the mesons through a geometrical cross section.
The mesonic bound states decay into glueballs and axions which eventually decay into the Standard Model particles. As a result, we found that thermal relic dark matter with a mass in the PeV range is possible, which is beyond the usual unitarity limit. 15 15 It should be emphasized that the present paper does not require any entropy production to dilute the dark matter density. For a heavy thermal relic dark matter scenario with entropy production see e.g. [42].
One caveat is that we assumed the same quark potential in the mesons and the baryons in our discussion. If the binding energies of the baryons by the Coulomb potential are much larger than the mesons in Eq. (15), the size of the baryons at T c can be much smaller. In this case, the baryon annihilation cross section is expected to be smaller than the one in Eq. (26), and hence, the upper limit on the dark matter mass should be lower. If the binding energies of the baryons are smaller than the mesons in Eq. (15), on the other hand, the upper limit on the dark matter mass can be weaker. To derive precise upper limit on the dark matter mass, we need to solve the strong gauge dynamics with heavy quarks precisely, which is quite challenging with the current techniques.
In the model presented in this paper, we have the axion which couples to both the dark matter sector and the Standard Model sector. It is an interesting question whether the axion in the present model can play the role of the axion which solves the strong CP -problem by identifying U (1) A with the Peccei-Quinn symmetry [43][44][45][46]. Since the U (1) A symmetry is not only broken by the QCD but also by SU (N c ) which possesses its own θ-term, it is apparently difficult for the axion in this model to solve the strong CP -problem. However, if the SU (N c = 3) can be regarded as a counterpart of the QCD in a mirror copy of the Standard Model, 16 the θ terms in SU (N c = 3) and the QCD are aligned, so that the axion in the present model might solve the strong CP -problem [47][48][49][50][51][52]. Such a possibility will be discussed elsewhere.
Finally, let us comment on a possible phenomenological application of the present model.
In recent years, the IceCube experiment [21][22][23] has reported the excess in the observed flux of extraterrestrial neutrinos in the PeV range. Dark matter with a mass in the PeV range is considered to be one of the attractive explanation of the excess [53][54][55]. For example, the excess can be accounted for by dark matter with spin 3/2 and a mass 2.4 PeV which decays into neutrinos via for M * 5 × 10 34 GeV (corresponding lifetime of dark matter of O(10 28 ) s) [53]. Here, L and H represent the lepton and Higgs doublets in the Standard Model and ψ ν is dark matter with spin 3/2, respectively, its relic density cannot be explained by thermal relic density due to the unitarity limit. As we have discussed, however, thermal relic density can be consistent with the observed dark matter even for PeV dark matter. In fact, ψ ν can be identified with the baryons N c = 3. 17 Therefore, the IceCube results can be interpreted by the decay of PeV thermal relic dark matter in the present model. In the main text, we assumed that the dark matter sector is connected to the Standard Model dominantly through the axion. In this appendix, we consider an alternative model to connect the dark matter sector to the Standard Model sector via the Higgs portal. 18 For that purpose, we introduce two additional flavors of the fermions in addition to the U -quarks, and assume that they form the doublet representation of the SU ( 17 If the operator in Eq. (33) is provided by a Planck suppressed operator of the quarks, M * is expected to be much larger than M * 5 × 10 34 GeV. To provide appropriate M * , we need further extension of the model at the energy scale much larger than M U such as the emergence of conformal dynamics. 18 In this model, the U -quarks annihilates not into φ's but into gluons and/or higgs at the perturbative freeze-out, which does not affect the thermal history after the confinement.
with a spin 3/2. Due to the spin-spin interaction, we expect that ∆ is heavier than N by Furthermore, the neutral baryon U DD is lighter due to the U (1) Y interaction, by, Therefore, in this case, the lightest baryon is expected to be U DD in N , which is neutral under U (1) Y and can be identified with dark matter. 19 To make unwanted charged particles in the dark matter sector decay, we introduce a light complex scalar field s which has a hypercharge −1 and the following coupling, L = y s UD + h.c.
Though this interaction, the mesons decay into s's (and glueballs) and the heavier baryons decay into the lightest baryon by emitting s's. Finally, s decays into the Standard Model sector via, for example, where M * denotes a dimensionful parameter. 20 One might be interested in a model where (U, D) and (Ū ,D) form the doublets of the 19 Due to the radiative corrections of U (1) Y gauge interactions, the mass of the D quark is expected to be smaller than that of the U quark. Here, we assume that the masses of U 's and D's are finely tuned so that the mass differences between the baryons are mainly given by Eqs. (B4) and (B5). 20 Since s is charged under U (1) Y , it should be heavy enough so that the constraints from the collider experiments are avoided.