Gravitational Vector Dark Matter

A new dark sector consisting of a pure non-abelian gauge theory has no renormalizable interaction with SM particles, and can thereby realise gravitational Dark Matter (DM). Gauge interactions confine at a scale $\Lambda_{\rm DM}$ giving bound states with typical lifetimes $\tau \sim M_{\rm Pl}^4/\Lambda^5_{\rm DM}$ that can be DM candidates if $\Lambda_{\rm DM} $ is below 100 TeV. Furthermore, accidental symmetries of group-theoretical nature produce special gravitationally stable bound states. In the presence of generic Planck-suppressed operators such states become long-lived: SU$(N)$ gauge theories contain bound states with $\tau \sim M_{\rm Pl}^8/\Lambda^9_{\rm DM}$; even longer lifetimes $\tau= (M_{\rm Pl}/\Lambda_{\rm DM})^{2N-4}/\Lambda_{\rm DM}$ arise from SO$(N)$ theories with $N \ge 8$, and possibly from $F_4$ or $E_8$. We compute their relic abundance generated by gravitational freeze-in and by inflationary fluctuations, finding that they can be viable DM candidates for $\Lambda_{\rm DM} \gtrsim 10^{10}$ GeV.


Introduction
A new quasi-stable particle with mass M , spin 0, 1/2 or 1 and gravitational interactions only is a phenomenologically viable DM candidate, dubbed 'gravitational DM' [1][2][3]. Such DM can be produced through gravitational scatterings or through fluctuations during inflation.
Can gravitational DM be realised in reasonable theories, or does it need ad hoc assumptions? For example a scalar S can always have non-gravitational renormalizable quartic interactions to the Higgs H, |H| 2 |S| 2 . A singlet fermion N can always have a renormalizable Yukawa interaction to the Higgs and left-handed leptons L, N LH. The models of vectors considered in [1][2][3] have the same problem, as they involve a scalar S to make vectors massive; furthermore abelian vectors can mix with hypercharge at renormalizable level. We thereby here consider a theory based on the SM plus a new non-abelian gauge group G and no scalars or fermions charged under it. The action is where M Pl = √ 8πM Pl = 1.2 × 10 19 GeV is the Planck mass, and L SM and L DM describe the renormalizable interactions in the SM and DM sectors. As the dark sector is a pure gauge theory with non-abelian gauge group G, at renormalizable level DM is decoupled from the SM. The most generic Lagrangian is where G a µν = ∂ µ G a ν −∂ ν G a µ −g DM f abc G b µ G c ν . The dark θ DM term is physical and non-perturbatively breaks P and CP in the dark sector; we assume that θ DM is of order unity, and it will not give qualitatively new effects. Given that the dark sector confines at a scale Λ DM , the glueball hadrons made of vectors are possible DM candidates provided that they are long lived enough. Different decay rates arise depending on the possible reasonable assumptions considered below.
1. In the most extreme case one could assume that L NRO = 0 i.e. that gravity is the only nonrenormalizable interaction. This is a consistent possibility in theories of renormalizable quantum gravity based on 4-derivatives, that however contain possibly problematic states with negative classical energy (see e.g. [4,5]). Under this assumption, most dark bound states undergo gravitational decays with rates Γ ∼ m 5 /M 4 Pl . They can be DM only if long-lived, which needs their mass m to be below ∼ 100 TeV. We will find that some gauge groups G predict other bound states that are exactly stable, in this limit.
2. In more general theories, gravity and SM interactions give rise, via perturbative quantum corrections such as renormalisation group (RG) running, to Planck-suppressed nonrenormalizable operators that link the SM and DM sectors and respect the accidental symmetries of the renormalizable Lagrangian L SM + L DM . One example is dimension 6 operators such as L DM |H| 2 /M 2 Pl [6]. Such operators correct the rates of gravitational processes by order unity factors. 3. As an alternative possibility, generic Planck-suppressed non-renormalizable operators might be present in some theories of quantum gravity, such as those with lots of new states around the Planck scale (e.g. string models), and those where gravity gets strongly coupled (non-perturbative quantum gravity, possibly dominated by black holes and wormholes, is expected to violate accidental symmetries [7,8]). Such operators imply decays of generic dark bound states. Some bound states remain long-lived enough to be DM candidates, even if heavy.
The paper is structured as follows.
In section 2 we study the bound states and their lifetimes, finding, in addition to the ordinary glueballs, special longer-lived glueballs if the gauge group is G = SU(N ) and very long-lived states if the gauge group SO(2N ) at large N > ∼ 8. Such states can be DM candidates even if heavier than about 10 10 GeV.
In section 3 we compute DM production via gravitational thermal freeze in. This is maximally efficient if the reheating temperature T RH is comparable to Λ DM . Our results differ from [1][2][3] that considered gravitational production of massive vectors with 3 degrees of freedom, as our massless vectors have 2 degrees of freedom, and acquire mass through confinement. Our study is similar to [9], that however only considered gravitational freeze-in of SU(N ) ordinary glueballs. 1 In section 4 we consider gravitational production during inflation. As gravitational production via scatterings is dominated by the highest reheating temperature after inflation, inflation itself can contribute more. Ignoring possibly large but model-dependent effects (production of DM from inflaton decay or from post-inflationary inflaton oscillations, if kinematically allowed) we find that pure gravitational production during inflation is sub-leading. Indeed our (purely transverse) vectors are conformally coupled at tree level, so they are not produced in conformally flat cosmological backgrounds. The conformal symmetry is broken by the quantum running of the gauge coupling g DM and by confinement, giving rise to computable effects.
Conclusions are given in section 5.

Bound states of vectors
where C G is the quadratic Casimir of the group G. Furthermore we define d G as the dimension of the group G, so that C G = N and d G = N 2 − 1 for G = SU(N ); and C G = 2(N − 2), d G = N (N − 1)/2 for SO(N ). So the running dark gauge coupling is related to the energy scale Λ DM at which the dark sector confines, α DM (Λ DM ) ∼ 4π, by

Ordinary glueballs
Analogy with QCD computations [11][12][13] suggests that the lightest dark glueball is the state with spin 0 and quantum numbers J P C = 0 ++ that corresponds to the gauge-invariant operator Figure 1: Possible gravitational decays of dark glueballs. SM denotes any Standard Model particle, including gravitons.
and T a are the generators in the fundamental. We assume M DG ≈ 2Λ DM . Such state decays gravitationally into two gravitons, or into SM particles via one-graviton exchange (see fig. 1). Since two gravitational couplings are involved, the decay amplitude arises at order 1/M 2 Pl , and thereby gives a decay width of order Order one couplings depend on non-perturbative factors. Furthermore, generic Planck-suppressed operators would contribute giving comparable decay rates. We thereby do not go beyond the estimate. Particles with DM-like relic abundance that decay into SM states with lifetime τ in the range 3 min < ∼ τ < ∼ 10 26 s are excluded. The lower bound is obtained from BBN, the upper bound from indirect detection (see e.g. [14]) and the intermediate range is covered by CMB observations (see e.g. [15][16][17]). Thereby, dark glueball DM with lifetime given by eq. (5) is excluded in the mass range 100 TeV < ∼ M DG < ∼ 10 10 GeV, see fig. 2.
QCD computations also indicate that various heavier glueballs (pseudo-scalars, resonances with spin 2, as well as GGG states with spin 1 and possibly 3), are light enough that they cannot undergo fast gauge decays into two lightest glueballs. Most of such states (the exceptions will be discussed in the next sections) undergo gravitational decays with rates comparable to Γ DG , possibly including decays that involve one lighter glueball and that can be faster, fig. 1b. The dark gauge group in general has a topological θ DM term, that violates space-time parity P and CP at non-perturbative level, such that P-even glueballs (such as Tr GG) mix with P-odd glueballs (such as Tr GG). The Tr GG glueball cannot decay gravitationally, but can decay via NRO with similar rate.
We next discuss the special glueballs, long-lived because of group-theoretical accidental symmetries.

Long-lived SU glueballs: group charge conjugation
The U(1), SU(N ), SO(2k) and E 6 groups with symmetric Dynkin diagrams have complex representations, and complex conjugation is a Z 2 outer automorphism of the group [18]. The explicit action of C can be obtained considering any complex representation with generators T a : gauge interactions are left invariant by G C −→ −G * . So C-even vectors are associated to purely imaginary generators T a . The overall minus sign arises because, in field theory, vectors couple to C-odd currents (of scalars, J µ ∼ φ * T a ∂ µ φ, of vectors and of fermions). A U(1) vector is In a basis where each vector G a µ is either Odd or Even under C the group structure constants f abc and the symmetric tensor d abc satisfy f OOO , f OEE = 0, f OOE , f EEE = 0 and d OOO , d OEE = 0, d OOE , d EEE = 0. So G a µν has the same parity as G a µ and The d tensor (related to anomalies, in theories with fermions) is non-vanishing only for SU(N ) groups, excluding SU(2) = SO(3) and including SU(4) = SO (6). For all other groups the C-odd GGG glueball vanishes. SO(2N ) and E 6 groups admit complex representations, but d abc vanishes.
The C-odd glueball predicted by SU(N ) groups is stable under gravitational decays, because gravitational interactions respect C, so charge conjugation acting on SM particles and on dark vectors are independent symmetries. On the other hand, the C-odd glueball can decay  (2), SO(5) = Sp(4), SU(4) = SO (6), gravitationally if the action contains the non-renormalizable operator Tr(G{G, G})/M 2 Pl , which explicitly breaks C-parity. The C-odd state can decay into SM particles in the presence of dimension-8 operators such as Tr( In both cases, the C-odd ball acquires a decay rate

Long-lived SO glueballs: group parity
Longer-lived glueball states arise if the gauge group is SO(N ) with even N ≥ 8. Recalling that SO(N ) vectors can be described by an anti-symmetric matrix G ij = G a T a ij , a long-lived special glueball with mass M ≈ N Λ DM /2 is formed by N/2 vectors and corresponds to the gauge-invariant operator where Lorentz indices have been omitted, and the Pfaffian is the square root of the determinant. Such state, being built with the Levi-Civita tensor, is odd under a Z 2 parity in internal SO(N ) space. Following [19,20] we dub it O-parity, in order not to confuse it with usual parity P. We dub the O-odd glueballs as 'odd-balls'. O-parity can be concretely represented, for example, as a reflection of the first component of the fundamental representation, such that it acts on vectors as G ij → (−1) δ 1i +δ 1j G ij , leaving the SO(N ) group algebra invariant [19,20]. More abstractly, the odd-ball transforms as Pf G → O-parity is an accidental symmetry of the renormalizable gauge action, as well as of its gravitational extension. Under assumptions 1. or 2. in the introduction, odd-balls are gravitationally stable. They decay gravitationally under assumption 3., namely if the action contains the higher-dimensional operator Pf G/M N −4

Pl
, which explicitly breaks O-parity. Odd-balls decay into SM fields in the presence of operators such as |H| 2 Pf G/M N −2
Odd-balls exist for SO(2k) groups because the fundamental representation is real, so that vectors can be written as a matrix G ij , with two lower indices, that turns out to be antisymmetric. 3 The list of Lie groups and their invariant tensors in table 1 [22,23], suggests that similar states exist if G is the exceptional group F 4 or E 8 . Concerning F 4 , its G ij matrix of vectors is a sub-group of SO(26) vectors with rank 24 (ignoring Lorentz indices and derivatives) [24]. Concerning E 8 , its 248 × 248 matrix of vectors G ij has rank 240; we don't know if the bound state built with the i 1 ···i 248 tensor is long-lived or can fragment into smaller bound states built with the A ijk and S i 1 ···i 8 invariant tensors of E 8 .

Hadronization of vectors
In the next section we will consider pairs of vectors produced by collisions with E = √ s Λ DM . We here estimate the particle multiplicity in vector jets, determined by soft-radiation doubly enhanced by IR logarithmic terms that can be resummed through evolution equations [26,27]. The numerical factor is the same for all groups, as the vector/vector splitting function, proportional to the Casimir C G , gets cancelled by the one loop running of the gauge coupling in it, with beta function proportional to C G . As a result We next need to estimate the relative amount of special glueballs with respect to ordinary glueballs. It's enough to consider G = SO(N ) and its special odd-balls. We model hadronization 2 A similar highly protected accidental global symmetry was used to propose a high-quality axion model in [21]. 3 The G ij matrix of a Sp(N ) group is symmetric so that its contractions with the tensor vanish. More in general, is not a fundamental tensor of Sp(N ), as it can be written as as follows. We assume that nearby vectors start forming bound states until achieving singlets under G. We perform a MonteCarlo simulation, that adds to a list of vectors one more vector G ij with random i < j indices until singlets are possible. Whenever k > 1 vectors can be contracted with δ ij tensors, we assume that these k vectors form a glueball with k vectors, that drop out of the bound state. Whenever the bound state contains N/2 vectors with all different indices, we assume that these N/2 vectors form an odd-ball. Running numerically up to N = 16 we find that the fraction of the dark vector energy that ends up in SO(N ) odd-balls is approximately 1.2 × 0.76 N , mildly suppressed at large N . The energy fraction that ends up in ordinary glueballs made of k vectors is approximatively given by the Poisson distribution e −µ µ k−1 /(k − 1)! with µ = (N + 4)/6.

DM production from thermal scatterings
We here consider DM produced through the freeze-in mechanism by which a hot bath of high temperature SM particles occasionally results in gravitational collisions that pair-produce dark matter. The 2 → 2 annihilation processes SM SM → DM DM occur through the s-channel exchange of a graviton h µν . Since the production is driven by non-renormalizable gravitational interactions, it is dominated by the highest available energy scales after the end of inflation, i.e the reheating temperature T RH . Defining the Hubble rate during inflation as H infl , the reheating temperature is • if T RH Λ DM DM is directly produced in form of hadronic bound states; as discussed in section 3.4.
In both cases, one needs to take into account the possibility that faster decays of ordinary glueballs reheat the Universe, see section 3.5.

Thermal production rate of massless dark vectors
We expand the metric as g µν = η µν + 2h µν /M Pl , whereM Pl = M Pl / √ 8π is the reduced Planck mass. The graviton propagator with quadri-momentum k is i(η µµ η νν + η µν η νµ )/2k 2 + · · · . One graviton h µν couples as h µν T µν /M Pl where T µν ≡ 2 δS/δg µν = T µν SM + T µν DM + · · · is the usual energy-momentum tensor. In our case it is traceless, as we consider massless vectors, and we neglect masses of SM particles at temperature much above the weak scale. Tree-level s-channel exchange of one graviton generates the effective amplitude A = −iT /M 2 Pl s where T = T µν T µν .
The differential cross sections for production of the d G massless gauge bosons from scatterings of particles with negligible masses and spin S = {0, 1/2, 1} are obtained as dσ S /dt = S f |A | 2 S /16πs 2 , where S f = 1/2 and |A | 2 is summed over the polarizations of the final states and averaged over the initial states. Most cross sections can be obtained from [25] 4 or from [2] (see also [1], where some cross sections differ) where s = (p 1 + p 2 ) 2 , t = (p 1 − q 1 ) 2 , u = (p 1 − q 2 ) 2 are the usual Mandelstam variables. The non-minimal coupling ξ of the spin-0 scalar h to gravity does not contribute to dσ 0 , because it adds to the SM energy momentum tensor a term proportional to (∂ µ ∂ ν − η µν ∂ 2 )h 2 and because T µν DM is trace-less and conserved.
The interaction rate density in thermal equilibrium at temperature T is where S i = 1/2 if the initial particles are identical (real scalars and gauge bosons) and S i = 1 otherwise (fermion/anti-fermion pair); g S are the polarization degrees of freedom of the initial particles (1 for scalars, 2 for fermions and massless vectors) and K 1 is the modified Bessel function. We find The total interaction rate density is γ eq (SM SM → GG) = N 0 γ eq 0 + N 1/2 γ eq 1/2 + N 1 γ eq 1 = where N S are the number of degrees of freedom, N 0 = 4, N 1/2 = 45 and N 1 = 12 in the SM.

Freeze-in abundance of massless vectors
We here compute the freeze-in gravitational production rate of dark vectors with negligible mass, Λ DM T RH . Assuming that the big-bang suddenly started at the maximal temperature T RH and that the number abundance of dark vectors n vanishes at T RH , its later evolution is dictated by the Boltzmann equationṅ Here a dot denotes d/dt, H R =ȧ/a = 8πρ R /3/M Pl is the expansion rate, ρ R = π 2 g * T 4 /30 with g * = 106.75 is the energy density of SM radiation at temperature T . It is convenient to rewrite the Boltzmann equation in terms of Y = n/s as function of z = T RH /T , where s = 4ρ R /3T is the SM entropy density. The resulting Boltzmann equation is sH R z dY /dz = 2γ eq , solved by A similar result is found with a more realistic definition of T RH , that assumes that SM particles are progressively reheated by the energy released by some non-relativistic energy density ρ φ that decays with width Γ φ into SM particles only. ρ φ could be due to the inflaton, or to some non-relativistic unstable particle. The Boltzmann equations now arė where H = 8π(ρ φ + ρ R )/3/M Pl . The reheating temperature T RH is defined in terms of Γ φ as the temperature at which What happens at T T RH gets diluted by the entropy release. Numerical solutions give

Cosmological evolution of dark vectors
After vectors are gravitationally pair-produced with small number density Y = n/s ∼ (T RH /M Pl ) 3 and energy density ρ ∼ T RH n, three qualitatively different cosmological evolutions can take place, depending on the value of Λ DM : • For large Λ DM (we determine below how large) vectors don't thermalize and their number or energy density are too small to thermally block confinement. Thereby pair-produced vectors immediately hadronise: each vector forms N DG (T RH ) dark glueballs with mass ∼ Λ DM , with N DG estimated previously in eq. (9). This happens when the number density ∼ N DG n is smaller than Λ 3 DM , and thereby for large Λ DM > ∼ N  For smaller Λ DM the vector energy density is large enough to prevent immediate hadronization: • For very low Λ DM < ∼M Pl (T RH /M Pl ) 15/4 screening of color charges blocks confinement and reduces initial-state radiation, that we neglect: N DG ∼ 1. Vector thermalization takes place when the SM sector has temperature T therm SM ∼M Pl (T RH /M Pl ) 3 , because at this point the scattering rate among vectors (with cross sections σ ∼ g 4 DM /T 2 ) exceeds the Hubble rate. Energy conservation implies that vectors thermalize to temperature T therm DM ∼ M Pl (T RH /M Pl ) 15/4 . Next vectors cool and hadronise when their temperature falls below T DM ∼ Λ DM , acquiring thermal abundance Y DM ∼ (T DM /T SM ) 3 ∼ (T RH /M Pl ) 9/4 . We can neglect cannibalistic 3 ↔ 2 processes that occur after hadronization, as they only give a logarithmic correction to Y DM .
• The intermediate rangeM Pl (T RH /M Pl ) 15/4 < ∼ Λ DM < ∼ T 2 RH /M Pl is similar to the previous case, except that vectors start interacting through non-perturbative interactions, that lead to their confinement when the universe expanded enough that n ∼ Λ 3 DM , corresponding to a temperature of the SM sector equal to T had SM ∼ Λ DM (M Pl /T RH ). At this point vectors form glueballs with mass ∼ Λ DM and thermal energy ∼ T had DM . Their abundance gets later corrected by cannibalistic processes, that we can again neglect.

Thermal production and decay of dark bound states
Finally, we consider thermal production at T RH < ∼ Λ DM of dark bound states. As these are unstable, the production rate via inverse decays is simply given in terms of their decay width as Model-independent bounds imply that DM freeze-in due to inverse decays negligibly contributes to the DM density. Therefore the production of DM for temperatures T < ∼ Λ DM is dominated by the exponentially suppressed tail of production of vectors.

Glueball decays and DM dilution
If ordinary dark glueballs have lifetime long enough that they dominate the energy density of the Universe while decaying into SM particles and gravitons, they substantially reheat the Universe and dilute the abundance of stabler DM odd-balls. The dilution effect is sizeable if Λ DM T 2 RH /M Pl and τ DG < T U ∼ 10 10 yr. In this region glueballs decay at temperature We can approximate decays as instantaneous, finding that they reheat the Universe up to T RH ∼ Λ DM (Λ DM /M Pl ) 3/2 . Then, the DM abundance gets diluted down to Y diluted DM = Y DM D DG where the dilution factor D DG is estimated as if glueballs are not thermal, if glueballs are thermal. (21) Combining all effects above, and considering a model with G = SO(10), gives the numerical result in fig. 3, that shows that gravitational freeze-in can reproduce the DM abundance in two different regimes: either as ordinary glueballs lighter than about 100 TeV (see also [9]), or as special long-lived odd-balls, that can thereby be much heavier, up to 10 14 GeV in the figure. The intermediate mass range around 10 7 GeV is excluded by bounds on the decaying ordinary glueballs, if the DM abundance is achieved. In this intermediate range dilution is sizeable and the DM abundance is proportional to Y DM /Y DG where Y DG is the abundance of the decaying ordinary glueballs. This ratio does not depend on T RH , so that requiring that the measured DM abundance is reproduced fixes the value of Λ DM (Λ DM ≈ 10 8 GeV for N = 10).
In the next section we will show that pure gravitational production during inflation gives a negligible extra contribution to the DM density.

DM production from inflationary fluctuations
Freeze-in gravitational production during the Big Bang is dominated by the highest temperatures around T RH , which means that earlier cosmology can contribute more if higher energies were present. We assume an earlier inflation phase, during which the energy density of the Universe was dominated by a scalar field φ, the inflaton. Its energy density is ρ =φ 2 /2 + V and its pressure is p =φ 2 /2−V . The inflaton potential V (φ) allows for a prolonged rapid expansion (usually on a plateau, the so-called slow-roll phase), followed by an oscillatory phase which sets the reheating temperature T RH , where the inflaton acquires its mass m φ . DM production can arise in multiple ways: In the above list, b) and c) are particle physics processes (see [28]) that depend on the inflaton model. In particular, if the inflaton decays gravitationally into both SM and DM, c) would overwhelm the scattering contribution computed in the previous section. We assume instead that m φ is small enough that b) and c) can be neglected, and focus instead on the purely inflationary production a). During inflation, the negative pressure of the inflation (thanks to its barotropic parameter w ≈ −1) drives the exponential expansion of the Universe as dark energy with scale factor a ∝ e H infl t with a roughly constant H infl = (8πV /3M 2 Pl ) 1/2 . The Universe, even if somewhat inhomogeneous at small scales, quickly ends up being described by the FLRW metric As such, we expect inflation to homogenize the field relatively quickly. The FLRW spacetime is conformally flat, so that gravity does not couple to massless classical vectors that enjoy a conformal symmetry. At quantum level, conformal invariance is broken Figure 4: Evolution during and after inflation of the relevant scales, assuming instantaneous reheating at inflation end, when a = a RH . The black continuous curve is the horizon, and the region beyond the horizon is shaded in gray. The half-wavy line is the transition from vector to hadron modes. The special mode k * is the one that, after inflation at a = a * , re-enters the horizon while transitioning from vector to hadron.
by the RG running of g DM , which generates a non-vanishing trace for the energy-momentum tensor, T µ µ = β g DM Tr G 2 µν /2g DM , and ultimately leads to confinement and dynamical generation of the mass scale Λ DM . However, we cannot perform lattice simulations of strongly-interacting vectors during Universe expansion, and inflationary production is usually computed in terms of free modes. We thereby approximate the system by making explicit its quasi-free degrees of freedom in the following unusual effective action for modes with k/a > ∼ Λ DM X [(∂X) 2 − (m 2 X + ξ X R)X 2 ]/2 for modes with k/a < ∼ Λ DM (23) where X are the multiple quasi-stable dark hadrons with masses m ∼ Λ DM and unspecified spins 0, 1, . . .. Eq. (23) means that, when expanding fields into modes, we only retain hadrons at low energy and vectors at high-energy; for these we approximate the scale anomaly by using the running gauge coupling g DM renormalized at their energy scale. As usual when computing in cosmology, we define k as the comoving momentum, so that k/a is the momentum. If H infl < ∼ Λ DM the vector description is never relevant (modes 'confine' early during inflation), and the computation reduces to inflationary production of massive dark hadrons X, to be studied in section 4.1. If H infl > ∼ Λ DM modes cross inflation end while being in the vector description and transition to the hadron description later: fig. 4 sketches how the relevant scales evolve, and the computation is performed in section 4.2. 5

Inflationary production at H < ∼ Λ DM
In this regime inflation can produce dark hadrons X. We focus on a massive scalar hadron bound state, as massive higher-spin bound states only exist in the low-energy regime where the peculiar behaviour of massive fundamental higher-spin states (e.g. longitudinal modes of spin 1, see [30][31][32][33]) does not take place. Then the bound-state field equations in the homogenous limit areẌ + 3HẊ + (m 2 + ξR)X = 0 (24) where R = 12H 2 infl during inflation. The non-minimal coupling to gravity ξ is, in principle, computable [34]. Non-perturbative dynamics is not expected to produce the special conformal value ξ = 1/6, as strong coupling breaks the conformal symmetry. In order to examine how dark matter is produced, we perturb the field around the homogeneous value X = 0 as [35] X( Neglecting time derivatives of H infl , the equations of motion for the perturbation components X k with comoving momentum k arë The equations for X k are solved in terms of Bessel functions of order iµ. Particle production during inflation can be computed at the end of inflation by using the in-out formalism: there are two classes of solution to the equation of motion, X in k and X out k , found by imposing as a boundary condition the Bunch-Davies vacuum at very early times aH k Since a depends on time, the vacuum is time-dependent, and so the mode functions evolve according to the Bogoliubov transformation The density of the final products is given by 0 in |a † out ( p)a out ( q)|0 in = |β k | 2 δ (3) ( p − q), which means The Bogoliubov factor that relates in and out states is |β k | 2 = 1/(e 2πµ − 1). The integral in eq. (29) is divergent. This is to be expected, since it corresponds to the total density of particles produced across all time. Moreover, the bulk of particle production occurs when the two terms in ω k are equal, corresponding to a vector mode with a wavelength set by its effective mass. We may then use the time at which production reaches its peak for each mode in order to write the integral for N X in terms of η. This is given by k/(aH infl ) = µ, or at conformal time η = −µ/k, and so

Inflationary production at H Λ DM
In this regime relevant modes exit from inflation while still being described by perturbative vectors, rather than by hadrons. The vector action that effectively accounts for the quantum RGE running is eq. (23). The vector equations of motion (we can drop non-abelian terms) are The coupling g DM runs with the length scale as where d ln a = Hdt = −dN corresponds to the number of e-foldings. We thereby have 1 in the perturbative regime, and ∼ 1 around confinement. As expected, the trace anomaly disappears when the running of g DM is neglected. Eliminating the temporal componentĜ 0 along with the divergence of the spatial componentsĜ i , the classical equation for the homogeneous transverse spatial components of vectors,Ĝ T , is We expand the transverseĜ T in terms of canonically normalized modes aŝ Neglecting time derivatives of H infl , the equations of motion for the perturbation components G k with comoving momentum k arë where the value of is evaluated around horizon exit, so that eq. (4) gives ≈ 1/ ln(a * /a RH ) in the notations of fig. 4. Eq. (37) has the same form as eq. (26) for scalar fluctuations X k , and is thereby similarly solved by Bessel functions of order iµ. Indeed, in the limit = 0, eq. (37) reduces to µ 2 = −1/4, the value that corresponds to a conformally coupled state that remains in its ground state despite Universe expansion. The term coming from the RG running of g DM makes µ 2 > −1/4, bringing the vector into the regime where it does not develop a coherent vacuum expectation, and its inflationary abundance is suppressed by a multiplicative factor of (see eq. (29) in [37]). As increases, the vector abundance remains small (|β k | 2 ∼ 1/(2πµ) and n infl ∼ µ 2 H 3 infl ) as long as barely exceeds 1/2, such that µ is small. Only in the deep non-perturbative regime ≈ 1 does the vector abundance become comparable to the bound-state abundance computed in section 4.1. The two regimes therefore smoothly connect.
In conclusion, the above demonstrates that the purely inflationary production of both free vectors and bound states, roughly given by a thermal density with T ∼ H infl or less, is negligible compared to freeze-in due to gravitational scatterings, As a result, the final result remains as shown in fig. 3, where only gravitational freeze-in is included.

Conclusions
We extended the Standard Model by adding a new 'dark' non-abelian gauge interaction and nothing else, in particular no matter fields charged under the new interaction. Then the new dark sector has no renormalizable interactions to the SM and confines at a scale Λ DM . We thereby obtained a dark sector that interacts with the SM sector only gravitationally. We also allowed for Planck-suppressed non-renormalizable operators. The dark sector forms various glueball bound states that can be long-lived enough to be Dark Matter candidates. The simplest glueballs (such as the lightest state Tr [G µν G µν ] with J P C = 0 ++ ) decay gravitationally with lifetime τ ∼ M 4 Pl /Λ 5 DM and can be DM candidates only if Λ DM is below 100 TeV.
We have shown that other special glueball states are odd under accidental symmetries of group-theoretical nature. Such special states are exactly stable in the limit where Plancksuppressed operators are neglected, and otherwise have Planck-enhanced lifetimes. SU(N ) gauge theories contain special states odd under group charge-conjugation with τ ∼ M 8 Pl /Λ 9 DM . Longer lifetimes τ = (M Pl /Λ DM ) 2N −4 /Λ DM are obtained in SO(N ) theories at N > ∼ 8 because some "odd-ball" states are odd under parity in group space. Similar states might exist for the F 4 and E 8 exceptional groups.
We compute the relic abundance generated by gravitational freeze-in and by purely inflationary fluctuations, assuming that model-dependent decays and scatterings of the inflaton into DM are negligible (for example because blocked by a too small inflaton mass). In particular, inflationary production of non-abelian vectors is possible because quantum corrections break the conformal symmetry of the classical action, necessitating an unusual computation. We found that gravitational freeze-in dominates for instantaneous reheating, and that a theory of non-abelian dark vectors provides successful gravitational DM candidates, either as ordinary glueballs for Λ DM < ∼ 100 TeV, or as special heavier glueballs, For example, fig. 3 considers a SO(10) theory, showing that odd-ball DM is allowed for 10 10 GeV < ∼ Λ DM < ∼ 10 15 GeV.
Concerning phenomenology, gravitational decays of ultra-heavy relics can give signals in indirect detection experiments. When a ultra-heavy particle decays into some SM particle, electro-weak and QCD log-enhanced quantum corrections generate a shower containing all SM particles. 6 Sub-leading inflationary production can give rise to iso-curvature inhomogeneities. These effects can easily be much below current bounds.
To conclude with a curiosity, we give a positive answer to the question raised by Dyson [39]: no fundamental reason prevents observing a graviton, and thereby confirming its quantum nature. Indeed, ultra-heavy DM particles could decay into high-energy gravitons, possibly giving a graviton flux as high as allowed by bounds on extra relativistic radiation. This is detectable in principle (although planet-scale detectors are needed) by gravitational scatterings on matter with cross section σ ∼ e 2 /M 2 Pl .