MIT-CTP/5412

††thanks: Email<EMAIL_ADDRESS>††thanks: Email<EMAIL_ADDRESS>††thanks: Email<EMAIL_ADDRESS>††thanks: Email:
<EMAIL_ADDRESS>††thanks: Email<EMAIL_ADDRESS>ORCID:
0000-0002-3082-0929

# Glueballs in a Thermal Squeezeout Model

Pouya Asadi  Center for Theoretical Physics, Massachusetts Institute of
Technology,
Cambridge, MA 02139, USA.  Eric David Kramer  Jerusalem College of Technology,
Jerusalem, Israel.  Racah Institute of Physics, Hebrew University of
Jerusalem, Jerusalem 91904, Israel.  Eric Kuflik  Racah Institute of Physics,
Hebrew University of Jerusalem, Jerusalem 91904, Israel.  Tracy R. Slatyer
Center for Theoretical Physics, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA.  Juri Smirnov  Stockholm University and The Oskar
Klein Centre for Cosmoparticle Physics,
Alba Nova, 10691 Stockholm, Sweden.

###### Abstract

It has been shown that a first order confinement phase transition can
drastically change the relic dark matter abundance in confining dark sectors
with only heavy dark quarks. We study the phenomenology of one such model with
a $Z^{\prime}$ portal to Standard Model. We find that dark glueballs are long-
lived in this setup and dilute the dark matter abundance after their decay to
Standard Model. With this effect, the correct relic abundance is obtained with
dark matter masses up to $\mathcal{O}(10^{6})$ TeV. We find that while a part
of the parameter space is already ruled out by direct detection and collider
searches, there is still a broad space of viable scenarios that can be probed
by future experiments.

##  I Introduction

The particle nature of dark matter (DM) and the interactions controlling its
dynamics are two of the major unanswered questions in the current standard
model (SM) of particle physics. While the original idea of Weakly Interacting
Massive Particles (WIMPs) Zeldovich:1965gev ; Lee:1977ua ; Steigman:1984ac ;
Steigman:2012nb  remains a viable possibility to date Leane:2018kjk , many
other intriguing dynamics have been studied as well.

If the ongoing experiments fail to find any DM particles in the conventional
mass range of thermal perturbative DM models with only DM-DM annihilations,
i.e. $\mathcal{O}(100)$ TeV or lower, it is only natural to consider next-to-
minimal dynamics that can open up new parameter space for higher masses (for
heavier DM and unitary bounds from other perturbative $2\to 2$ freezeout
processes see Refs. Kim:2019udq ; Kramer:2020sbb ) . An interesting class of
next-to-minimal DM models are confining dark sector scenarios, see for
instance Refs. Gudnason:2006yj ; Alves:2009nf ; Kilic:2009mi ; Hambye:2009fg ;
Kribs:2009fy ; Alves:2010dd ; Bai:2010qg ; Feng:2011ik ; Fok:2011yc ;
Lewis:2011zb ; Frigerio:2012uc ; Buckley:2012ky ; Appelquist:2013ms ;
Bhattacharya:2013kma ; Cline:2013zca ; Boddy:2014yra ; Boddy:2014qxa ;
Hochberg:2014kqa ; Appelquist:2015yfa ; Antipin:2015xia ; Soni:2016gzf ;
Harigaya:2016nlg ; Mitridate:2017oky ; DeLuca:2018mzn ; Contino:2018crt ;
Gross:2018zha ; Geller:2018biy ; Beylin:2019gtw ; Dondi:2019olm ;
Buttazzo:2019mvl ; Landini:2020daq ; Brower:2020mab ; Baldes:2020kam ;
Contino:2020god ; Beylin:2020bsz ; Asadi:2021yml ; Asadi:2021pwo ;
Gross:2021qgx ; Arakawa:2021wgz ; Reichert:2021cvs ; Cheng:2021kjg ;
Baldes:2021aph ; Bensalem:2021qtj , and Refs. Kribs:2016cew ; Cline:2021itd
for recent reviews.

Using a novel approach, Refs. Asadi:2021yml ; Asadi:2021pwo  showed that the
dark confinement phase transition (PT) of such models can substantially affect
the DM abundance calculation, as we summarize in the next section. As a result
of this “thermal squeezeout” effect, the dark matter mass consistent with the
relic density increases to masses far beyond the perturbative unitarity bound
on thermal DM mass Griest:1989wd ; vonHarling:2014kha ; Smirnov:2019ngs .

The thermal squeezeout effect is independent of the details of the portal
between SM and the dark sector. Nonetheless, a portal is still required
between the two sectors to allow the dark glueballs to decay to SM particles.
If these glueballs were stable, they would have become the dominant fraction
of the DM in such setups Garani:2021zrr , and in the mass range where thermal
squeezeout becomes relevant would lead to an overproduction of the DM
abundance.

In this work we consider a specific model for the thermal squeezeout scenario
found in Refs. Asadi:2021yml ; Asadi:2021pwo , that includes a single flavor
of heavy dark quarks charged under a confining SU(3) gauge group. To enable a
decay of all dark glueballs to SM particles, we must either have particles
charged both under SM and dark gauge groups (e.g. see Ref. Mitridate:2017oky
), or a vector portal between the two sectors Juknevich:2009gg . Owing to its
simplicity, we will focus on the latter and use a
$\mathrm{U(1)}_{\mathrm{B-L}}$ portal in this work.

We find that, for the ranges of dark quark masses and dark confinement scales
we are studying, it is natural for the dark glueballs to be long-lived enough
to eventually dominate the energy density of the universe after the PT, giving
rise to an early matter-dominated epoch. Once the glueballs decay to SM
particles, they can inject substantial entropy into the SM thermal bath that
can further suppress the DM relic abundance. This, in turn, pushes the relic
abundance line to even higher masses than predicted in Refs. Asadi:2021yml ;
Asadi:2021pwo .

Figure 1:  The schematic evolution of different relics abundances in our
setup. We show the schematic contours for SM (red - neglecting changes in
number of relativistic degrees of freedom), dark quarks and dark baryons
(dashed and solid black lines), and dark gluons and glueballs (dashed and
solid green lines). The PT epoch is marked by the dashed orange line. We have
assumed the glueballs are long-lived enough that they take over the energy
budget of the universe at some point (marked by the dashed blue line) and
eventually decay to the SM at a lower temperature (dashed purple line). As a
result, they inject significant amount of entropy into the SM bath that
further suppresses the dark baryons abundance. Due to this injection, SM
comoving abundance dilutes more slowly during this epoch.

A schematic evolution of different relics’ abundances in our setup is shown in
Fig. 1. The abundance of dark quarks is depleted during the PT due to the
thermal squeezeout effect; the glueballs’ abundance will also be suppressed
during the confinement PT due to complicated dynamics of the PT as we will
explain in this work. The glueballs then freeze out and, if they are long-
lived enough, eventually take over the energy budget of the universe at
$T=T_{i}$. They will eventually decay to SM particles at
$T(t)=T(\tau_{0^{++}})\equiv T_{\tau_{0^{++}}}$ where $\tau_{0^{++}}$ is the
lifetime of the lightest glueball, injecting substantial amount of entropy
into the thermal bath that dilutes the DM abundance further. Our entropy
injection setup can be thought of as a special case of WIMPs Without Weakness
Asadi:2021bxp .

This entropy injection has the most extreme effect when the glueballs decay
roughly around the onset of BBN. We will focus on this scenario. We will also
study other signals of this specific setup in ongoing experiments and find
that despite large DM masses, there can still be potentially observable
signals in a range of experiments.

In Sec. II we introduce the model and discuss the squeezeout effect and
entropy injection that dilute DM abundance. Phenomenology of the model is
studied in Sec. III. We conclude in Sec. IV. In App. A we use simple
thermodynamic relations to provide an analytical approximation for the PT
timescale. A dynamical approximation for the glueball abundance during the PT
is presented in App. B. Finally, in App. C we justify why keeping only the
lightest glueball state in our analysis is a good approximation.

##  II The Model

We extend the Lagrangian of the SM with a confining SU(3), new dark quarks
$q$, and a gauged $\mathrm{U(1)}_{\mathrm{B-L}}$

$\displaystyle\mathcal{L}$ $\displaystyle\supset$
$\displaystyle\bar{q}(i\not{D}-m_{q})q-\frac{1}{4}Z^{{}^{\prime}\mu\nu}Z^{{}^{\prime}}_{\mu\nu}-\frac{1}{4}G^{\mu\nu}G_{\mu\nu},$
(1)

where $Z^{{}^{\prime}\mu\nu}$ ($G^{\mu\nu}$) is the field strength of the new
$\mathrm{U(1)}_{\mathrm{B-L}}$ gauge group (dark SU(3)).111We assume right-
handed neutrinos exist to ensure anomaly cancellation; they do not affect the
rest of our study in this work. The covariant derivative of the SM fields now
includes the interaction with the new $\mathrm{U(1)}_{\mathrm{B-L}}$ gauge
boson $Z^{\prime}$. The mass and gauge coupling of this new gauge field are
denoted by $m_{Z^{\prime}}$ and $g_{Z^{\prime}}$, respectively. We focus on
the parameter space with $m_{Z^{\prime}}\leqslant\Lambda$, where $\Lambda$ is
the dark sector’s confinement scale, to simplify the glueballs’ decay
calculation. It is shown that a vector boson portal like the $Z^{\prime}$ in
our model enables the decay of all glueballs Juknevich:2009ji ;
Juknevich:2009gg . Another intriguing scenario arises in the case that the
$Z’$ mass becomes small, which leads to long range interactions diyann .

With a $\mathrm{U(1)}_{\mathrm{B-L}}$ portal the couplings of all SM particles
to the $Z^{\prime}$ are determined. The dark quarks are vector-like fermions;
thus, the anomaly constraints are guaranteed to be satisfied. We explicitly
set the charges of dark quarks under $\mathrm{U(1)}_{\mathrm{B-L}}$ to
$q_{q}=1$. The remaining degrees of freedom are (1) the heavy dark quark mass
$m_{q}$, (2) the dark confinement scale $\Lambda$, (3) the $Z^{\prime}$ mass
$m_{Z^{\prime}}$, and (4) its gauge coupling $g_{Z^{\prime}}$. Our goal is to
explore the space of these parameters and find the parts that explain the
observed DM abundance.

Fig. 2 (top row), shows the processes that keep the dark quarks in kinetic
equilibrium with the SM and the dark gluons before the PT. The dark quarks
remain in equilibrium with SM (via the top-left diagram in Fig. 2) when the
heat exchange rate is larger than the Hubble rate:

$\displaystyle\frac{\Gamma_{h}}{H(T)}\simeq\frac{n_{f}\langle\sigma v\rangle
T/m_{q}}{g_{*}^{1/2}T^{2}/M_{\mathrm{pl}}}\gtrsim 1,$ (2)

where $\Gamma_{h}$ is the heat transfer rate between dark quarks and SM
fermion $f$, $n_{f}$ is SM fermion number density, $\langle\sigma v\rangle$ is
the relevant cross section, the $T/m_{q}$ factor is inverse of number of
interactions required for enough heat transfer Rubakov:2017xzr , $g_{*}$ is
the relativistic degrees of freedom at temperature $T$, and
$M_{\mathrm{pl}}=1.2\times 10^{19}$ is the Planck mass. Assuming
$m_{Z^{\prime}}$ is negligible compared to transferred energy in the top-left
diagram of Fig. 2, we find that at $T=\Lambda$

$\displaystyle\frac{\Gamma_{h}}{H(\Lambda)}$ $\displaystyle\approx$
$\displaystyle
10\times\alpha^{2}_{Z^{\prime}}\frac{m_{q}}{\Lambda}\frac{M_{\mathrm{pl}}}{\Lambda}$
$\displaystyle\approx$ $\displaystyle
10^{16}\times\alpha_{Z^{\prime}}^{2}\times\left(\frac{m_{q}/\Lambda}{10^{2}}\right)\times\left(\frac{\Lambda}{10^{3}~{}\mathrm{TeV}}\right)^{-1}.$

For the range of parameters studied in this work (see below), this equation
suggests that for $\alpha_{\rm Z^{\prime}}\gtrsim 10^{-8}$ the kinetic
equilibrium will be maintained at the time of confinement PT. Similarly, we
can check that, since dark gluons are massless, with an elastic scattering
rate factor that scales as $\langle\sigma v\rangle\sim T^{-2}$, they maintain
their equilibrium with dark quarks (and subsequently with SM particles) until
the PT as well.

Figure 2: Diagrams giving rise to heat exchange between different particles in
the spectrum before (top) and after (bottom) the confinement PT. One can
easily check that the top diagrams keep the dark quarks $q$ in equilibrium
with SM fermions $f$ (top left) and with dark gluons $g$ (top right). This
suggests the SM and dark gluons remain in kinetic equilibrium as well. After
the PT the abundance of all the relics in the dark sector get substantially
suppressed and heat exchange rates become much lower. The most efficient
diagrams for heat exchange after the PT are shown between the dark baryons and
SM fermions (bottom left), and dark glueballs and baryons (bottom right).
These diagrams can not maintain the kinetic equilibrium between SM and dark
glueballs after the PT, see text for further details.

The gluon abundance tracks the equilibrium value while the quarks eventually
freeze out and their freeze-out abundance can be calculated using usual
perturbative techniques Asadi:2021pwo . The process controlling their
abundance is $qq\rightarrow gg$ and the relevant cross section was summarized
in the appendix of Ref. Asadi:2021pwo .222The annihilation to SM via a
$Z^{\prime}$ is also allowed but we checked that it is sub-dominant to the
aforementioned process.

We consider heavy enough dark quarks such that (1) they freeze out before the
confinement PT, and (2) the confinement transition is of first order
Svetitsky:1982gs ; Kaczmarek:1999mm ; PhysRevD.60.034504 ; Aoki:2006we ;
Saito:2011fs . In Refs. Asadi:2021yml ; Asadi:2021pwo  it was argued that the
PT has a profound effect on the dark matter abundance. We start by summarizing
this effect.

###  II.1 Thermal squeezeout of dark matter - a summary

At the onset of the dark confinement PT, bubbles of the confined phase
nucleate inside the ambient sea of the deconfined phase. For these bubbles to
nucleate efficiently slight supercooling is required. After their nucleation,
bubbles grow since they are energetically favored. During this growth, latent
heat is released, and the vicinity of the bubble heats back to
$T\approx\Lambda$. It has been argued Witten:1984rs ; Asadi:2021pwo  that the
bubble growth rate is limited by the heat diffusion rate, which in turn is
proportional to the amount of supercooling during the PT. For the small
supercoolings that happen in our confinement PT the bubble walls reach a non-
relativistic terminal velocity. Further numerical studies are required to
improve our understanding of the bubble wall velocity.

The evolution of the PT and its effect on quarks are summarized in Fig. 3. As
the bubbles grow, they run into heavy dark quarks. An isolated quark sees an
effectively infinite energy barrier as it encounters a bubble of the confined
phase, so it gets pushed back into the deconfined phase and remains trapped
there.

Figure 3: Different stages of the confinement PT in the dark sector and its
effect on the dark quarks. Figure is taken from Refs. Asadi:2021yml ;
Asadi:2021pwo . The PT starts with bubbles of the confined phase (dark purple)
nucleating inside the ambient sea of the deconfined phase (light purple). As
the bubbles grow, they push aside isolated quarks (black dots) and eventually
sweep them inside isolated pockets of the deconfined phase after the
percolation. As each pocket contracts, it brings the trapped quarks into
interaction again. The quarks annihilate against each other into
gluons/glueballs very efficiently and only a small fraction of them survive in
the form of stable baryons (orange dots).

After bubbles percolate, the deconfined phase survives only in the form of
isolated pockets surrounded by a sea of confined phase. In Ref. Asadi:2021pwo
, we calculate the number density of quarks in these pockets and argue that
they are still very diffused after the percolation. As these pockets contract,
the quarks are brought back together and go through a second stage of
annihilation that depletes their abundance. We find that only a small fraction
of these quarks form stable baryons during this squeezeout and move into the
confined phase (the baryons are color-neutral and do not feel an infinite
potential barrier when encountering the pocket wall). The surviving baryons
eventually account for the DM observed today.

In Refs. Asadi:2021yml ; Asadi:2021pwo  it was also pointed out that the
accidental net baryon number in the pockets can give rise to a lower bound on
the survival rate from the pockets. In this work we assume this bound on the
survival factor is saturated; this allows us to find the highest viable DM
masses in this setup.

It is worth pointing out that the calculation of Ref. Asadi:2021pwo  is merely
the initiation of a closer look at the effect of a confinement PT on the DM
abundance and more careful studies are in order to obtain more precise
predictions for the DM abundance. In particular, the inferred relation between
the DM mass and confinement scale in Ref. Asadi:2021pwo  has a roughly one-
order-of-magnitude error bar due to a number of systematic uncertainties. In
this work we use the central value of this relation as a benchmark.

###  II.2 Glueballs in the early universe

In SU(3) Yang-Mills, there are twelve different glueballs (with different
$J^{PC}$ numbers) that are stable in the absence of other interactions
Juknevich:2009gg . The lightest of these glueballs, $J^{PC}=0^{++}$, has a
mass of $m_{0^{++}}\approx 7\Lambda$ (see Ref. Mathieu:2008me  and the
references therein). The mass of these states spans a range of $\sim
7-20~{}\Lambda$. In our analysis we only keep the lightest glueball
($0^{++}$); in App. C we argue that this is a good approximation.

In Fig. 2 (bottom row), we show the interactions between the dark sector and
the SM particles that could potentially bring them into kinetic equilibrium
after the PT. The glueballs’ heat exchange with SM goes through dark baryons,
whose abundance is significantly suppressed during the PT. As a result, the
glueballs do not maintain equilibrium with the rest of the bath after the PT.
Instead of calculating the full dark glueball-baryon energy exchange rate, we
will simply show that the typical glueball scatters less than once per Hubble
time, and so thermalization would not occur even if an $\mathcal{O}(1)$
fraction of the particles’ energy was transferred per scattering. We expect
that for very heavy dark baryons the energy transferred per scattering will be
parametrically small, further suppressing thermalization.

Keeping in mind that the dark baryons eventually make up the DM today, we can
write $n_{\mathrm{B}}$ as

$n_{\mathrm{B}}\sim\frac{10^{-12}}{\xi}\times\left(\frac{1~{}\mathrm{TeV}}{m_{q}}\right)\times
T^{3},$ (4)

where $\xi\leqslant 1$ is any further suppression in baryons number density
between the PT and today. Thus, at $T=\Lambda$ the scattering rate per Hubble
can be rewritten as

$\frac{n_{\mathrm{B}}\langle\sigma
v\rangle}{H(T)}\sim\frac{10^{-1}}{\xi}\times\alpha_{\mathrm{GB}}\alpha_{\mathrm{GB-B}}\times\frac{\left(1~{}\mathrm{TeV}\right)^{2}}{m_{\mathrm{DM}}\Lambda},$
(5)

where $\alpha_{\mathrm{GB}}$ ($\alpha_{\mathrm{GB-B}}$) is the coupling
between glueballs (glueballs and baryons) and is at most $\mathcal{O}(1)$. In
the upcoming sections we show that $\xi\gtrsim 10^{-4}$ is possible in our
setup and we focus on $m_{\mathrm{DM}}\gtrsim 10^{4}$ TeV and $\Lambda\gtrsim
10$ TeV; for this range of variables Eq. (5) indicates that after the PT
glueballs and SM particles are not in kinetic equilibrium.

During the PT, the gluon and the glueball gases coexist. The former can
transfer its energy and entropy to either the latter or to the SM bath. The
expansion of the universe allows the SM and the glueball baths to absorb the
released latent heat while their temperature remains roughly constant. In App.
A we provide an argument that this implies the PT takes around a few percent
of the Hubble time, which is in agreement with the numerical results of Ref.
Asadi:2021pwo .

An important quantity for the rest of this work is the relative entropy
density stored in the glueball gas after the PT

$R\equiv\frac{s_{\mathrm{GB}}}{s_{\mathrm{SM}}},$ (6)

where $s_{\mathrm{GB}}$ is the glueball and $s_{\mathrm{SM}}$ the SM entropy
density. This ratio remains constant after the PT since the glueballs and SM
are not in kinetic equilibrium. In what follows we discuss our prediction for
the viable range of $R$.

If the glueball number density after the PT is high enough, glueballs can
still maintain kinetic and chemical equilibrium amongst themselves. The
chemical equilibrium can be maintained via $3\rightarrow 2$ interactions with
large cross sections PhysRevLett.103.153201 ; Forestell:2016qhc ;
Bhatia:2020itt .

The $3\rightarrow 2$ glueball interaction being in equilibrium suggests that
their chemical potential is zero. Furthermore, if the number density is high
enough to maintain the glueball gas in kinetic equilibrium with itself, the
highest physically meaningful temperature it can have is $T\approx\Lambda$
(for higher temperatures the glueballs melt into gluons). This temperature can
be used to find an upper bound on the entropy of the glueballs. Assuming an
equilibrium distribution for $0^{++}$ gas with $m_{0^{++}}\approx 7\Lambda$ at
$T=\Lambda$ and comparing it to the SM bath entropy density at the same
temperature333In doing so we are ignoring possible percent level deviations
from $T=\Lambda$ for both the glueball bath and the SM bath., we find

$R_{\mathrm{\max}}\sim 2.5\times 10^{-4}.$ (7)

Note that this upper bound was derived with the assumption that enough
glueballs form after the PT to maintain their kinetic and chemical
equilibrium. Accurately calculating the value of $R$ would require detailed
numerical studies and is beyond the scope of this work.

In App. B we provide further dynamical arguments that imply this bound is
actually saturated during the PT. The analytic estimate we provide in App. B
suggests that the resulting glueball number density is of the order

$\displaystyle n_{\rm GB}\approx 2\cdot 10^{-6}\,T_{c}^{3}\left(\frac{M_{\rm
pl}}{T_{c}}\right)^{0.2}\,,$ (8)

which implies that the glueball number density is above the value expected
from the thermal abundance (the number density corresponding to
$R_{\mathrm{max}}$) at $T\approx T_{c}$ for $T_{c}<500\,\text{ TeV}$ and
varies very slowly with $T_{c}$. For the range of $T_{c}$ studied in this work
(see upcoming Figs. 5-8, this approximation suggests that number of glueballs
formed after the PT nearly saturates the bound in Eq. (7).

In addition, we investigate whether the number changing interactions are in
equilibrium given the above conditions. As discussed in Forestell:2016qhc ,
the cross section for the glueball $3\rightarrow 2$ interaction is given by

$\displaystyle\langle\sigma_{3\rightarrow
2}v^{2}\rangle\approx\frac{1}{(4\pi)^{3}}\left(\frac{4\pi}{N_{c}}\right)^{6}\frac{1}{m_{\rm
GB}^{5}}\,.$ (9)

Using the freezeout condition for this interaction

$\displaystyle\left(n_{\rm GB}^{\rm min}\right)^{2}\langle\sigma_{3\rightarrow
2}v^{2}\rangle=H_{\rm PT}\approx T_{c}^{2}/M_{\rm pl}\,,$ (10)

we find the threshold glueball abundance

$\displaystyle n_{\rm GB,\rm eq.}^{3\rightarrow 2}\geq
10^{2}\,T_{c}^{3}\,\sqrt{\frac{T_{c}}{M_{\rm pl}}}\,,$ (11)

above which $3\rightarrow 2$ interactions are in equilibrium. Comparing with
Eq. (8), we find that after the PT this equilibrium condition is satisfied for
$T_{c}<10^{5}\text{ TeV}$. The parameter space studied in this work is within
this range, see the upcoming Figs. 5-8. Therefore, in the parameter region we
will investigate, those number changing processes will promptly reduce the
number density of glueballs right after the PT to the thermally predicted
value, saturating $R_{\mathrm{max}}$.

All in all, we see that even if we start from $R\geqslant R_{\mathrm{max}}$
during the PT, $3\rightarrow 2$ freeze-out interactions bring $R$ down to
$R_{\mathrm{max}}$ immediately after the PT. These glueballs are then
decoupled from the SM bath and their entropy is conserved. While the analysis
of App. B suggests $R\sim R_{\mathrm{max}}$ after the PT, it overlooks many
complications of the PT and needs to be corroborated in explicit numerical
simulations. As a result, we will use $R\leqslant R_{\mathrm{max}}$ as an
input in our calculation.

We use the analytic technique developed in Ref. Forestell:2016qhc  to
calculate the asymptotic glueball abundance after the PT; the error introduced
by using the analytic approximation of Ref. Forestell:2016qhc  (as opposed to
the full numerical treatment) is subdominant to the uncertainty introduced by
the PT dynamics; see Ref. Asadi:2021pwo  for further details.

We assume that $0^{++}$, the lightest glueball, is the only one that survives
after the PT (see App. C for an explanation why this is a good approximation).
Following the calculation of Ref. Forestell:2016qhc , one can show that

$x_{0^{++},fo}^{5/2}e^{2x_{0^{++},fo}}\approx\frac{g_{*S,fo}}{180\pi}R\left(\frac{m_{0^{++}}^{4}M_{pl}\langle\sigma_{32}v^{2}\rangle}{\sqrt{\frac{4}{5}\pi^{3}g_{*,fo}}}\right)^{3/2}$
(12)

where $x_{0^{++}}=m_{0^{++}}/T_{0^{++}}$ with $m_{0^{++}}\approx 7\Lambda$
Morningstar:1999rf ; Mathieu:2008me  ($T_{0^{++}}$) being the lightest
glueball mass (temperature of the glueball bath),
$\langle\sigma_{32}v^{2}\rangle$ is the interaction rate for the $3\rightarrow
2$ number-changing process that controls glueballs’ abundance after the PT,
subscript $fo$ refers to the time of glueball freezeout after the PT, and
$g_{*(S),fo}$ is the relativistic degrees of freedom (for entropy) when the
glueball number-changing process freezes out. We can solve this equation
iteratively to find $x_{0^{++},fo}$. The asymptotic glueball abundance then
can be written as Forestell:2016qhc

$Y_{0^{++}}\approx\frac{R}{x_{0^{++},fo}}.$ (13)

In this calculation we are neglecting the finite glueball lifetime. For large
enough values of $R$, glueballs can eventually dominate the energy budget of
the universe before they decay, giving rise to an early matter-dominant epoch.
Using Eq. (13) we can calculate the energy density of the glueballs after
their freezeout, which in turn can be used to calculate the temperature
$T_{i}$ at which the universe enters an early matter-dominant epoch. Assuming
$g_{*}=g_{*S}$ at $T=T_{i}$ (which is a good approximation since the only
particle with different temperature than SM is the non-relativistic glueball),
we can show

$T_{i}=\frac{4}{3}m_{0^{++}}Y_{0^{++}}.$ (14)

Combining Eqs. (12)-(14) we can calculate $T_{i}$ as a function of $R$.

Eventually, these glueballs have to decay away to SM before BBN, see for
instance Refs. Jedamzik:2006xz ; Kawasaki:2017bqm . In our model, this decay
is facilitated via the interaction with a $Z^{\prime}$ and proceeds through a
dark quark loop (see App. C for further details of various glueballs decay
rates).

The decay channel of $0^{++}$ glueballs is shown in Fig. 4. For
$m_{Z^{\prime}}\leqslant\Lambda$, the $Z^{\prime}$ mass gives rise to less
than $10\%$ correction in $0^{++}$ lifetime; thus for simplicity we neglect
the effect of $m_{Z^{\prime}}$ and find the following simple formula for the
glueballs decay rate

$\Gamma=\frac{\alpha_{D}^{2}\alpha_{Z^{\prime}}^{2}}{4\pi}\left(\frac{1}{60}\right)^{2}\frac{m_{0^{++}}^{3}F_{0^{++}}^{2}}{m_{q}^{8}},$
(15)

where $\alpha_{D}$ ($\alpha_{Z^{\prime}}$) is the dark confining gauge
symmetry’s (new $\mathrm{U(1)}_{\mathrm{B-L}}$’s) structure constant, $m_{q}$
is the dark quark mass, and $F_{0^{++}}$ is the glueball’s decay constant; we
use Chen:2005mg ; Juknevich:2009gg

$F_{0^{++}}\approx\frac{m_{0^{++}}^{3}}{4\alpha_{D}(m_{0^{++}})},$ (16)

with $\alpha_{D}(m_{0^{++}})$ denoting the dark structure constant evaluated
at the mass of the lightest glueball.

Figure 4:  The 1-loop diagram giving rise to the $0^{++}$ glueball decay. For
the range of masses we study, the final $Z^{\prime}$s are on-shell. The heavy
quark loop gives rise to a suppression of $\left(\Lambda/m_{q}\right)^{8}$ in
the glueball’s lifetime Juknevich:2009ji ; Juknevich:2009gg ; for the heavy
quarks in our model this can naturally give rise to long-lived glueballs.

Note that for $m_{Z^{\prime}}\geqslant m_{0^{++}}/2$ the decay channel of Fig.
4 is forbidden and $0^{++}$ glueballs can only decay to SM via off-shell
$Z^{\prime}$ particles. By focusing on the range
$m_{Z^{\prime}}\leqslant\Lambda$ we avoid this extra complication.

Figure 5: Contours of the ratio of SM entropy after ($S_{f}$) and before
($S_{i}$) the glueball decay, for $\alpha_{Z^{\prime}}=10^{-2}$ (left column),
$\alpha_{Z^{\prime}}=10^{-4}$ (right columns), $R=R_{\mathrm{max}}$ (top row),
and $R=10^{-5}$ (bottom row). As $\alpha_{Z^{\prime}}$ increases, the lifetime
of the glueballs ($\tau_{0^{++}}$) decreases, which in turn leads to smaller
entropy injection for every point in the parameter space. The contours show
substantial entropy injection over a large part of the parameter space. Every
point in the red region is ruled out by constraints on the amount of glueballs
at the onset of BBN ($t=0.01$s), see the text for more details.

###  II.3 Entropy injection

The glueballs’ energy density can be calculated as a function of $R$ using
Eqs. (12)-(14). They will eventually decay to SM particles at
$T=T_{\tau_{0^{++}}}$, injecting a substantial amount of entropy into the
thermal bath that dilutes the DM abundance further. The calculation for the
amount of this entropy injection and the dilution it gives rise to can be
found in the literature, e.g. see Refs. Kolb:1990vq ; Berlin:2016gtr ;
Contino:2018crt ; Asadi:2021bxp . In the case of a single unstable glueball
the result is

$\xi\equiv\left(\frac{S_{f}}{S_{i}}\right)^{-1}=\left(1+1.65<g_{*S}^{1/3}>\left(\frac{T_{i}^{4}}{(\Gamma
M_{pl})^{2}}\right)^{1/3}\right)^{-3/4}.$ (17)

where $S_{f}$ ($S_{i}$) is the SM bath total entropy after the decay (before
the decay), $T_{i}$ is the temperature at which the universe enters the early
matter-dominant epoch, $\Gamma$ is the decay rate of the unstable particle,
and $\langle g_{*}^{1/3}\rangle$ is an average of the relativistic degrees of
freedom in the early universe (defined in the appendix of Ref. Asadi:2021bxp
).

Equation (17) clearly shows that the entropy injection is larger when the
glueballs have a larger lifetime, i.e. smaller decay rate $\Gamma$. There are,
however, bounds on the lifetime of a long-lived relic in the early universe
from BBN, see for instance Refs. Jedamzik:2006xz ; Kawasaki:2017bqm . For
simplicity, in this work we use the bound of

$\Omega_{0^{++}}\leqslant
10^{3}~{}~{}\mathrm{at}~{}t=0.1~{}\mathrm{s},~{}\&~{}\tau_{0^{++}}\leqslant
10^{-2}~{}\mathrm{s}$ (18)

(see App. C for more details). This puts a bound on the amount of entropy
injection and DM dilution that can take place in the early universe.

In Fig. 5 we show the contours of $S_{f}/S_{i}$ on a plane of $\Lambda-
m_{q}/\Lambda$ for two different values of $\alpha_{Z^{\prime}}$ and two
different values of $R$. It is evident that we can indeed have non-negligible
entropy injection for a large part of the parameter space.

To understand the shape of the contours in Fig. 5 we should keep in mind that
for a fixed $\alpha_{Z^{\prime}}$, as we go to higher $m_{q}/\Lambda$ ratios,
Eq. (15) suggests a smaller $\Gamma$, which in turns suggests more entropy
injection (see Eq. (17)). Furthermore, going to larger $\Lambda$ (with fixed
$m_{q}/\Lambda$) means $T_{i}$ increases, which again according to Eq. (17)
implies larger entropy injection as well.

Comparing the two columns of Fig. 5 we see that when the gauge coupling
decreases, a smaller part of the parameter space gives rise to glueball decays
fast enough to avoid the BBN bounds. On the other hand, this prolonged
lifetime means a larger entropy injection for a fixed point on the $\Lambda-
m_{q}/\Lambda$ plane. The figures also suggest that larger values of $R$ imply
a larger entropy injection for a fixed point in the parameter space, as
expected.

Figure 6: Contours of maximum $S_{f}/S_{i}$ (green) and the required
$\alpha_{Z^{\prime}}$ (red) for $R=R_{\mathrm{max}}$ from Eq. (7) (top) and
$R=10^{-5}$ (bottom). For each point on the plane $\alpha_{Z^{\prime}}$ is set
to the value that maximizes the glueball lifetime $\tau_{0^{++}}$ given the
BBN constraints discussed in the text. In the red region the required
$\alpha_{Z^{\prime}}>\sqrt{4\pi}$ and so our perturbative treatment breaks
down. We find substantial entropy injection that can further dilute the DM
abundance for most of the parameter space.

If instead of fixing $\alpha_{Z^{\prime}}$, we fix the glueball lifetime to
the maximum value allowed by BBN bounds Jedamzik:2006xz ; Kawasaki:2017bqm ,
we get the maximum amount of entropy injection for any point in the parameter
space. This condition can be used to calculate $\alpha_{Z^{\prime}}$ for every
point on the parameter space. In the rest of this work we focus on this
“maximum dilution” scenario to find the upper bound on how much the glueball
entropy injection can push the relic abundance line to higher masses. The
contours of constant entropy injection ratio $S_{f}/S_{i}$ in this scenario
are shown in Fig. 6 for two different values of $R$.

The $m_{q}/\Lambda$ ratio enters the entropy injection calculation only via
$\Gamma$, see Eq. (17). Thus, since in Fig. 6 we fix $\Gamma$ for any given
$\Lambda$, the $S_{f}/S_{i}$ contours are independent of $m_{q}/\Lambda$.
Notice that to keep $\Gamma$ fixed as we increase $m_{q}/\Lambda$,
$\alpha_{Z^{\prime}}$ increases as well. As a result, eventually for large
enough values of $m_{q}/\Lambda$ the required $\alpha_{Z^{\prime}}$ exceeds
$\sqrt{4\pi}$ and our perturbative treatment breaks down. This puts an upper
bound on $m_{q}/\Lambda$, denoted by the red region on the plane of Fig. 6.
The value of $\alpha_{Z^{\prime}}$ everywhere on Fig. 6 is large enough to
guarantee $Z^{\prime}$ is in kinetic equilibrium with SM before it decays.

Due to the entropy injection, DM can have masses even higher than predicted by
the squeezeout effect of Refs. Asadi:2021yml ; Asadi:2021pwo . In the upcoming
section we combine the suppression factors from the squeezeout during the PT
and the glueball decay to find the relic abundance line on our parameter
space. As mentioned earlier, in our calculation we used the central prediction
of Refs. Asadi:2021yml ; Asadi:2021pwo  for the pocket survival factor during
the squeezeout; various sources of uncertainty in modeling of the PT in Refs.
Asadi:2021yml ; Asadi:2021pwo  give rise to around an order of magnitude
uncertainty in the DM mass in each direction in the upcoming results.

##  III The Phenomenology

In the previous section we found that the glueball decay in our model can
significantly change the viable parameter space compared to Refs.
Asadi:2021yml ; Asadi:2021pwo . In this section we study the signals of our
model in various current experiments. To avoid cluttering the discussion, we
only focus on the scenario where the glueball lifetime $\tau_{0^{++}}$ is set
to its maximum allowed by BBN, see Fig. 6 and the discussions around it. In
this scenario the DM mass is pushed to the highest values possible in our
setup. We will find that while some parts of the parameter space are already
probed and ruled out, there is still viable parameter space left that can be
probed in future experiments.

There are five free parameters in our setup ($m_{q}$, $\Lambda$,
$\alpha_{Z^{\prime}}$, $m_{Z^{\prime}}$, and $R$) of which one
($\alpha_{Z^{\prime}}$) is fixed by the maximum glueball lifetime assumption
for every point in the parameter space. In Fig. 7, the bounds from different
experiments are summarized for different values of $m_{Z^{\prime}}$ and $R$ on
a plane of $\Lambda-m_{q}/\Lambda$. In what follows we review different
searches that can probe our model and explain different regions in the figure.

Figure 7: The relic abundance line (solid orange), and constraints from direct
detection (green) and collider (orange) searches on the parameter space of our
model for different values of $m_{Z^{\prime}}$ and $R$. For any point below
the dashed green line the scattering rate at the direct detection experiments
will be below the neutrino floor. Contours of constant DM mass are denoted by
dashed blue lines. In each plot the purple region corresponds to where the
universe is overclosed. The relic abundance line when the glueball decay is
neglected is denoted by the dashed orange line. We set $\alpha_{Z^{\prime}}$
(red lines) to the value that minimizes $\Gamma$ as allowed by the BBN bounds
(similar to Fig. 6). The red region is ruled out since the required
$\alpha_{Z^{\prime}}$ there exceeds $\sqrt{4\pi}$ and so we expect
perturbatively calculated rates to break down. Depending on the $Z^{\prime}$
mass, bounds on the model can change significantly. The value of $R$ can
substantially affect the final relic abundance line location; for lower values
of $R$ the relic abundance line approaches the dashed orange line. While a
large part of the parameter space is already ruled out, viable parameter space
still exists.

Before that, we should note that Fig. 7 suggests that as the value of $R$
increases, the relic abundance line gets pushed to higher DM masses, since
more entropy will be injected into the SM after the glueballs decay. For low
enough $R$ values the injected entropy can be neglected and our prediction for
the relic abundance line goes back to the prediction of Ref. Asadi:2021yml ;
Asadi:2021pwo  (the dashed orange line in Fig. 7). The case of
$R=R_{\mathrm{max}}$ corresponds to the largest amount of entropy injected,
and thus the highest DM mass consistent with today’s observed DM abundance.

###  III.1 Direct detection

Our DM particles scatter elastically and coherently off nuclei via a t-channel
exchange of $Z^{\prime}$. Direct detection experiments can probe this
scattering rate and put a bound on the DM-nucleon scattering rate. The cross
section for this process can be calculated straightforwardly

$\sigma_{\mathrm{DD}}\approx
144\pi\alpha_{Z^{\prime}}^{2}\frac{\mu_{N}^{2}}{m_{Z^{\prime}}^{4}},$ (19)

where $\mu_{N}=m_{\mathrm{DM}}m_{\mathrm{N}}/(m_{\mathrm{DM}}+m_{\mathrm{N}})$
is the reduced mass of the incoming DM particle and the target nucleon with
mass $m_{\mathrm{N}}$. In writing this equation we have assumed the
transferred momentum $q\sim\mu_{N}v_{\mathrm{DM}}$ (where $v_{\mathrm{DM}}$ is
the incoming DM velocity) is much smaller than $m_{Z^{\prime}}$; we can check
that for the allowed $Z^{\prime}$ masses due to the collider constraints
(discussed in the next section), this is always satisfied.

The rate in Eq. (19) should be compared to the bounds on the DM-nucleon cross
section. In Fig. 7 we show the current bounds from the XENON1T experiment
XENON:2018voc  (green regions), as well as the contour with the same rate as
the neutrino floor in the future XENONnT experiment XENON:2020kmp . For any
point below the neutrino floor contour (dashed green) the neutrino scattering
gives rise to a large background and novel search techniques are needed, e.g.
see Refs. Hochberg:2016ntt ; Essig:2016crl ; Kadribasic:2017obi ;
Budnik:2017sbu ; Rajendran:2017ynw ; Griffin:2018bjn ; Coskuner:2019odd ;
Blanco:2021hlm .

Depending on $R$ and $m_{Z^{\prime}}$, we find that direct detection
experiments can probe part of the relic abundance line in our model.

Note that according to Eq. (19), since $\mu_{N}$ goes to $m_{\mathrm{N}}$ in
the limit of $m_{\mathrm{DM}}\gg m_{\mathrm{N}}$, the rate will be independent
of the DM mass. As we go to higher $\Lambda$, i.e. higher $m_{Z^{\prime}}$,
the scattering rate is suppressed and we can only probe larger values of
$\alpha_{Z^{\prime}}$. This explains the shape of green regions in Fig. 7.

As discussed earlier, we only consider $m_{Z^{\prime}}\leqslant\Lambda$. From
Fig. 7 it is clear that if we had gone to higher values of the $Z^{\prime}$
mass the direct detection experiments would lose their sensitivity.

###  III.2 Indirect detection

Direct annihilation of dark sector baryons is highly suppressed Witten:1979kh
; Geller:2018biy , and takes place in a rearrangement reaction, where the
quarks in the baryon and anti-baryons are arranged to form unstable meson
states that subsequently decay to SM particles Dondi:2019olm .

The capture cross section is dominated by the large angular momentum modes,
that are limited by $\ell\approx R_{B}\times p$, where
$R_{B}=(\alpha_{D}(m_{q})\,m_{q})^{-1}$ is the baryon Bohr radius and $p$ its
momentum. As found in Refs. Mitridate:2017oky ; Geller:2018biy  the cross
section is close to the geometric bound (with an additional suppression factor
at large kinetic energies $E_{\rm kin}$) and is given by

$\displaystyle\langle\sigma_{\rm B\bar{B}}v_{\rm rel}\rangle$ $\displaystyle=$
$\displaystyle\frac{\pi R_{B}^{2}\,\langle v_{\rm rel}\rangle/2}{\sqrt{E_{\rm
kin}/E_{B}}}=\frac{\sqrt{3}\pi}{4\alpha_{D}m_{q}^{2}}$ $\displaystyle\sim$
$\displaystyle
10^{-30}\frac{\text{cm}^{3}}{\text{s}}\,\left(\frac{10^{4}\,\text{TeV}}{m_{\mathrm{DM}}}\right)^{2}\left(\frac{0.1}{\alpha_{D}}\right)\,,$

where $v_{\rm rel}$ is the dark baryons’ relative velocity, and $E_{B}$ their
binding energy.

In principle residual interactions between the baryons can lead to a
Sommerfeld enhancement Hisano:2004ds ; Hisano:2006nn ; Cirelli:2007xd ;
Arkani-Hamed:2008hhe ; Cassel:2009wt ; Slatyer:2009vg  of the annihilation. In
the regime of $\Lambda/m_{q}\ll 1$, the baryon coupling to glueballs is very
small, preventing a large Sommerfeld enhancement. However, for small enough
$m_{Z^{\prime}}$, interactions with $Z^{\prime}$ give rise to an enhancement
of the annihilation rate in Eq. (III.2) by $\alpha_{Z^{\prime}}/v_{\rm
rel}\lesssim 10^{2}-10^{3}$.

Nonetheless, the indirect detection experiments such as Fermi-LAT Fermi-
LAT:2017opo  or IceCube IceCube:2017rdn  are not sensitive enough to test
cross sections below $10^{-26}\text{cm}^{3}/\text{s}$ at such high energies.
Thus, even if Sommerfeld enhancement takes place in a concrete model, the
sensitivity of present day experiments is insufficient to probe our viable
parameter space.

However, there is still some potential for interesting indirect-detection
signals, requiring further study. The early matter-dominant epoch in our setup
gives rise to the formation of dense DM substructures at small scales. These
dense clumps could boost the DM annihilation rate today Erickcek:2011us ;
StenDelos:2019xdk . Furthermore, excitation and de-excitation of dark baryons
in dense regions could enhance the signal in indirect detection experiments.
DM particles could also accumulate in the potential well of astrophysical
objects today. This can enhance indirect detection signal or affect the
evolution of such bodies, e.g. see Refs. PhysRevD.40.3221 ; Leane:2017vag ;
Baryakhtar:2017dbj ; Leane:2021ihh . We leave a detailed study of these
effects for future works.

###  III.3 Collider bounds

Since DM candidates in our model are heavier than the unitarity bound, it is
impossible to produce them at collider facilities in the foreseeable future.
However, the $Z^{\prime}$ portal can potentially have signals at the LHC.
There are various searches looking for a heavy $Z^{\prime}$ beyond the SM at
the LHC. In Fig. 7 we use the most recent non-resonance searches at ATLAS
ATLAS:2020yat  and CMS CMS:2021ctt  to derive the collider bounds on the
model. These searches study the interference with SM Drell-Yan process with
final electrons or muons. The bounds from these models can be approximated as

$\frac{m_{Z^{\prime}}}{\sqrt{\alpha_{Z^{\prime}}}}\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle>}}{{\sim}}\;$}\mathcal{B},$
(21)

with $\mathcal{B}$ read from the results of Refs. ATLAS:2020yat ; CMS:2021ctt
. The exact number used on the right-hand side can change depending on the the
final state ($ee$, $\mu\mu$, or $ee+\mu\mu$) and is at most in the
$\mathcal{B}\approx 30$ TeV ballpark:

$m_{Z^{\prime}}\lower
3.01385pt\hbox{$\;\stackrel{{\scriptstyle\textstyle>}}{{\sim}}\;$}30\sqrt{\alpha_{Z^{\prime}}}~{}\mathrm{TeV}.$
(22)

The value of $\alpha_{Z^{\prime}}$ in Fig. 7 is set from the calculation of
the glueballs’ lifetime. In Fig. 7 we use different values of
$m_{Z^{\prime}}/\Lambda$ in each plot and use Eq. (22) to put bounds on the
model. To explain the shape of the excluded region (orange region), note that
in the figure as we go to lower values of $\Lambda$ the $Z^{\prime}$ mass
decreases as well, while as we go to higher $m_{q}/\Lambda$ values
$\alpha_{Z^{\prime}}$ increases.

###  III.4 Other signals

Our setup can potentially give rise to many other signals in different
experiments, see Ref. Asadi:2021pwo . Nonetheless, many of these signals are
suppressed in the model and range of parameters studied here. Below we briefly
comment on a few such suppressed signals that can not be detected in the
foreseeable future.

The first order confinement PT in our model can give rise to gravitational
waves (GWs). It should be noted that the GW spectrum is independent of the
heavy quark mass as the quarks are decoupled by the time of the PT and do not
affect it. The GW signals from PTs are a function of a few thermodynamic
parameters and have been studied extensively in the literature, e.g. see Refs.
Flanagan:2005yc ; Schwaller:2015tja ; Caprini:2015zlo ; Jaeckel:2016jlh ;
Katz:2016adq ; Bartolo:2016ami ; Weir:2017wfa ; Croon:2018erz ; Ellis:2018mja
; Breitbach:2018ddu ; Helmboldt:2019pan ; Caprini:2019egz ; Brower:2020mab ;
Huang:2020mso .

One of these thermodynamic quantities is the velocity of the bubbles during
the percolation. Unfortunately, for the non-relativistic velocities in our
model of the PT the signals will be strongly suppressed and will not be
detected by any of the proposed detectors in the foreseeable future
Caprini:2015zlo .

Furthermore, the residual dark QCD interactions give rise to some self-
interaction among DM particles. Nonetheless, for the heavy DM candidates in
our setup we can easily show that the self-interaction rate is far below the
current upper bounds Markevitch:2003at ; Feng:2009mn ; Buckley:2009in ;
10.1111/j.1365-2966.2012.21182.x ; 10.1093/mnrasl/sls053 ; Tulin:2017ara ;
Bondarenko:2020mpf  and do not give rise to any detectable signals.

The squeezeout mechanism opened up parameter space above the unitarity bound
on thermal DM models through inhomogeneous DM distribution during the PT. In
Ref. Asadi:2021pwo  the possibility of detecting such inhomogeneities was
studied for the range of $\Lambda$ studied here and it was shown that this
effect is too small to be detected in on-going surveys.

Another intriguing scenario arises in the case that the $Z^{\prime}$ mass
becomes small, which leads to long range interactions. This can substantially
increase the capture probability of the dark baryons in compact objects, such
as exoplanets Leane:2020wob , leading either to heat signatures in the case of
prompt $Z^{\prime}$ decay in the objects, or to new high energy signals at the
surface of those objects Leane:2021ihh ; Leane:2021tjj . A detailed study of
this signal is outside the scope of this work and is left for future studies.

Finally, the possibility of forming stable macroscopic objects has been
studied in similar setups Gross:2021qgx  as well. However, it has been shown
that such objects can only form when both $m_{q}/\Lambda$ and $m_{q}$ are much
larger than the range considered in our setup.

Figure 8: The available parameter space of the model. We use the maximum
dilution assumption and fix the relic abundance to the observed value to
calculate $\alpha_{Z^{\prime}}$ and dark matter mass as a function of
$\Lambda$ and $R$. Contours of $\alpha_{Z^{\prime}}$ (red) and dark matter
mass (blue) are shown. In the gray region we fall outside the $m_{q}/\Lambda$
window where the analysis of Ref. Asadi:2021pwo  is viable. The red region
denotes where values of $\alpha_{Z^{\prime}}$ required by the maximum dilution
assumption exceed $\sqrt{4\pi}$ and so the assumption of perturbativity is
violated. Bounds from direct detection (green) and collider searches (orange)
are shown as well, see the previous sections for further details; the dashed
green line denotes the boundary of neutrino floor in direct detection
searches; novel search techniques are required to probe the parameter space to
the right on this line in direct detection experiments. These bounds can
change as $m_{Z^{\prime}}$ varies. We find a large available parameter space.
Depending on the value of $R$, whose precise calculation would require further
extensive numerical work including lattice simulations, we can get the right
DM abundance with DM masses roughly in the $\left(10^{5},~{}10^{6}\right)$ TeV
range.

###  III.5 Summary of signals

As argued before, we have five free parameters in our setup: $m_{q}$,
$\Lambda$, $m_{Z^{\prime}}$, and $\alpha_{Z^{\prime}}$ are the parameters of
the model, while $R$ parameterizes our lack of understanding of the details of
the PT. With the “maximum dilution” assumption we can fix the value of
$\alpha_{Z^{\prime}}$ as a function of the other parameters. Furthermore, as
evident from Fig. 7 for any fixed value of $\Lambda$ and $R$ the DM mass (or
equivalently dark quark mass $m_{q}$) can be determined such that we get the
right DM relic abundance.

In Fig. 8 we show the viable parameter space on the plane of $\Lambda-R$ for
different values of $m_{Z^{\prime}}$. The maximum dilution assumption and
getting the right DM relic abundance are used to calculate $m_{q}$ and
$\alpha_{Z^{\prime}}$ as a function of $\Lambda$ and $R$. Details of the
direct detection and collider constraints are discussed in the previous
sections. These bounds strongly depend on the value of $m_{Z^{\prime}}$ as
depicted by the figure.

Figure 8 shows that as we go to higher values of $R$ the relic abundance line
moves to higher masses; this was expected since $R$ is an indicator of how
much entropy is stored in the glueball bath and can be injected into the SM
after its decay. The highest value of $R$ shown in the figure is
$R_{\mathrm{\max}}$ defined in Eq. (7); recall that we derived this bound by
assuming enough glueballs are produced such that the glueball $3\rightarrow 2$
process attains equilibrium, and that glueballs then follow a thermal
distribution with $T=\Lambda$. If we consider lower values of $R$, eventually
the injected entropy becomes so small that it can not affect the relic
abundance of dark matter and we will recover the prediction of Ref.
Asadi:2021yml ; Asadi:2021pwo  for the DM mass.

As argued earlier (see also App. B), there are some arguments suggesting
$R_{\mathrm{max}}$ is saturated after the PT. For this value of $R$, Fig. 8
suggests the DM mass is pushed to values as high as $\sim 10^{6}$ TeV.

##  IV Conclusions

In this work we provided a specific model of “thermal squeezeout”
Asadi:2021yml ; Asadi:2021pwo  with a $Z^{\prime}$ portal between the dark and
the visible sector. We pointed out that in the range of parameters where the
squeezeout effect becomes relevant, dark glueballs can naturally be long-lived
enough that they can give rise to an early matter-dominant epoch and inject
significant entropy into the SM, which can further push the viable DM mass to
higher values.

The final DM mass range is subject to uncertainty in calculating the survival
factor during the PT Asadi:2021yml ; Asadi:2021pwo , as well as uncertainty in
the amount of entropy in the glueball bath after the PT, captured by the
parameter $R$ (see Eq. (6)). These uncertainties arise from our lack of
understanding of the details of a confinement PT in the early universe, and
motivate future work on the subject.

All in all, we found that if we neglect the uncertainty in determining the
survival factor during the PT, in the most extreme case, DM masses as high as
$\sim 10^{6}$ TeV are viable in this setup. Including the uncertainties in
calculating the survival factor (see Ref. Asadi:2021pwo  for further details)
could potentially push this upper bound to around $\sim 10^{7}$ TeV.

We identified the parameter space of the model that is already ruled out by
various experiments. In particular, direct detection and collider bounds can
rule out some parts of the parameter space for sufficiently light $Z^{\prime}$
masses. The model still has a large viable parameter space that can be probed
in future experiments. The fact that the rich dynamics in minimal confining
dark sectors, such as our setup here, can open up parameter space for DM
masses far above the often-quoted unitarity bound Griest:1989wd ;
vonHarling:2014kha ; Smirnov:2019ngs  should motivate developing new searches
that can probe this part of the parameter space.

The model can also have further interesting signals whose in-depth study is
left for future works. Firstly, it could give rise to interesting signals when
a DM particle excites to other hadrons and, subsequently, radiates energy
during its de-excitation. Accumulation of the DM around various celestial
bodies could lead to novel indirect detection signals or modify the properties
and evolution of the bodies in question PhysRevD.40.3221 ; Leane:2017vag ;
Baryakhtar:2017dbj ; Leane:2021ihh . Furthermore, the early matter-dominant
epoch in our setup, which is required for significant entropy injection after
the glueballs decay, can also give rise to interesting signatures in the
matter power spectrum and indirect detection experiments Erickcek:2011us ;
StenDelos:2019xdk ; Barenboim:2021swl , as well as modifying the spectrum of
gravitational waves from other sources Cui:2018rwi ; Auclair:2019wcv ;
Gouttenoire:2019kij ; Hook:2020phx ; Chang:2021afa ; Ertas:2021xeh .

## Acknowledgments

We thank Cari Cesarotti, William Detmold, Yann Gouttenoire, Di Liu, Patrick
Meade, Ofri Telem, and Ken Van Tilburg for useful discussions. The work of PA
and TRS was supported by the U.S. Department of Energy, Office of Science,
Office of High Energy Physics, under grant Contract Number DE-SC0012567. EK
and EDK are supported by the Israel Science Foundation (grant No. 1111/17) and
by the Binational Science Foundation (grants No. 2020220). JS is supported by
the European Research Council under grant Number 742104.

##  Appendix A Thermodynamic Approach to the Phase Transition

In our previous work Asadi:2021pwo , a numerical analysis was developed for
studying how long the PT takes. Here we provide an estimate of the length of
the phase transition based purely on thermodynamic considerations. We will
show that, throughout the phase transition, entropy is approximately
conserved, except for the initial bubble formation, and for effects that are
proportional to the amount of supercooling $\epsilon\equiv(T_{c}-T)/T_{c}$.
The total entropy production will be given by $\int\\!\epsilon(x)\,dx$, where
$x$ is the phase fraction converted to confined phase. To prove that the
entropy is approximately conserved, we assume that the temperature
$1/T\equiv\left({\partial S}/{\partial E}\right)_{V}$ and that the pressure
$P\equiv T\left({\partial S}/{\partial V}\right)_{E}$ are well-defined during
the phase transition. The condition for this to be valid is that local kinetic
equilibrium be maintained. This will be the case if the scattering time
$t_{\rm sc}\equiv 1/(n\sigma v)$ be much less than the time over which the
energy density changes by an order 1 fraction, i.e. the total time for the
phase transition, since the latent heat $L$ is an order $10^{-1}$ fraction of
the standard model density $\rho_{\rm SM}$. In this case, we can assume the
first law of thermodynamics during the phase transition:

$\displaystyle dS=\frac{dE}{T}+\frac{P}{T}\,dV\;,$ (23)

for a given comoving volume $V$. We also note that the the second Friedman
equation, describing cooling and dilution due to Hubble expansion,

$\displaystyle\dot{\rho}=-3H(\rho+P)$ (24)

can be written, defining a comoving volume $V=a^{3}(t)$, in the form

$\displaystyle dE=-P\,dV$ (25)

where $E\equiv\rho V$ and $dV=3HVdt$. Here, $P$ is the internal pressure of
that gas, which means the gas is doing the maximum work allowed by its
pressure. Such maximum work processes are special and imply conservation of
entropy. Plugging Eq. (25) back into the first law (23), we indeed find that
$dS=0$ identically for any comoving volume.

The question that remains is whether (25) remains valid during the phase
transition, when the metric and stress tensor are both inhomogeneous. In fact,
we can derive a version of (25) that does not assume homogeneity. Assuming the
metric has the form $g_{\mu\nu}={\rm diag}(1,-g^{(3)}_{ij}(\vec{x},t))$, with
$g^{(3)}$ the metric of the spatial splices, we can write covariant
conservation of stress-energy $\nabla_{\mu}T^{\mu 0}=0$ as:

$\displaystyle\partial_{t}\big{[}\sqrt{g^{(3)}}\rho(\vec{x},t)\big{]}=-\partial_{t}[\sqrt{g^{(3)}}\,]P(\vec{x},t)\;.$
(26)

Where we assumed

$\displaystyle T_{\mu\nu}={\rm
diag}\Big{(}\rho(\vec{x},t)\,,\,-g^{(3)}_{ij}(\vec{x},t)P(\vec{x},t)\Big{)}$
(27)

and we used $\partial_{t}\,{\rm
det}(g^{(3)})=\dot{g}^{(3)}{}^{ij}g^{(3)}_{ij}$. If we further assume the
pressure $P(\vec{x},t)$ is homogeneous, as it should be for a first-order
phase transition, and integrate over comoving coordinates $\vec{x}$, we
recover (25).

We can estimate the amount of entropy generated by slight deviations from
equilibrium. The first such process is the initial rapid growth of the bubbles
after their nucleation in the supercooled deconfined phase, reheating it back
to temperature $T_{c}$ where equilibrium is restored. As a bubble of confined
phase within the supercooled deconfined phase expands by a volume $\Delta V$,
the entropy in that volume drops by an amount $-L\Delta V/T_{c}$, where $L$ is
the latent heat per unit volume of the phase transition. The corresponding
entropy drop per unit volume is $-L/T_{c}$. Since the confined phase is also
lower in energy, the bubble will also release an amount of energy in the form
of thermal energy, equal to $L\Delta V$. This will increase the total entropy
by an amount $L\Delta V/T$, giving a total change in entropy $\Delta
S=(L/T-L/T_{c})\Delta V\simeq\epsilon L\Delta V/T_{c}$.

A more careful calculation, assuming that the temperature rises continuously
as the energy is released, gives a factor of a half

$\displaystyle\Delta S\simeq\frac{\epsilon L\Delta V}{2T_{c}}$ (28)

where $\epsilon$ is the initial supercooling at the onset of bubble
nucleation. Of course, the above calculation assumes that latent heat has only
weak dependence on temperature away from $T_{c}$. For the remainder of the
phase transition, where the supercooling is both much smaller and slowly
varying, the calculation is much the same, and we find that if a fraction
$\delta x$ of phase is converted to confined phase, an amount $\Delta
S=\epsilon L\delta xV/T_{c}$ of entropy is produced within a comoving volume
$V$.

A formula for the entropy production that gives both these effects together is

$\displaystyle\Delta
S\simeq\frac{L}{T_{c}}\int_{0}^{1}\\!\epsilon(x)V(x)dx\;.$ (29)

if the volume $V(x)$ is assumed to vary very little during the phase
transition, it can be pulled out of the integral, giving $\Delta
S=LV\bar{\epsilon}/T_{c}$, where $\bar{\epsilon}$ is the average supercooling
during the PT. This finding is in agreement with the observation of Ref.
DeGrand:1984uq .

Using this fact, and the fact that the entropy density of a system with zero
chemical potential is given by

$\displaystyle s=\frac{\rho+P}{T}\;,$ (30)

we can calculate the ratio of initial and final volumes of a given comoving
patch by using

$\displaystyle TS_{f}$
$\displaystyle=(\rho_{f}+P_{f})V_{f}=(\rho_{i}+P_{i})V_{i}+T\Delta S$
$\displaystyle\simeq(\rho_{i}+P_{i}+\bar{\epsilon}L)V_{i}$ (31)

we get

$\displaystyle\frac{V_{f}}{V_{i}}\simeq\frac{\rho_{i}+P_{i}+\bar{\epsilon}L}{\rho_{f}+P_{f}}\;.$
(32)

In the confined phase, we can approximate that $\rho_{f}\simeq\rho_{\rm SM}$,
since the glueball density $\rho_{\rm GB}\ll\rho_{\rm SM}$. In the deconfined
phase, the density would then need to be $\rho_{i}=\rho_{\rm SM}+L$, where $L$
is the latent heat density between the two phases, computed from the lattice
to be $L\sim 1.413\,T_{c}^{4}$ Lucini:2005vg . The energy in latent heat is
roughly a third of what would be expected from the Stefan-Boltzmann law,
$\rho_{g}=\pi^{2}/30\times 2(N_{c}^{2}-1)T^{4}_{c}$, the difference being due
to dynamical effects.

It is possible that these same dynamical effects near the phase transition
cause the glueball density to be slightly higher by the same amount than its
non-interacting equilibrium value, and the gluon density to be slightly higher
such that the difference is still $L$, but we assume these effects to be
small.

For the pressure in both phases, we use the fact that the pressure is constant
throughout the phase transition. In particular, lattice computations
Giusti_2017  indicate that the pressure of the strong sector at the phase
transition is much less than the energy density, $P\ll L$. We are therefore
justified in taking $P\simeq P_{\rm SM}=\rho_{\rm SM}/3$ both before and after
the phase transition. Using this, and assuming the supercooling
$\bar{\epsilon}$ to be small Asadi:2021pwo , we find a simple expression for
the scale factor at the end of the phase transition:

$\displaystyle\frac{a_{f}}{a_{i}}\simeq\left(1+\frac{3L}{4\rho_{\rm
SM}}\right)^{1/3}\\!\\!\\!\\!=\left(1+\frac{3\cdot 1.413}{4\cdot
106.75\,\pi^{2}/30}\right)^{1/3}\\!\\!\\!\\!\simeq 1.01\;\,.$ (33)

Furthermore, assuming a deSitter type evolution, we find that the time scale
of the PT is $t_{\rm PT}\approx 0.01/H$, i.e. of order $1\%$ of the Hubble
time, which is consistent with the numerical findings in Ref. Asadi:2021pwo .

The crucial insight we gained is that the PT is a relatively slow process,
which means that processes that are not frozen out at $T=T_{c}$, and therefore
fast at the time of the PT, will be efficient in maintaining local thermal
equilibrium throughout the PT.

##  Appendix B Dynamic Approach to the Phase Transition

In addition to the pure thermodynamic discussion, we provide an approximation
that allows us to estimate the number density of glueballs that are produced
in the PT from the dynamical point of view. Given the complicated strong
dynamics during the PT, the arguments in this appendix should be understood to
yield only a rough approximation for the value of $R$, while a detailed
numerical simulation is in order for determining its exact value.

The basic picture is that glueballs are produced in $2\rightarrow 2$ reactions
from dark sector gluons, at the boundary between the two phases. The process
can be viewed as a two step reaction. First, at the wall between the
deconfined and confined phase, two gluons are converted into a glueball. This
part of the reaction is a non-perturbative process, but its rate factor scales
as $1/T_{c}^{2}$. In the second step, in order to conserve momentum, a
glueball splits into two glueball final states. As known from the effective
Lagrangian for glueballs, the coupling strength of the three glueball
interaction is large, and we only pay a price from the kinematic suppression
of producing two massive particles.

Thus this process has a reaction rate density

$\displaystyle\gamma_{\rm GB}=n_{g}^{2}\langle\sigma_{\rm GB}v\rangle$ (34)

where the reaction rate factor is given by

$\displaystyle\langle\sigma_{\rm GB}v\rangle\approx
1/T_{c}^{2}\exp{\left[-2m_{\rm GB}/T_{c}\right]}\,.$ (35)

As discussed in Ref. Asadi:2021bxp , the bubbles nucleate and grow in size
until the bubble walls collide, at a typical radius $R_{1}$, called the
percolation radius. After that, the pockets of deconfined phase shrink,
therefore, the largest relative volume that is occupied by bubble walls, of
the region that contains one bubble, is around the time of percolation.
Therefore, the most efficient glueball production takes place when the bubbles
have a size $R\sim R_{1}$, and we focus here on that phase.

The relative volume that is occupied by the walls in a region that contains
one bubble at this time is thus given by

$\displaystyle f_{V}=\frac{V_{\rm walls}}{V_{\rm bubble}}=\frac{4\pi
T_{c}^{-1}R_{1}^{2}}{\frac{4\pi}{3}\,R_{1}^{3}}=3\left(\frac{10^{4}\,T_{c}}{M_{\rm
pl}}\right)^{2/3}\,.$ (36)

here the expression for $R_{1}$ of Ref. Asadi:2021bxp  was used. The total
number of glueballs that are produced on average per volume in the PT is thus

$\displaystyle n_{\rm GB}=\gamma_{\rm GB}\,f_{V}\,t_{\rm prod}\,,$ (37)

where $t_{\rm prod}$ is the time scale over which production is efficient.
Since the glueball production process is dominated around the bubble
percolation, the production time scale is given by the time the system spends
around the percolation period $t_{\rm prod}\approx R_{1}/v_{w}$. This gives us
an estimate for the produced glueball number density of

$\displaystyle n_{\rm GB}\approx T_{c}^{3}\frac{\exp{\left[-2m_{\rm
GB}/T_{c}\right]}}{v_{w}}\,,$ (38)

where $v_{w}\approx 0.2\left(T_{c}/M_{\rm pl}\right)^{0.2}$, as shown in Ref.
Asadi:2021bxp .

We find that while the mass of the glueball leads to a suppressed production
probability, the relatively slow advance of the wall (which produces the
glueballs) partially compensates that suppression. The glueball number
densities that are achieved are comparable or larger than the thermally
expected abundances, used to derive Eq. (7).

Comparing the energy density $m_{\rm GB}n_{\rm GB}$ implied by Eq. (38) to the
latent heat per volume $L$ released in the PT, one can think of the excess
glueballs (with respect to their equilibrium value) as being produced by a
small fraction $f\sim 10^{-2}$ of the total latent heat. In chemical
equilibrium, this latent heat would be entirely released in the form of heat.
In this case, however, where the rate of $3\rightarrow 2$ interactions is fast
but finite, we can think of the PT happening in steps, where the excess
glueballs are initially produced, and subsequently eliminated by $3\rightarrow
2$ interactions, driving the glueball number density to its equilibrium value
and releasing heat. Comparing the rate of $3\to 2$ (at densities at or above
the equililbrium value $n_{\rm GB}\simeq(m_{\rm
GB}T_{c}/2\pi)^{3/2}\exp(-m_{\rm GB}/T_{c})\,$) to the total time for the PT,
$\mathcal{O}(10^{-2}-10^{-1})H^{-1}$, shows that the rate of $3\rightarrow 2$
is indeed fast, and that we can therefore assume the equilibrium value for
$n_{\rm GB}$.

##  Appendix C Heavier Glueballs

Throughout this work, we only kept the lightest glueball in our treatment and
claimed heavier ones do not affect our analysis in a significant way. In this
appendix we justify this claim.

In Ref. Forestell:2016qhc  it is argued that immediately after the PT, if
enough glueballs are produced so that various number changing processes are in
equilibrium, the abundance of heavier CP-even glueballs deplete very rapidly.
To see this, we should keep in mind that the heavier CP-even glueballs
coannihilate via $2\rightarrow 2$ processes to the $0^{++}$ state, while the
$0^{++}$ abundance goes through freezeout via $3\rightarrow 2$ processes with
longer timescales.

The CP-odd glueballs also have a far lower abundance compared to the $0^{++}$
glueball after the PT Forestell:2016qhc . As a result, before their decay, the
lightest glueballs ($0^{++}$) dominate the dynamics of the glueball bath such
that we can study a simplified model with only these glueballs in the
spectrum.

Neglecting the $Z^{\prime}$ mass, we can use the results of Ref.
Juknevich:2009gg  to show that in our model

$\displaystyle\tau_{0^{++}}\approx 1.766\tau_{0^{-+}}\approx
0.040\tau_{2^{++}}\approx 0.004\tau_{2^{-+}},$ (39)

while other glueballs decay to one of these glueballs with a far shorter
lifetimes when $\alpha_{Z^{\prime}}$ is perturbative. Combined with the above
argument, we find that we can only keep $0^{++}$ and $2^{\pm+}$ (that are more
long-lieved than $0^{++}$) throughout our study.

When glueballs decay to SM, the entropy they inject into SM can be calculated
as

$S=\int\frac{dQ}{T},$ (40)

where $dQ$ is the heat they inject into SM during an infinitesimal time
interval and $T$ is the temperature at which this injection takes place. Since
the $0^{++}$ states have a much higher energy density than $2^{\pm+}$ ones,
the only way the latter can have significant entropy injection into SM is that
the injection takes place at lower temperatures. Equation (39) suggests that
indeed $2^{\pm+}$ entropy injection happens at $T_{\tau_{2^{\pm+}}}\sim
10^{-1}T_{\tau_{GB}}$; nonetheless, since the total energy $Q$ stored in
$2^{\pm+}$ states is much lower than the $0^{++}$ states Forestell:2016qhc ,
$T_{\tau_{2^{\pm+}}}\sim 10^{-1}T_{\tau_{GB}}$ is not low enough to change the
total injected entropy into SM perceptibly. Hence, keeping only the $0^{++}$
state is a good approximation for calculating the injected entropy into the SM
bath.

The only relevant effect of the heavier glueballs is that if they decay too
late, they can be in conflict with stringent bounds from BBN Jedamzik:2006xz ;
Kawasaki:2017bqm . These bounds are avoided by demanding that all the
glueballs decay before the onset of BBN, i.e. $t\sim 1$s. Since the longest-
living glueball is $2^{-+}$, this translates into $\tau_{2^{-+}}\lesssim 1$s;
when combined with Eq. (39), this suggests

$\tau_{0^{++}}\lesssim 10^{-2}~{}\mathrm{s}.$ (41)

We use this upper bound on the lightest glueball lifetime in our analysis. A
more careful treatment can allow for $\mathcal{O}(1)$ change in this bound,
which in turn merely gives rise to $\mathcal{O}(1)$ change in the dilution
factor from the entropy injection. This is sub-dominant to various sources of
uncertainty in our calculation of DM relic abundance during the PT
Asadi:2021pwo  and is neglected in this work.

## References

  * (1) Y. b. Zeldovich, “Survey of Modern Cosmology,” Adv. Astron. Astrophys. 3 (1965) 241–379.
  * (2) B. W. Lee and S. Weinberg, “Cosmological Lower Bound on Heavy Neutrino Masses,” Phys. Rev. Lett. 39 (1977) 165–168.
  * (3) G. Steigman and M. S. Turner, “Cosmological Constraints on the Properties of Weakly Interacting Massive Particles,” Nucl. Phys. B 253 (1985) 375–386.
  * (4) G. Steigman, B. Dasgupta, and J. F. Beacom, “Precise Relic WIMP Abundance and its Impact on Searches for Dark Matter Annihilation,” Phys. Rev. D 86 (2012) 023506, arXiv:1204.3622 [hep-ph].
  * (5) R. K. Leane, T. R. Slatyer, J. F. Beacom, and K. C. Y. Ng, “GeV-scale thermal WIMPs: Not even slightly ruled out,” Phys. Rev. D 98 no. 2, (2018) 023016, arXiv:1805.10305 [hep-ph].
  * (6) H. Kim and E. Kuflik, “Superheavy Thermal Dark Matter,” Phys. Rev. Lett. 123 no. 19, (2019) 191801, arXiv:1906.00981 [hep-ph].
  * (7) E. D. Kramer, E. Kuflik, N. Levi, N. J. Outmezguine, and J. T. Ruderman, “Heavy Thermal Dark Matter from a New Collision Mechanism,” Phys. Rev. Lett. 126 no. 8, (2021) 081802, arXiv:2003.04900 [hep-ph].
  * (8) S. B. Gudnason, C. Kouvaris, and F. Sannino, “Dark Matter from new Technicolor Theories,” Phys. Rev. D 74 (2006) 095008, arXiv:hep-ph/0608055.
  * (9) D. S. M. Alves, S. R. Behbahani, P. Schuster, and J. G. Wacker, “Composite Inelastic Dark Matter,” Phys. Lett. B 692 (2010) 323–326, arXiv:0903.3945 [hep-ph].
  * (10) C. Kilic, T. Okui, and R. Sundrum, “Vectorlike Confinement at the LHC,” JHEP 02 (2010) 018, arXiv:0906.0577 [hep-ph].
  * (11) T. Hambye and M. H. G. Tytgat, “Confined hidden vector dark matter,” Phys. Lett. B 683 (2010) 39–41, arXiv:0907.1007 [hep-ph].
  * (12) G. D. Kribs, T. S. Roy, J. Terning, and K. M. Zurek, “Quirky Composite Dark Matter,” Phys. Rev. D 81 (2010) 095001, arXiv:0909.2034 [hep-ph].
  * (13) D. Spier Moreira Alves, S. R. Behbahani, P. Schuster, and J. G. Wacker, “The Cosmology of Composite Inelastic Dark Matter,” JHEP 06 (2010) 113, arXiv:1003.4729 [hep-ph].
  * (14) Y. Bai and R. J. Hill, “Weakly Interacting Stable Pions,” Phys. Rev. D 82 (2010) 111701, arXiv:1005.0008 [hep-ph].
  * (15) J. L. Feng and Y. Shadmi, “WIMPless Dark Matter from Non-Abelian Hidden Sectors with Anomaly-Mediated Supersymmetry Breaking,” Phys. Rev. D 83 (2011) 095011, arXiv:1102.0282 [hep-ph].
  * (16) R. Fok and G. D. Kribs, “Chiral Quirkonium Decays,” Phys. Rev. D 84 (2011) 035001, arXiv:1106.3101 [hep-ph].
  * (17) R. Lewis, C. Pica, and F. Sannino, “Light Asymmetric Dark Matter on the Lattice: SU(2) Technicolor with Two Fundamental Flavors,” Phys. Rev. D 85 (2012) 014504, arXiv:1109.3513 [hep-ph].
  * (18) M. Frigerio, A. Pomarol, F. Riva, and A. Urbano, “Composite Scalar Dark Matter,” JHEP 07 (2012) 015, arXiv:1204.2808 [hep-ph].
  * (19) M. R. Buckley and E. T. Neil, “Thermal dark matter from a confining sector,” Phys. Rev. D 87 no. 4, (2013) 043510, arXiv:1209.6054 [hep-ph].
  * (20) Lattice Strong Dynamics (LSD) Collaboration, T. Appelquist et al., “Lattice Calculation of Composite Dark Matter Form Factors,” Phys. Rev. D 88 no. 1, (2013) 014502, arXiv:1301.1693 [hep-ph].
  * (21) S. Bhattacharya, B. Melić, and J. Wudka, “Pionic Dark Matter,” JHEP 02 (2014) 115, arXiv:1307.2647 [hep-ph].
  * (22) J. M. Cline, Z. Liu, G. Moore, and W. Xue, “Composite strongly interacting dark matter,” Phys. Rev. D 90 no. 1, (2014) 015023, arXiv:1312.3325 [hep-ph].
  * (23) K. K. Boddy, J. L. Feng, M. Kaplinghat, and T. M. P. Tait, “Self-Interacting Dark Matter from a Non-Abelian Hidden Sector,” Phys. Rev. D 89 no. 11, (2014) 115017, arXiv:1402.3629 [hep-ph].
  * (24) K. K. Boddy, J. L. Feng, M. Kaplinghat, Y. Shadmi, and T. M. P. Tait, “Strongly interacting dark matter: Self-interactions and keV lines,” Phys. Rev. D 90 no. 9, (2014) 095016, arXiv:1408.6532 [hep-ph].
  * (25) Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky, and J. G. Wacker, “Model for Thermal Relic Dark Matter of Strongly Interacting Massive Particles,” Phys. Rev. Lett. 115 no. 2, (2015) 021301, arXiv:1411.3727 [hep-ph].
  * (26) T. Appelquist et al., “Stealth Dark Matter: Dark scalar baryons through the Higgs portal,” Phys. Rev. D 92 no. 7, (2015) 075030, arXiv:1503.04203 [hep-ph].
  * (27) O. Antipin, M. Redi, A. Strumia, and E. Vigiani, “Accidental Composite Dark Matter,” JHEP 07 (2015) 039, arXiv:1503.08749 [hep-ph].
  * (28) A. Soni and Y. Zhang, “Hidden SU(N) Glueball Dark Matter,” Phys. Rev. D 93 no. 11, (2016) 115025, arXiv:1602.00714 [hep-ph].
  * (29) K. Harigaya, M. Ibe, K. Kaneta, W. Nakano, and M. Suzuki, “Thermal Relic Dark Matter Beyond the Unitarity Limit,” JHEP 08 (2016) 151, arXiv:1606.00159 [hep-ph].
  * (30) A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, “Dark Matter as a weakly coupled Dark Baryon,” JHEP 10 (2017) 210, arXiv:1707.05380 [hep-ph].
  * (31) V. De Luca, A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, “Colored Dark Matter,” Phys. Rev. D 97 no. 11, (2018) 115024, arXiv:1801.01135 [hep-ph].
  * (32) R. Contino, A. Mitridate, A. Podo, and M. Redi, “Gluequark Dark Matter,” JHEP 02 (2019) 187, arXiv:1811.06975 [hep-ph].
  * (33) C. Gross, A. Mitridate, M. Redi, J. Smirnov, and A. Strumia, “Cosmological Abundance of Colored Relics,” Phys. Rev. D 99 no. 1, (2019) 016024, arXiv:1811.08418 [hep-ph].
  * (34) M. Geller, S. Iwamoto, G. Lee, Y. Shadmi, and O. Telem, “Dark quarkonium formation in the early universe,” JHEP 06 (2018) 135, arXiv:1802.07720 [hep-ph].
  * (35) V. Beylin, M. Y. Khlopov, V. Kuksa, and N. Volchanskiy, “Hadronic and Hadron-Like Physics of Dark Matter,” Symmetry 11 no. 4, (2019) 587, arXiv:1904.12013 [hep-ph].
  * (36) N. A. Dondi, F. Sannino, and J. Smirnov, “Thermal history of composite dark matter,” Phys. Rev. D 101 no. 10, (2020) 103010, arXiv:1905.08810 [hep-ph].
  * (37) D. Buttazzo, L. Di Luzio, P. Ghorbani, C. Gross, G. Landini, A. Strumia, D. Teresi, and J.-W. Wang, “Scalar gauge dynamics and Dark Matter,” JHEP 01 (2020) 130, arXiv:1911.04502 [hep-ph].
  * (38) G. Landini and J.-W. Wang, “Dark Matter in scalar Sp($\mathcal{N}$) gauge dynamics,” JHEP 06 (2020) 167, arXiv:2004.03299 [hep-ph].
  * (39) Lattice Strong Dynamics Collaboration, R. C. Brower et al., “Stealth dark matter confinement transition and gravitational waves,” Phys. Rev. D 103 no. 1, (2021) 014505, arXiv:2006.16429 [hep-lat].
  * (40) I. Baldes, Y. Gouttenoire, and F. Sala, “String Fragmentation in Supercooled Confinement and implications for Dark Matter,” arXiv:2007.08440 [hep-ph].
  * (41) R. Contino, A. Podo, and F. Revello, “Composite Dark Matter from Strongly-Interacting Chiral Dynamics,” arXiv:2008.10607 [hep-ph].
  * (42) V. Beylin, M. Khlopov, V. Kuksa, and N. Volchanskiy, “New physics of strong interaction and Dark Universe,” Universe 6 no. 11, (2020) 196, arXiv:2010.13678 [hep-ph].
  * (43) P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer, and J. Smirnov, “Accidentally Asymmetric Dark Matter,” arXiv:2103.09822 [hep-ph].
  * (44) P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer, and J. Smirnov, “Thermal Squeezeout of Dark Matter,” arXiv:2103.09827 [hep-ph].
  * (45) C. Gross, G. Landini, A. Strumia, and D. Teresi, “Dark Matter as dark dwarfs and other macroscopic objects: multiverse relics?,” JHEP 09 (2021) 033, arXiv:2105.02840 [hep-ph].
  * (46) J. Arakawa, A. Rajaraman, and T. M. P. Tait, “Annihilogenesis,” arXiv:2109.13941 [hep-ph].
  * (47) M. Reichert, F. Sannino, Z.-W. Wang, and C. Zhang, “Dark Confinement and Chiral Phase Transitions: Gravitational Waves vs Matter Representations,” arXiv:2109.11552 [hep-ph].
  * (48) H.-C. Cheng, L. Li, and E. Salvioni, “A Theory of Dark Pions,” arXiv:2110.10691 [hep-ph].
  * (49) I. Baldes, Y. Gouttenoire, F. Sala, and G. Servant, “Supercool Composite Dark Matter beyond 100 TeV,” arXiv:2110.13926 [hep-ph].
  * (50) W. Bensalem and D. Stolarski, “Flavor and $CP$ Violation from a QCD-like Hidden Sector,” arXiv:2111.05515 [hep-ph].
  * (51) G. D. Kribs and E. T. Neil, “Review of strongly-coupled composite dark matter models and lattice simulations,” Int. J. Mod. Phys. A 31 no. 22, (2016) 1643004, arXiv:1604.04627 [hep-ph].
  * (52) J. M. Cline, “Dark atoms and composite dark matter,” in Les Houches summer school on Dark Matter.  8, 2021.  arXiv:2108.10314 [hep-ph].
  * (53) K. Griest and M. Kamionkowski, “Unitarity Limits on the Mass and Radius of Dark Matter Particles,” Phys. Rev. Lett. 64 (1990) 615.
  * (54) B. von Harling and K. Petraki, “Bound-state formation for thermal relic dark matter and unitarity,” JCAP 12 (2014) 033, arXiv:1407.7874 [hep-ph].
  * (55) J. Smirnov and J. F. Beacom, “TeV-Scale Thermal WIMPs: Unitarity and its Consequences,” Phys. Rev. D 100 no. 4, (2019) 043029, arXiv:1904.11503 [hep-ph].
  * (56) R. Garani, M. Redi, and A. Tesi, “Dark QCD matters,” JHEP 12 (2021) 139, arXiv:2105.03429 [hep-ph].
  * (57) J. E. Juknevich, “Pure-glue hidden valleys through the Higgs portal,” JHEP 08 (2010) 121, arXiv:0911.5616 [hep-ph].
  * (58) P. Asadi, T. R. Slatyer, and J. Smirnov, “WIMPs Without Weakness: Generalized Mass Window with Entropy Injection,” arXiv:2111.11444 [hep-ph].
  * (59) J. E. Juknevich, D. Melnikov, and M. J. Strassler, “A Pure-Glue Hidden Valley I. States and Decays,” JHEP 07 (2009) 055, arXiv:0903.0883 [hep-ph].
  * (60) D. Liu, Y. Gouttenoire, and E. Kuflik, “Thermal squeezeout of charged relics,” work in progress.
  * (61) V. A. Rubakov and D. S. Gorbunov, Introduction to the Theory of the Early Universe: Hot big bang theory.  World Scientific, Singapore, 2017.
  * (62) B. Svetitsky and L. G. Yaffe, “Critical behavior at finite-temperature confinement transitions,” Nuclear Physics B 210 no. 4, (1982) 423 – 447. http://www.sciencedirect.com/science/article/pii/0550321382901729.
  * (63) O. Kaczmarek, F. Karsch, E. Laermann, and M. Lutgemeier, “Heavy quark potentials in quenched QCD at high temperature,” Phys. Rev. D 62 (2000) 034021, arXiv:hep-lat/9908010.
  * (64) C. Alexandrou, A. Boriçi, A. Feo, P. de Forcrand, A. Galli, F. Jegerlehner, and T. Takaishi, “Deconfinement phase transition in one-flavor qcd,” Phys. Rev. D 60 (Jun, 1999) 034504. https://link.aps.org/doi/10.1103/PhysRevD.60.034504.
  * (65) Y. Aoki, G. Endrodi, Z. Fodor, S. Katz, and K. Szabo, “The Order of the quantum chromodynamics transition predicted by the standard model of particle physics,” Nature 443 (2006) 675–678, arXiv:hep-lat/0611014.
  * (66) WHOT-QCD Collaboration, H. Saito, S. Ejiri, S. Aoki, T. Hatsuda, K. Kanaya, Y. Maezawa, H. Ohno, and T. Umeda, “Phase structure of finite temperature QCD in the heavy quark region,” Phys. Rev. D 84 (2011) 054502, arXiv:1106.0974 [hep-lat]. [Erratum: Phys.Rev.D 85, 079902 (2012)].
  * (67) E. Witten, “Cosmic Separation of Phases,” Phys. Rev. D 30 (1984) 272–285.
  * (68) V. Mathieu, N. Kochelev, and V. Vento, “The Physics of Glueballs,” Int. J. Mod. Phys. E 18 (2009) 1–49, arXiv:0810.4453 [hep-ph].
  * (69) N. P. Mehta, S. T. Rittenhouse, J. P. D’Incao, J. von Stecher, and C. H. Greene, “General theoretical description of $n$-body recombination,” Phys. Rev. Lett. 103 (Oct, 2009) 153201. https://link.aps.org/doi/10.1103/PhysRevLett.103.153201.
  * (70) L. Forestell, D. E. Morrissey, and K. Sigurdson, “Non-Abelian Dark Forces and the Relic Densities of Dark Glueballs,” Phys. Rev. D 95 no. 1, (2017) 015032, arXiv:1605.08048 [hep-ph].
  * (71) D. Bhatia and S. Mukhopadhyay, “Unitarity limits on thermal dark matter in (non-)standard cosmologies,” JHEP 03 (2021) 133, arXiv:2010.09762 [hep-ph].
  * (72) C. J. Morningstar and M. J. Peardon, “The Glueball spectrum from an anisotropic lattice study,” Phys. Rev. D 60 (1999) 034509, arXiv:hep-lat/9901004.
  * (73) K. Jedamzik, “Big bang nucleosynthesis constraints on hadronically and electromagnetically decaying relic neutral particles,” Phys. Rev. D 74 (2006) 103509, arXiv:hep-ph/0604251.
  * (74) M. Kawasaki, K. Kohri, T. Moroi, and Y. Takaesu, “Revisiting Big-Bang Nucleosynthesis Constraints on Long-Lived Decaying Particles,” Phys. Rev. D 97 no. 2, (2018) 023502, arXiv:1709.01211 [hep-ph].
  * (75) Y. Chen et al., “Glueball spectrum and matrix elements on anisotropic lattices,” Phys. Rev. D 73 (2006) 014516, arXiv:hep-lat/0510074.
  * (76) E. W. Kolb and M. S. Turner, “The Early Universe,” Front. Phys. 69 (1990) 1–547.
  * (77) A. Berlin, D. Hooper, and G. Krnjaic, “Thermal Dark Matter From A Highly Decoupled Sector,” Phys. Rev. D 94 no. 9, (2016) 095019, arXiv:1609.02555 [hep-ph].
  * (78) XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 no. 11, (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * (79) XENON Collaboration, E. Aprile et al., “Projected WIMP sensitivity of the XENONnT dark matter experiment,” JCAP 11 (2020) 031, arXiv:2007.08796 [physics.ins-det].
  * (80) Y. Hochberg, Y. Kahn, M. Lisanti, C. G. Tully, and K. M. Zurek, “Directional detection of dark matter with two-dimensional targets,” Phys. Lett. B 772 (2017) 239–246, arXiv:1606.08849 [hep-ph].
  * (81) R. Essig, J. Mardon, O. Slone, and T. Volansky, “Detection of sub-GeV Dark Matter and Solar Neutrinos via Chemical-Bond Breaking,” Phys. Rev. D 95 no. 5, (2017) 056011, arXiv:1608.02940 [hep-ph].
  * (82) F. Kadribasic, N. Mirabolfathi, K. Nordlund, A. E. Sand, E. Holmström, and F. Djurabekova, “Directional Sensitivity In Light-Mass Dark Matter Searches With Single-Electron Resolution Ionization Detectors,” Phys. Rev. Lett. 120 no. 11, (2018) 111301, arXiv:1703.05371 [physics.ins-det].
  * (83) R. Budnik, O. Chesnovsky, O. Slone, and T. Volansky, “Direct Detection of Light Dark Matter and Solar Neutrinos via Color Center Production in Crystals,” Phys. Lett. B 782 (2018) 242–250, arXiv:1705.03016 [hep-ph].
  * (84) S. Rajendran, N. Zobrist, A. O. Sushkov, R. Walsworth, and M. Lukin, “A method for directional detection of dark matter using spectroscopy of crystal defects,” Phys. Rev. D 96 no. 3, (2017) 035009, arXiv:1705.09760 [hep-ph].
  * (85) S. Griffin, S. Knapen, T. Lin, and K. M. Zurek, “Directional Detection of Light Dark Matter with Polar Materials,” Phys. Rev. D 98 no. 11, (2018) 115034, arXiv:1807.10291 [hep-ph].
  * (86) A. Coskuner, A. Mitridate, A. Olivares, and K. M. Zurek, “Directional Dark Matter Detection in Anisotropic Dirac Materials,” Phys. Rev. D 103 no. 1, (2021) 016006, arXiv:1909.09170 [hep-ph].
  * (87) C. Blanco, Y. Kahn, B. Lillard, and S. D. McDermott, “Dark Matter Daily Modulation With Anisotropic Organic Crystals,” Phys. Rev. D 104 (2021) 036011, arXiv:2103.08601 [hep-ph].
  * (88) E. Witten, “Baryons in the 1/n Expansion,” Nucl. Phys. B 160 (1979) 57–115.
  * (89) J. Hisano, S. Matsumoto, M. M. Nojiri, and O. Saito, “Non-perturbative effect on dark matter annihilation and gamma ray signature from galactic center,” Phys. Rev. D 71 (2005) 063528, arXiv:hep-ph/0412403.
  * (90) J. Hisano, S. Matsumoto, M. Nagai, O. Saito, and M. Senami, “Non-perturbative effect on thermal relic abundance of dark matter,” Phys. Lett. B 646 (2007) 34–38, arXiv:hep-ph/0610249.
  * (91) M. Cirelli, A. Strumia, and M. Tamburini, “Cosmology and Astrophysics of Minimal Dark Matter,” Nucl. Phys. B 787 (2007) 152–175, arXiv:0706.4071 [hep-ph].
  * (92) N. Arkani-Hamed, D. P. Finkbeiner, T. R. Slatyer, and N. Weiner, “A Theory of Dark Matter,” Phys. Rev. D 79 (2009) 015014, arXiv:0810.0713 [hep-ph].
  * (93) S. Cassel, “Sommerfeld factor for arbitrary partial wave processes,” J. Phys. G 37 (2010) 105009, arXiv:0903.5307 [hep-ph].
  * (94) T. R. Slatyer, “The Sommerfeld enhancement for dark matter with an excited state,” JCAP 02 (2010) 028, arXiv:0910.5713 [hep-ph].
  * (95) Fermi-LAT Collaboration, M. Ackermann et al., “The Fermi Galactic Center GeV Excess and Implications for Dark Matter,” Astrophys. J. 840 no. 1, (2017) 43, arXiv:1704.03910 [astro-ph.HE].
  * (96) IceCube Collaboration, M. G. Aartsen et al., “Search for Neutrinos from Dark Matter Self-Annihilations in the center of the Milky Way with 3 years of IceCube/DeepCore,” Eur. Phys. J. C 77 no. 9, (2017) 627, arXiv:1705.08103 [hep-ex].
  * (97) A. L. Erickcek and K. Sigurdson, “Reheating Effects in the Matter Power Spectrum and Implications for Substructure,” Phys. Rev. D 84 (2011) 083503, arXiv:1106.0536 [astro-ph.CO].
  * (98) M. Sten Delos, T. Linden, and A. L. Erickcek, “Breaking a dark degeneracy: The gamma-ray signature of early matter domination,” Phys. Rev. D 100 no. 12, (2019) 123546, arXiv:1910.08553 [astro-ph.CO].
  * (99) I. Goldman and S. Nussinov, “Weakly interacting massive particles and neutron stars,” Phys. Rev. D 40 (Nov, 1989) 3221–3230. https://link.aps.org/doi/10.1103/PhysRevD.40.3221.
  * (100) R. K. Leane, K. C. Y. Ng, and J. F. Beacom, “Powerful Solar Signatures of Long-Lived Dark Mediators,” Phys. Rev. D 95 no. 12, (2017) 123016, arXiv:1703.04629 [astro-ph.HE].
  * (101) M. Baryakhtar, J. Bramante, S. W. Li, T. Linden, and N. Raj, “Dark Kinetic Heating of Neutron Stars and An Infrared Window On WIMPs, SIMPs, and Pure Higgsinos,” Phys. Rev. Lett. 119 no. 13, (2017) 131801, arXiv:1704.01577 [hep-ph].
  * (102) R. K. Leane, T. Linden, P. Mukhopadhyay, and N. Toro, “Celestial-Body Focused Dark Matter Annihilation Throughout the Galaxy,” Phys. Rev. D 103 no. 7, (2021) 075030, arXiv:2101.12213 [astro-ph.HE].
  * (103) ATLAS Collaboration, G. Aad et al., “Search for new non-resonant phenomena in high-mass dilepton final states with the ATLAS detector,” JHEP 11 (2020) 005, arXiv:2006.12946 [hep-ex]. [Erratum: JHEP 04, 142 (2021)].
  * (104) CMS Collaboration, A. M. Sirunyan et al., “Search for resonant and nonresonant new phenomena in high-mass dilepton final states at $\sqrt{s}$ = 13 TeV,” JHEP 07 (2021) 208, arXiv:2103.02708 [hep-ex].
  * (105) E. E. Flanagan and S. A. Hughes, “The Basics of gravitational wave theory,” New J. Phys. 7 (2005) 204, arXiv:gr-qc/0501041 [gr-qc].
  * (106) P. Schwaller, “Gravitational Waves from a Dark Phase Transition,” Phys. Rev. Lett. 115 no. 18, (2015) 181101, arXiv:1504.07263 [hep-ph].
  * (107) C. Caprini et al., “Science with the space-based interferometer eLISA. II: Gravitational waves from cosmological phase transitions,” JCAP 1604 no. 04, (2016) 001, arXiv:1512.06239 [astro-ph.CO].
  * (108) J. Jaeckel, V. V. Khoze, and M. Spannowsky, “Hearing the signal of dark sectors with gravitational wave detectors,” Phys. Rev. D94 no. 10, (2016) 103519, arXiv:1602.03901 [hep-ph].
  * (109) A. Katz and A. Riotto, “Baryogenesis and Gravitational Waves from Runaway Bubble Collisions,” JCAP 1611 no. 11, (2016) 011, arXiv:1608.00583 [hep-ph].
  * (110) N. Bartolo et al., “Science with the space-based interferometer LISA. IV: Probing inflation with gravitational waves,” JCAP 1612 no. 12, (2016) 026, arXiv:1610.06481 [astro-ph.CO].
  * (111) D. J. Weir, “Gravitational waves from a first order electroweak phase transition: a brief review,” Phil. Trans. Roy. Soc. Lond. A376 no. 2114, (2018) 20170126, arXiv:1705.01783 [hep-ph].
  * (112) D. Croon, V. Sanz, and G. White, “Model Discrimination in Gravitational Wave spectra from Dark Phase Transitions,” JHEP 08 (2018) 203, arXiv:1806.02332 [hep-ph].
  * (113) J. Ellis, M. Lewicki, and J. M. No, “On the Maximal Strength of a First-Order Electroweak Phase Transition and its Gravitational Wave Signal,” JCAP 04 (2019) 003, arXiv:1809.08242 [hep-ph].
  * (114) M. Breitbach, J. Kopp, E. Madge, T. Opferkuch, and P. Schwaller, “Dark, Cold, and Noisy: Constraining Secluded Hidden Sectors with Gravitational Waves,” JCAP 1907 no. 07, (2019) 007, arXiv:1811.11175 [hep-ph].
  * (115) A. J. Helmboldt, J. Kubo, and S. van der Woude, “Observational prospects for gravitational waves from hidden or dark chiral phase transitions,” arXiv:1904.07891 [hep-ph].
  * (116) C. Caprini et al., “Detecting gravitational waves from cosmological phase transitions with LISA: an update,” JCAP 03 (2020) 024, arXiv:1910.13125 [astro-ph.CO].
  * (117) W.-C. Huang, M. Reichert, F. Sannino, and Z.-W. Wang, “Testing the Dark Confined Landscape: From Lattice to Gravitational Waves,” arXiv:2012.11614 [hep-ph].
  * (118) M. Markevitch, A. Gonzalez, D. Clowe, A. Vikhlinin, L. David, W. Forman, C. Jones, S. Murray, and W. Tucker, “Direct constraints on the dark matter self-interaction cross-section from the merging galaxy cluster 1E0657-56,” Astrophys. J. 606 (2004) 819–824, arXiv:astro-ph/0309303.
  * (119) J. L. Feng, M. Kaplinghat, H. Tu, and H.-B. Yu, “Hidden Charged Dark Matter,” JCAP 07 (2009) 004, arXiv:0905.3039 [hep-ph].
  * (120) M. R. Buckley and P. J. Fox, “Dark Matter Self-Interactions and Light Force Carriers,” Phys. Rev. D 81 (2010) 083522, arXiv:0911.3898 [hep-ph].
  * (121) M. Vogelsberger, J. Zavala, and A. Loeb, “Subhaloes in self-interacting galactic dark matter haloes,” Monthly Notices of the Royal Astronomical Society 423 no. 4, (07, 2012) 3740–3752. https://doi.org/10.1111/j.1365-2966.2012.21182.x.
  * (122) J. Zavala, M. Vogelsberger, and M. G. Walker, “Constraining self-interacting dark matter with the Milky Way’s dwarf spheroidals,” Monthly Notices of the Royal Astronomical Society: Letters 431 no. 1, (02, 2013) L20–L24. https://doi.org/10.1093/mnrasl/sls053.
  * (123) S. Tulin and H.-B. Yu, “Dark Matter Self-interactions and Small Scale Structure,” Phys. Rept. 730 (2018) 1–57, arXiv:1705.02358 [hep-ph].
  * (124) K. Bondarenko, A. Sokolenko, A. Boyarsky, A. Robertson, D. Harvey, and Y. Revaz, “From dwarf galaxies to galaxy clusters: Self-Interacting Dark Matter over 7 orders of magnitude in halo mass,” arXiv:2006.06623 [astro-ph.CO].
  * (125) R. K. Leane and J. Smirnov, “Exoplanets as Sub-GeV Dark Matter Detectors,” Phys. Rev. Lett. 126 no. 16, (2021) 161101, arXiv:2010.00015 [hep-ph].
  * (126) R. K. Leane and T. Linden, “First Analysis of Jupiter in Gamma Rays and a New Search for Dark Matter,” arXiv:2104.02068 [astro-ph.HE].
  * (127) G. Barenboim, N. Blinov, and A. Stebbins, “Smallest Remnants of Early Matter Domination,” arXiv:2107.10293 [astro-ph.CO].
  * (128) Y. Cui, M. Lewicki, D. E. Morrissey, and J. D. Wells, “Probing the pre-BBN universe with gravitational waves from cosmic strings,” JHEP 01 (2019) 081, arXiv:1808.08968 [hep-ph].
  * (129) P. Auclair et al., “Probing the gravitational wave background from cosmic strings with LISA,” JCAP 04 (2020) 034, arXiv:1909.00819 [astro-ph.CO].
  * (130) Y. Gouttenoire, G. Servant, and P. Simakachorn, “Beyond the Standard Models with Cosmic Strings,” JCAP 07 (2020) 032, arXiv:1912.02569 [hep-ph].
  * (131) A. Hook, G. Marques-Tavares, and D. Racco, “Causal gravitational waves as a probe of free streaming particles and the expansion of the Universe,” JHEP 02 (2021) 117, arXiv:2010.03568 [hep-ph].
  * (132) C.-F. Chang and Y. Cui, “Gravitational Waves from Global Cosmic Strings and Cosmic Archaeology,” arXiv:2106.09746 [hep-ph].
  * (133) F. Ertas, F. Kahlhoefer, and C. Tasillo, “Turn up the volume: Listening to phase transitions in hot dark sectors,” arXiv:2109.06208 [astro-ph.CO].
  * (134) T. A. DeGrand and K. Kajantie, “Supercooling, Entropy Production and Bubble Kinetics in the Quark - Hadron Phase Transition in the Early Universe,” Phys. Lett. B 147 (1984) 273–278.
  * (135) B. Lucini, M. Teper, and U. Wenger, “Properties of the deconfining phase transition in SU(N) gauge theories,” JHEP 02 (2005) 033, arXiv:hep-lat/0502003.
  * (136) L. Giusti and M. Pepe, “Equation of state of the su(3) yang–mills theory: A precise determination from a moving frame,” Physics Letters B 769 (Jun, 2017) 385–390, arXiv:1612.00265 [hep-lat]. http://dx.doi.org/10.1016/j.physletb.2017.04.001.