GUTzilla Dark Matter

Motivated by gauge coupling unification and dark matter, we present an extension to the Standard Model where both are achieved by adding an extra new matter multiplet. Such considerations lead to a Grand Unified Theory with very heavy WIMPzilla dark matter, which has mass greater than ~10^7 GeV and must be produced before reheating ends. Naturally, we refer to this scenario as GUTzilla dark matter. Here we present a minimal GUTzilla model, adding a vector-like quark multiplet to the Standard Model. Proton decay constraints require the new multiplet to be both color and electroweak charged, which prompts us to include a new confining SU(3) gauge group that binds the multiplet into a neutral composite dark matter candidate. Current direct detection constraints are evaded due to the large dark matter mass; meanwhile, next-generation direct detection and proton decay experiments will probe much of the parameter space. The relic abundance is strongly dependent on the dynamics of the hidden confining sector, and we show that dark matter production during the epoch of reheating can give the right abundance.


I. INTRODUCTION
Grand unified theories (GUTs) [1,2] are one of the most attractive and well-studied scenarios for physics beyond the Standard Model (SM). With just the particle content of the SM, the three gauge couplings run tantalizingly close to one another at around 10 15 GeV. However, they do not meet at a single scale. The possibility of a GUT thus motivates additional new matter below the GUT scale, which can modify the running and allow unification. In principle, such new matter may be present anywhere above the weak scale up to the GUT scale, and there are limitless possibilities.
On the other hand, the existence of dark matter (DM) provides one of the strongest signs of physics beyond the SM. The existing searches for DM have dominantly focused on weak-scale thermal relics within the Weakly Interacting Massive Particle (WIMP) paradigm [3]; however, the lack of definitive signals from (in)direct detection experiments [4,5] and at the Large Hadron Collider (LHC) [6,7] have placed increasingly stringent constraints on WIMP models. Therefore, it is prudent to re-examine our theory assumptions and explore alternative DM beyond the WIMP, including DM at much higher mass scales.
In this paper, we propose an extension to the SM that gives both gauge coupling unification and a very heavy DM candidate, with mass well above the weak scale. We extend the SM with an additional matter multiplet χ, which is part of a larger split GUT multiplet. In order for χ to give successful gauge coupling unification, the multiplet must be both electroweak and color charged. Due to the color charge, we are led to consider χ charged under an additional confining hidden gauge group. The DM is then a composite state with electroweak interactions, which can evade direct detection bounds for masses above 10 7 GeV.
We refer to such GUT-motivated heavy DM as GUTzilla DM, by analogy with the WIMPzilla DM scenario [8][9][10][11][12]. In general, very heavy DM cannot be produced thermally during a radiation-dominated era; if the DM were in thermal equilibrium, then a large annihilation rate would be required to avoid overclosing the universe, which runs into unitarity bounds for DM masses above the 100 TeV scale [13]. Instead, in the WIMPzilla scenario, the relic abundance is set before the end of reheating. Then the relic abundance is naturally suppressed if the reheating temperature is smaller than the DM mass, and it is possible for the DM mass to span many orders of magnitude. Specifically, we consider DM production in inflaton direct decay and production from the SM thermal bath during the reheating epoch.
Other models that accomodate both unification and DM with extensions to the SM have been studied before [14][15][16][17], including SO(10) unification [18,19]. These models typically require multiple new particles and hierarchies of scales. Here, we will consider the scenario with only one new hierarchy associated with the χ multiplet. The χ multiplet will invariably be part of a split GUT multiplet, and we will simply assume that the splitting is accomplished by a minimal amount of fine-tuning. We will comment on a possible connection between the fine-tuning and the anthropic principle later.
Our paper is structured as follows: In Sec. II, we outline the requirements on the new matter multiplet χ such that gauge coupling unification is achieved; we show that in order to satisfy proton decay constraints, χ needs to have both electroweak and color charge. A viable DM candidate would then need to be a composite state composed of the χ, which we assume is the result of a new confining gauge interaction. In Sec. III, we construct a minimal model of GUTzilla DM, where χ is a fundamental of a confining SU (3) H . For simplicity, our discussion will only focus on the scenario where the confinement scale, Λ H , is smaller than m χ , such that non-perturbative physics is less important. Our model is also viable for larger Λ H , and we briefly discuss this possibility. In Sec. IV, we present the predictions for direct detection and proton decay signals, finding that current constraints require DM masses of at least 10 7 GeV. We also calculate the hidden-sector contributions to the Higgs potential, finding an improved stability of the electroweak vacuum. We turn to the cosmology of such heavy DM in Sec. V, where we discuss DM production by inflaton decay in Sec. V A and production from the SM thermal bath in Sec. V B. We summarize our findings and comment on future directions in Sec. VI. An additional mechanism for DM production during reheating is discussed in Appendix A.

II. GAUGE COUPLING UNIFICATION
The SM has gauge group SU (3) C × SU (2) L × U (1) Y , where the coupling strength for each of these gauge interactions receives scale-dependent quantum corrections. At one loop, they are given by 1 where b SM a is the running coefficient in the SM, Fig. 1 illustrates the SM running coupling at one-loop. The three couplings approach one another at around µ ∼ 10 15 GeV but do not unify. Gauge coupling unification then requires new matter lighter than ∼ 10 15 GeV or a huge threshold correction around the unification scale.
The simplest way to achieve unification is by adding a new matter field χ with mass m χ 10 15 GeV, as illustrated in the right panel of Fig. 1. The addition of χ modifies the running couplings in Eq. 1 as 1 The hypercharge coupling strength is defined as α 1 = 5g 2 Y /12π, where the extra factor of 5 /3 comes from embedding U (1) Y ⊂ SU (5). The coefficient b χ a can be written as b χ a = −N χ s χ c a , where N χ is the overall multiplicity of the χ multiplet, c a is the sum of the Dynkin index, and s χ is the spin factor defined as Unification requires that all three couplings meet 2 at some scale Λ GUT , i.e. α a (Λ GUT ) = α GUT . Given that the gauge coupling unification is only dependent on the combination b χ a log mχ Λ GUT at one loop, Λ GUT and α(Λ GUT ) are therefore invariant under the transformation, where n is an arbitrary constant. Using the transformation in Eq. 5, m χ can be raised arbitrarily close to Λ GUT .
More generally, since χ is part of a unified multiplet, the conditions for unification will be modified depending on the mass scale of the multiplet. Let χ ⊕ χ be a complete GUT multiplet and assume that the masses of χ are all of the same order, m χ . The requirement of coupling unification using Eq. 3 implicitly assumes that m χ = Λ GUT , which can be relaxed. For fermionic χ, m χ is naturally below Λ GUT . Above the scale m χ , the extra running from χ ⊕ χ comes in a complete multiplet and does not affect unification. Gauge coupling unification then depends only on b χ a log mχ m χ at one loop, and thus will be invariant under the 2 We do not consider the strong-coupling unification scenario [20][21][22] nor the unification involving higherdimensional operators [23,24].
following set of transformations: The transformations in Eq. 6 preserve both α(Λ GUT ) and Λ GUT . However, Eq. 7 keeps Λ GUT the same, but modifies α(Λ GUT ) due to the extra running between (m χ , Λ GUT ). So while the full running is a function of (m χ , m χ , N χ s χ ), the transformations above demonstrate that Λ GUT is solely dependent on the SM representation of χ at one loop. A preliminary analysis of allowed Λ GUT can thus immediately place restrictions on allowed representations for χ. Perturbatively, the presence of additional matter generally causes the gauge coupling to run larger. If the χ multiplet is SU (2) L neutral, then 2π/α 1 and 2π/α 3 are pushed to smaller values; comparing with Fig. 1, we see that no solution with unification is possible in this case. Next, consider a χ multiplet which is SU (3) C neutral: while unification is now possible, the maximum possible Λ GUT is where α SM 1 and α SM 3 intersect, or Λ GUT < 3.4 × 10 14 GeV. However, Λ GUT directly controls the proton decay rate, and current bounds require a GUT scale of at least O(10 15 ) GeV [25,26]. We conclude that χ must be charged with respect to both SU (3) C and SU (2) L . 3 How does a color-charged χ give rise to a DM candidate? Without additional structure, the color-charged χ particles will form bound states with light quarks, resulting in strong interactions with ordinary matter. Such strongly interacting DM has already been ruled out by earth heating and direct detection bounds [4,27,28]. On the other hand, the DM-nucleon interaction will be suppressed if the Bohr radius of the bound state is much smaller than 1/Λ QCD , which can be achieved by adding a new confining gauge group G H with a large confinement scale Λ H . The DM is then a neutral composite bound state of χ, though it generally has nonzero hypercharge. Stringent direct detection bounds due to Z exchange then leads us to consider a very heavy, non-thermal DM candidate. We turn to specifics of this scenario in the following section.

III. MINIMAL GUTZILLA
In the minimal GUTzilla DM model, we add to the SM an extra Dirac fermion multiplet χ ⊕ χ , where m χ < m χ . In Sec. II, we have shown that χ needs to be charged under both SU (3) C and SU (2) L . The smallest such representation for χ is (3, 2)1 /6 , which is a subset of the 10 and 15 representations of an SU (5) GUT. 4 (While our considerations do not depend on the unification group, we will use the language of SU (5) for simplicity.) In order to form color-neutral DM, we introduce a hidden sector gauge group SU (N H ) (with N H = 3), which confines at Λ H . 5 The χ ⊕ χ multiplet transforms as the fundamental representation of SU (3) H . Then the GUTzilla DM is a stable baryonic state composed of three χ fermions. Depending on Λ H , the composite sector also contains new meson or glueball states, which decay quickly into SM particles. Depending on the hierarchy between Λ H , m χ and m χ , our previous analysis of gauge coupling unification may be modified. Our model has three distinctive physical regimes: • Λ H < m χ < m χ : The gauge running computation is simplest for this hierarchy, with heavy χ. The SM gauge couplings receive new contributions at the scale m χ , where the hidden sector coupling is perturbative, and it is straightforward to determine the running at one-loop. The hidden baryons are composites of the χ fermions, such that m DM ≈ 3m χ .
• m χ < Λ H < m χ : In this case, hidden sector pions π H are present and the SM gauge couplings are modified at the scale m π H . The running between m π H and Λ H can be calculated in chiral perturbation theory. Non-perturbative physics comes in around the confinement scale and will introduce extra threshold corrections. For scales larger than Λ H , the perturbative one-loop analysis applies again. One can estimate the correction to the running in the chiral regime. Since the SM gauge group explicitly breaks the chiral flavor symmetry, the pion masses are only smaller than Λ H by a loop factor. Then the change to 2π/α SM is at most of order | log(α SM /4π)| 5. Given the small running coefficient due to scalars, such a contribution is subdominant to potential threshold corrections near Λ GUT . We will then treat this scenario in the same way as the heavy-χ case, keeping in mind that the running calculation applies with the substitution m χ → Λ H and a large uncertainty exists from extra running due to pions and other non-perturbative composite states.
• m χ < m χ < Λ H : The changes to the SM gauge coupling running mainly arise from the light pions in the hidden sector. Again, the pion masses have large contributions due to SM gauge interactions and are one-loop suppressed compared to Λ H . The resulting pion spectrum has small mass splittings, and thus the modification of the SM running is too small to achieve unification.
For concreteness, we discuss below the case where Λ H < m χ < m χ , such that we may follow the perturbative one-loop analysis. As noted above, for m χ < Λ H < m χ , the gauge coupling unification is very similar as long as we make the identification m χ Λ H in the running calculation. For a Dirac fermion χ in the representation (3, 2)1 /6 , the contribution to the running in Eq. 3 is 5 There are other possibilities for the hidden sector gauge group, such as SO(2N ) H . We briefly comment on this in Sec. VI.
Assuming coupling unification, the GUT scale is given by Λ GUT = 3 × 10 15 GeV, and the mass hierarchy between m χ , m χ is given by where N χ s χ = 3 × 4/3 for our model. Such a small m χ /m χ can be achieved by tuning a Yukawa coupling of the χ ⊕ χ multiplet with a GUT-breaking Higgs field. Unification can also be achieved for scalar χ ⊕ χ , which we do not discuss here. The confinement scale of the hidden sector Λ H is in general a free parameter. For example, for χ ⊕ χ transforming as a 10 multiplet of SU (5), the renormalization group equation of the gauge coupling of SU (3) H is given by The gauge group SU (3) H will remain asymptotically free and Λ H can range all the way from Λ H 1 GeV to 10 14 GeV for moderate coupling at GUT scale, α H (Λ GUT ) ∈ (0.01, 1).
For the regime we are interested in, the inverse radius of the composite particle is much larger than the electroweak scale. Then the doublet (DM, DM + ) is essentially an elementary particle at low energies, and the DM-nucleon scattering rate is dominated by Z-exchange. Electroweak symmetry breaking effects will induce a mass splitting for the doublet, which is independent of m DM [30]: The charged DM + particle can decay through an off-shell W + , which can lead to a soft pion or leptons. The two-body decay DM + → DM + π + dominates, with a rate given by such that the DM + easily decays away before Big-Bang Nucleosynthesis (BBN).
Meanwhile, the stability of the neutral DM state can be guaranteed by symmetries. One could simply impose a Z 2 charge (−1) χ , or a continuous U (1) χ ⊃ (−1) χ . These symmetries can also be obtained within SO(10) unification. For example, one can embed the χ within a 45 or 54 multiplet of SO(10), so that χ has a U (1) B−L charge of 2/3. When SO(10) is broken into SU (5) × U (1) X by a 126 Higgs, a discrete (−1) 3(B−L) remains unbroken. Then (−1) χ can be identified with (−1) 3(B−L)+F , where F is fermion number. Geometrically, such a parity can be thought of as the spinor parity in SO (10). Analogous to the Lorentz group, representations of the Lie algebra so(10) are actually representations of the universal cover Spin (10), where Spin(10)/Z 2 = SO(10), and all spinor representations (SM fermions) are charged under the extra Z 2 . Then (−1) χ can be identified with In addition to the GUTzilla DM, there are additional composite states arising in the hidden sector. The physics of these states depend on Λ H and m χ . While our main focus is on the heavy-χ scenario (Λ H < m χ ), we will also discuss the alternative QCD-like case (Λ H > m χ ) for completeness. Both scenarios provide a stable DM candidate and similar low-energy phenomenology.

B. Heavy-χ scenario
When the hidden confinement scale Λ H is much smaller than m χ , the lightest hidden sector states are glueballs with various spin and quantum numbers [31]. There are additional heavier meson states which decay rapidly into glueballs and SM gauge particles. The lightest glueball is a scalar and can decay back into the SM through dimension-8 operators obtained from integrating out the χ ⊕ χ [32]. These operators can be written schematically as where the c i are O(1) coefficients. The higher spin terms include non-trivial tensor contractions between the SM field strengths and the higher spin glueball fields. The dimension-8 operators induce decay of the scalar glueballs into SM gauge bosons, with a rate of order As long as Λ H is sufficiently large, the glueball will decay well before BBN (∼ 1 sec): If CP is conserved, there are additional higher spin states that can only decay radiatively, which will lead to a stronger bound for Λ H . Then Eq. 15 will serve as a conservative bound for the hidden sector confinement scale.

C. QCD-like scenario
When the confining scale of the hidden sector is larger than m χ , the hidden sector undergoes chiral symmetry breaking. The light degrees of freedom are pseudo-Nambu-Goldstone bosons, or pions π H . From the perspective of the hidden sector, there is an approximate SU (3×2) L ×SU (3×2) R global symmetry explicitly broken by SM gauge interactions. Below the confinement scale, this flavor symmetry is spontaneously broken to the diagonal group There are a total of 35 pion fields, and they reside in SM representations given by (8,3) The DM again is a baryon doublet, but its mass is dominated by the confinement scale Λ H instead of the masses of its constituents, where N H = 3.
The pions can decay through dimension-5 operators in chiral perturbation theory, which easily satisfies BBN constraints for the DM masses considered.

IV. PHENOMENOLOGY
In this section, we consider the main phenomenological implications of GUTzilla DM, namely DM direct detection and proton decay. Strong constraints on direct detection experiments require a large DM mass, m DM 10 8 GeV. Additionally, we comment on the modification to the Higgs potential and vacuum stability, finding that the inclusion of the hidden sector improves stability.

A. Direct Detection
Given that the DM has a non-vanishing U (1) Y charge, it interacts with a nucleus via tree-level exchange of a Z boson. For a given nucleus N , the average per-nucleon scattering cross section is given by where G F is the Fermi constant, µ n is the reduced mass of the nucleon and DM, Y is the hypercharge of DM, and A N and Z N are the atomic number and charge of the nucleus, respectively. Such an interaction is highly constrained by direct detection experiments; in the high-mass limit, the tightest bounds come from the LUX experiment [4]: This constraint on the DM mass translates to constraints on the hidden sector. In the case Λ H m χ , the DM mass bound leads to a bound on m χ m DM /3. For Λ H m χ , the DM mass bound leads to a bound on the hidden sector confinement scale Λ H m DM /N H . Together with gauge coupling unification, m DM can roughly be in the range 10 8 to 10 12 GeV.
In Fig. 2 we show the sensitivity of direct detection and proton decay experiments to our model. We show the constraint from LUX (2015) [4] by the vertical solid line, assuming the m DM = 3 × m χ . We also show the projected sensitivity of the LZ experiment [33], and a direct detection experiment whose sensitivity is limited by the neutrino background [34]. It can be seen that a large portion of the parameter space can be tested by future direct detection experiments. With multiple target nuclei, it is also possible to test whether a DM candidate interacts via Z exchange [11], which would point towards very heavy DM as in Eq. 20.

B. Proton Decay
In a generic GUT, there are new heavy particles at the unification scale that can mediate proton decays (for a review, see Ref. [35]). In an SU (5) GUT, a proton can decay into a meson and a lepton via the exchange of gauge bosons charged under both SU (3) C ×SU (2) L ⊂ SU (5). The most stringent proton decay constraint comes from the p → π 0 + e + channel. In the effective field theory below the GUT scale, such a decay arises from the following dimension-6 operator, where M XY is the mass of GUT gauge bosons, g GUT is the gauge coupling constant at the GUT scale, and A R,L is the Wilson coefficient from renormalization group running from the GUT scale down to the hadron scale. Here we have neglected the effects of the quark mixing angles since this depends on how quark masses are unified. We have also ignored potential contributions from the Higgs sector, which are Yukawa-coupling suppressed. Then the decay rate of the proton is given by where W 0 is the quantity encoding the form factor of a pion and a proton. Lattice calculations show that |W 0 | = 0.103 GeV 2 at the renormalization scale of 2 GeV, with uncertainty of 40% [36].
In Fig. 2, we show the prediction for the decay rate of the proton in our model as a function of m χ . The lower green line shows the prediction for M XY = Λ GUT = 3 × 10 15 GeV. The upper green line shows the prediction for M XY = 7×10 15 GeV, which can be achieved by a moderate threshold correction around the GUT scale, ∆(2π/α) = 5, as we discuss below. In fixing α GUT , we assume that χ is embedded into a 10 of SU (5). If χ is embedded into a larger representation, α GUT is larger and the proton decay rate becomes larger. The light blue shaded region is excluded by Super-Kamiokande, Γ −1 (p → π 0 e + ) > 1.7 × 10 34 years (90%CL) [26]. We show the expected sensitivity of Hyper-Kamiokande, Γ −1 (p → π 0 e + ) > 1.3 × 10 35 years (90%CL) [37]. It can be seen that the entire parameter space can be covered by Hyper-Kamiokande. The calculations for M XY , A L,R and the treatment of threshold corrections are described below.
a. Estimation of M XY The masses of the X/Y gauge bosons are typically of order Λ GUT . If threshold contributions to the running are present, then M XY can be raised and the proton lifetime can be increased. Generally, an accurate estimate for M XY requires taking into account any additional split multiplets around Λ GUT and/or higher order corrections to the running couplings. In order to account for these model-dependent corrections, we simply relax the coupling unification requirement by varying the mass ratio m χ /m χ and allowing the couplings to differ by some amount. Thus, we can determine M XY by demanding that for all the SM couplings, and where ∆ max parameterizes the deviation from unification. With only the SM particle content, the minimal value of ∆ max is 25 around M XY ∼ 10 14 GeV, which is ruled out. With the additional χ multiplet, the case ∆ max = 0 corresponds to the scenario of no threshold corrections with M XY = Λ GUT = 3 × 10 15 GeV. Allowing ∆ max = 5, M XY can be in the range 1.5 × 10 15 to 7.0 × 10 15 GeV. The green curves in Fig. 2 show the proton decay lifetime given two values of M XY , corresponding to ∆ max = 0 and the upper bound on M XY for ∆ max = 5.
b. Estimation of A R,L The dimension-6 operators in Eq. 21 obtain anomalous dimensions from gauge interactions. Under renormalization group evolution, the Wilson coefficients receive significant multiplicative corrections. The coefficient A R,L at different scales is then related by [38,39]: Taking A L,R (Λ GUT ) = 1, and M Λ GUT M XY = 10 15 GeV, one obtains A SM R (2 GeV) 3.0 and A SM L (2 GeV) 3.4 in the SM. In our GUTzilla model, the introduction of χ ⊕ χ only modifies A L,R at the percent level. Since possible threshold corrections dominate the uncertainty in the proton lifetime, we simply take A L,R to be the SM values in our calculation.

C. Vacuum Stability
In the SM, the Higgs quartic coupling receives a large negative contribution from the top Yukawa coupling, which can lead to a meta-stable or unstable electroweak vacuum [40][41][42]. Given the current Higgs mass and top mass measurements, an NNLO calculation for the Higgs potential has firmly excluded SM vacuum stability at the 2σ level (see Ref [43][44][45] and references therein). The SM Higgs quartic becomes negative at around 10 11 GeV; the presence of the additional χ multiplet with m χ 10 11 GeV could increase the gauge coupling and improve the stability of the Higgs potential. For our minimal GUTzilla model, a small m χ 6 × 10 7 GeV is needed to stabilize the Higgs potential within 1σ. The Higgs quartic running is illustrated in Fig. 3, which shows the quartic coupling including the leading-order effect of χ, χ . We use the central value for the Higgs mass, and show the effect of varying the top-quark pole mass within ±1σ.

V. GUTZILLA COSMOLOGY
A massive DM candidate in thermal equilibrium during a radiation-dominated era is easily overproduced; unitarity limits on the DM annihilation cross section require m DM 300 TeV [13]. Instead, processes before the end of reheating can set the abundance of GUTzilla DM and thus get around this bound, as long as the reheating temperature T RH is less than the DM mass. In this section we describe the various possibilities.
Superheavy DM may be produced gravitationally during the transition from an inflationary phase to a matter-dominated era [10,46,47]. This mechanism is sufficient for producing the correct relic abundance of DM if both T RH and the Hubble scale at the end of inflation, If m DM < m φ , production during reheating is also possible [8][9][10][11][12][49][50][51]. There are three possible mechanisms in play: inflaton decay, thermal production, and inelastic scattering between inflaton decay products and the SM plasma. DM production from inflaton decay, shown in Fig. 4, will be important as long as it is kinematically accessible. For heavy GUTzilla DM with m DM 10 8 GeV, overproduction of DM will then place constraints on the reheating temperature [12], which we will discuss in Sec. V A. However, these constraints can be evaded if the inflaton dominantly decays to SM singlets which are lighter than the DM, thus shutting off the previous mechanism. Then other possibilities such as thermal production and/or inelastic scattering of inflaton decay products will become important. We will discuss the thermal production channel in Sec. V B. The inelastic scattering case is highly model-dependent, and we give a simplified treatment in Appendix A.
In each of these scenarios, GUTzilla DM production is further complicated by the hidden sector dynamics. If the temperature of the hidden sector at the time of DM production is smaller than Λ H , GUTzilla DM will be directly produced. Otherwise, the constituents χ will first be produced and the DM bound states are formed only after the hidden sector confining phase transition. As before, we will primarily focus on the heavy-χ scenario, m DM Λ H . Note that in the QCD-like scenario, non-perturbative processes produce a substantial amount of DM, leading to a strong constraint on T RH unless m DM > m φ . Fig. 5 illustrates the parameter space to produce GUTzilla DM, depending on the production mechanism, T RH , and m DM . The derivation of these bounds can be found in the remainder of this section. For direct inflaton decay, a low T RH is required in order to avoid over-producing DM; this constraint is shown as the dotted-black line. If the inflaton only couples to SM singlets, and the singlet has suppressed coupling to SM and hidden sector states, a looser constraint from thermal production applies and is shown in the red region. In GeV. The dashed-black line shows the upper limit on T RH coming from direct decays of the inflaton (Sec. V A). When direct decays of the inflaton to DM are turned off, the thermal production of DM during reheating can give the dominant relic abundance (Sec. V B) -in this case, the red-shaded region indicates where Ω DM h 2 is greater than the observed value. For the constraint from direct decays, Eq. 29 is used for an estimate of the average DM multiplicity N DM in inflaton decay, α is fixed at 0.05, and we include an additional factor of (1 − 4m 2 DM /m 2 φ ) 1/2 to take into account phase-space suppressions.
this case, the correct DM relic abundance can readily be obtained for high T RH . In addition, the constraints can be evaded when DM becomes heavy enough such that it is kinematically inaccessible. In the direct decay case, this happens when 2m DM > m φ , shown as the dashed black line.

A. Inflaton Decay
When the inflaton directly decays to SM-charged particles, production of DM can proceed by gauge boson emissions and subsequent splitting into χ particles. 6 An illustrative diagram is shown in Fig. 4. Due to the large entropy production during the reheating period, DM production from inflaton decay is most prominent at the end of reheating, around T = T RH . The DM relic density from inflaton decay can be estimated as where T RH is the reheating temperature, m φ is the inflaton mass, and N DM is the average number of composite baryonic DM per inflaton decay. Generally, N DM depends on the inflaton coupling. We will take the most conservative approach and assume there is no direct coupling of the inflaton with χ. However, as long as the inflaton primarily decays into SM charged particles, the decay products can undergo showering and radiate hidden sector particles, which can eventually hadronize into DM. At high energies, these showering processes are perturbative and can be calculated systematically [52].
Consider the average number of χ particles produced in the shower. Aχχ splitting is necessarily preceded by a gauge boson emission. Then at leading order, the average number of χ from an inflaton decay is where N a gauge is a splitting kernel given by and a denotes SM and hidden sector gauge bosons. N a gauge is the average number of gauge bosons at leading order, and α a and C a F are the corresponding gauge coupling and Casimir for gauge boson emission. Note the form of N a gauge is valid only for a non-Abelian gauge group, which is assumed to dominate the shower. Performing the integration in Eq. 26 gives For m φ m χ , resummation of the large logarithm results in an exponential enhancement [52]. We find that as long as m 2 φ /4m 2 χ 10 10 , the perturbative estimate in Eq. 26 is sufficient.
For GUTzilla models, the DM is composed of three χ, which will require three separate gauge boson splittings toχχ in addition to a suppression factor to form a baryon, which we take to be ∼ 1/N 2 H . The average DM multiplicity is then estimated to be For a typical SM interaction, α ∼ 0.05, and for m φ /m χ ∼ 10 5 , we have N DM ∼ 10 −4 .
So far we have only included perturbative contributions from showering and Eq. 29 ignores contributions from non-perturbative processes, which is valid in the heavy-χ limit. In the opposite QCD-like regime where m χ /Λ H is small, non-perturbation fragmentation and hadronization can also produce baryons, leading to a large DM multiplicity.
In the Lund string model, baryon fragmentation can be thought of as breaking of the gluon string by a diquark/anti-diquark pair, with a fragmentation function of the form ∼ exp(−4m 2 χ /Λ 2 H ) [53,54]. In the light quark regime, diquark fragmentation is not significantly suppressed and the DM multiplicity will be of order the hidden gluon multiplicity N DM ∼ 2α log 2 (m 2 φ /4m 2 χ )/27π ∼ 0.01. In the heavy-χ limit m χ Λ H , diquark fragmentation is exponentially suppressed and does not contribute to baryon production.
Given our very conservative estimate of N DM , in order to avoid overproducing the DM relic density ρ DM /s 10 −9 GeV, we find the reheating temperature is constrained to be The black dashed curve in Fig. 5 shows the regions of parameter space excluded from overproduction of the DM, where Eq. 29 is used for an estimate of the average DM multiplicity from inflaton decay. We see that even a suppressed N DM can lead to tight constraints on the reheating temperature. The constraint can be relaxed, however, if no direct coupling between the inflaton and SM charged particles exists; then sub-dominant processes become important, as we discuss below and in Appendix A.

B. Thermal Production
If the inflaton only couples to SM singlets S (e.g., a right-handed neutrino), the decay into DM will have to proceed through the coupling between S and other SM particles, which can be highly suppressed. Subsequent decay from S → DM can be forbidden as long as m S < m DM . Then production of DM from the SM thermal bath is relevant, and the right relic density can be achieved thanks to dilution from entropy production by the inflaton.
During reheating, the inflaton gradually transfers energy to the SM plasma. The SM bath will heat up to a maximum temperature T max > T RH , while the energy density of the universe is still dominated by that of the inflaton. The energy density then becomes dominated by the relativistic SM bath at T RH . As long as m χ < T max , then the χs can be pair produced from the SM thermal bath via gauge interactions. The comoving number density of χ freezes out when T ∼ m χ , and they can later be bound up into DM baryons when the temperature drops below Λ H .
For light DM, the χ particles are in thermal equilibrium during the inflaton-dominated era. Following Ref. [55], the DM relic density is given by where x f = m χ /T f , and T f is the freeze-out temperature. The annihilation cross section of χ is σv 4πα 2 /m 2 χ , and we have x f 10 + log T RH 80 TeV To avoid over-production of DM, the reheating temperature is bounded by The red shaded region in Fig. 5 shows the reheat temperatures excluded for thermal production. When the DM mass is larger than ∼ 5 × 10 10 GeV, χ is not in chemical equilibrium. Then out-of-equilibrium production and inelastic scattering (Appendix A) processes may contribute to the DM abundance; since these are much more model-dependent, we have not shown these constraints. Lastly, here we have assumed that T max is always larger than m DM and that kinetic equilibrium is established. A simple estimate gives T max ∼ (m φ M Pl ) 1/4 T 1/2 RH [55], which is well above the DM mass for the parameter space shown here. However, a more detailed recent analysis shows that thermalization can be slower, with a lower T max ∼ α [56]; depending on the specifics of this thermalization process, we expect the excluded region will be modified somewhat.

VI. CONCLUSION
In this paper, we investigated a new class of models linking gauge coupling unification and DM through the introduction of a single multiplet. In order to achieve unification, the new multiplet must be a split GUT multiplet and the lighter component must be SU (3) C and SU (2) L charged. This prompted us to include a hidden confining sector to screen these interactions and leads to a composite baryonic DM. These DM can be very heavy and thus evade direct detection constraints, and provides a new motivation for considering the heavy DM WIMPzilla scenario. We refer to this as the GUTzilla DM scenario.
We presented a minimal implementation of GUTzilla DM by adding a split Dirac fermion χ ⊕ χ multiplet, where χ transforms as (3, 2)1 /6 under the SM gauge group and as a 3 under an SU (3) H hidden gauge group. The DM is then a baryon state made of three χ. While our considerations do not explicitly depend on the hierarchy between Λ H and m χ , we focused on the heavy-χ case (Λ H < m χ ) for simplicity.
Phenomenologically, the most prominent signatures of GUTzilla DM are direct detection and proton decay. The current direct detection limit points to a GUTzilla DM with masses at least of order 10 8 GeV, which will be readily tested at future LZ and Hyper-Kamiokande experiments. We also show that the addition of GUTzilla DM can improve the stability of the Higgs potential to within 1σ as long as the DM is not too heavy.
The relic abundance of the DM is set before the end of reheating. For DM mass larger than the Hubble scale at the end of the inflation, the abundance is saturated by gravitational production if both the Hubble scale and the reheating temperature are large enough. For DM mass smaller than the inflaton mass, production of DM during reheating is possible and can put a tight constraint on the reheating temperature. If the inflaton directly decays into SM charged particles, DM is easily overproduced unless the reheating temperature is very low. In the heavy-χ scenario we are considering, suppression of baryon production helps to alleviate these constraints. On the other hand, we show that a large reheating temperature is still possible if the inflaton decays to SM singlets, assuming these singlets do not have large direct coupling to the hidden sector. In this case, the DM abundance is saturated by thermal production during reheating.
There are many variations on the minimal GUTzilla DM that could be considered. In this paper, we introduced an SU (3) confining gauge group to obtain an electromagnetic and color neutral baryon from χ's. One possibility is to introduce an SO(2N ) gauge group, where χ is in a fundamental representation. The lightest baryon is expected to be composed of χ N and χ †N , and would be neutral under the SM gauge group. The DM can then be much lighter than the GUTzilla mass range. However, a possible problem of this model is an existence of a meson composed of two χ, which may be stable due to accidental χ number conservation and hence cause cosmological problems. This problem could be avoided by introducing a higher dimensional operator breaking the accidental χ number conservation. We defer further discussion of this model to future work.
Finally, one may ask whether the addition of split multiplets for gauge coupling unification introduces additional fine-tuning. This could be addressed by the anthropic principle [57,58], by attributing the fine-tuning m χ /m χ 1 to the necessity of obtaining enough DM for structure formation [59,60]. For a fixed reheating temperature and an inflaton mass, m χ should be small enough to obtain the DM density. The mass splitting is explained if m χ is biased toward a high energy scale. Note that any further mass splitting within χ is disfavored, as it requires unnecessary fine-tuning as far as the DM abundance is concerned. From the landscape point of view [61][62][63], this explanation of the splitting requires that there is no habitable vacuum with a less fine-tuned parameter set. To put it the other way around, if GUTzilla DM is present in our universe, we may infer restrictions on the landscape of parameters related with the abundance of DM, e.g. the inflaton mass, the reheating temperature, and the decay constant of a QCD axion [64][65][66]. is the decay rate of S including the Lorentz boost factor. An extra factor m φ /E ψ is included as a rough estimate of the multiplicity factor. In the limit that Γ S H, Γ φ , the number density of S will reach an equilibrium density with Γ S n S ∼ Γ φ n φ . Effectively, one can ignore the intermediate S state and treat the inflaton as a source for production of ψ. A steady state solution will be reached with n ψ m φ Γ φ n φ /(E ψ Γ split ) (for Γ split H). Taking coherence effects [67,68] into account, the splitting rate roughly follows Γ split ∼ α 2 T 3 /E ψ [56,69].
To obtain the DM density at a temperature T , we use Eq. A2 and substitute in n SM g * T 3 for the thermal bath. For the cross section to produce DM, we take σv ∼ α 2 N DM (ŝ) /(E ψ T ), where N DM (ŝ) denotes the average DM multiplicity per inelastic scattering event atŝ = E ψ T (see Eq. 29). Note that this is valid for T < Λ H , such that the baryons are directly produced in the collision. Assuming that DM is produced primarily at a single temperature T , we then rescale the number density at T to that at T RH , below which the comoving DM density freezes out.
Keeping only the leading power-law dependence, the resulting DM abundance is given by In general, the full DM production must be integrated over the allowed range of T ∈ [T RH , T max ] and the allowed energy range of E ψ ∈ [m 2 DM /T, m φ ]. Note that due to the log-enhancement in N DM (ŝ) favoring larger E ψ , the production rate typically peaks at an intermediate energy. Maximizing over T and E ψ in Eq. A4 to estimate the DM relic abundance, we have found that the inelastic scattering gives similar or somewhat lower relic abundance compared to thermal production. While the thermal production mechanism suffers from a larger entropy dilution, the inelastic scattering has a large suppression from N DM 10 −10 .