Thermal Squeezeout of Dark Matter

We carry out a detailed study of the confinement phase transition in a dark sector with a $SU(N)$ gauge group and a single generation of dark heavy quark. We focus on heavy enough quarks such that their abundance freezes out before the phase transition and the phase transition is of first-order. We find that during this phase transition the quarks are trapped inside contracting pockets of the deconfined phase and are compressed enough to interact at a significant rate, giving rise to a second stage of annihilation that can dramatically change the resulting dark matter abundance. As a result, the dark matter can be heavier than the often-quoted unitarity bound of $\sim100~$TeV. Our findings are almost completely independent of the details of the portal between the dark sector and the Standard Model. We comment briefly on possible signals of such a sector. Our main findings are summarized in a companion letter, while here we provide further details on different parts of the calculation.


Introduction
The cosmic abundance of dark matter (DM) is comparable to the abundance of Standard Model (SM) particles up to an O(1) factor [1]. The similarity of these two ostensibly unrelated abundances raises the suspicion that the two sectors may have been in chemical equilibrium at some point in their history, implying some sort of interaction portal between the SM and DM. Numerous experimental efforts to look for such a portal are under way. Nonetheless, the particle nature of DM and any potential portals to the SM remain unknown at present.
In addition to probing interactions between DM and the SM, many experimental and theoretical efforts aim to probe possible dynamics within the dark sector itself. Often simplified dark sectors with only a single DM particle are considered. Yet, the rich gauge structure of the SM offers no particular reason to believe that the dark sector will be significantly simpler. A wide range of more involved dark sectors have been studied, especially scenarios with a new confining force in the dark sector; see for instance Refs. . Depending on the details of the sector, different hadronic states can be the DM candidate in different theories and the stabilizing symmetry and the DM mass scale can vary widely [32]. Such a sector can further give rise to rich dynamics that can potentially solve other problems in the SM as well, e.g. see Refs. [38,6,39] where the observed baryon asymmetry is tied to the DM abundance.
An interesting class of confining dark sector models is the scenario where all the dark quarks are substantially heavier than the dark confinement scale, Λ. These models and their experimental signals have been studied extensively, see for instance Refs. [6,[26][27][28]30]. For sufficiently heavy quarks, lattice calculations have shown that the phase transition in such a sector is of first-order for SU (N ) with N = 3 [40][41][42][43][44] or N > 3 [45,46]. There has been a recent surge of interest in the study of the potential effects of first-order phase transitions in other DM models, e.g. see Refs. [47][48][49][50][51][52] 1 , but the effects of the phase transition on the relic abundance of dark matter in confining dark sector models are mostly unexplored, with the exception of a recent study of dark sectors with only light quarks (m q Λ) [56].
In this work, we consider the simplest such confining model -an SU (3) gauge theory with one heavy quark in the fundamental representation -and focus on the effects of the first-order phase transition on the DM relic abundance calculation. Similar to the arguments put forward in Ref. [57], we will argue that toward the end of the phase transition we will be left with pockets of the high temperature, i.e. deconfined, phase submerged in a sea of low temperature, i.e. confined, phase. We will argue that the heavy quarks are all initially trapped inside these contracting pockets. 2 To determine important properties of these pockets such as their initial characteristic size and contraction rate, we will develop a simplified model to numerically simulate the phase transition.
As a pocket contracts, the dark quarks within it are compressed, allowing them to recouple and go through a second stage of annihilation. We calculate the fraction of the quarks that survive this new annihilation stage and escape the pockets in the form of dark baryons. We refer to this process as "thermal squeezeout", as the dark quarks are squeezed within the pockets and eventually leak out, in contrast to the standard "thermal freezeout". We find a dramatic suppression in the final abundance thanks to this phenomenon, which points to much heavier dark matter parameter space than was previously thought. In particular, the fact that the local DM density is much larger than the globally averaged DM density during this second stage of annihilation invalidates the homogeneity assumption made in the unitarity bound argument of Ref. [58]. 3 As a result, this model allows the DM candidate to be heavier than the perturbative unitarity bound on Weakly-Interacting Massive Particles (WIMPs), despite being thermal.
Since this second stage of annihilation is controlled by the dynamics within the dark sector and not the interaction between the dark sector and SM, our results are largely independent of the portal to the SM. In fact, we do not constrain ourselves to any specific portal in this paper. The only assumptions we make about such a portal is that (i) it exists, (ii) it keeps the SM and the dark sector in thermal contact during the phase transition, and (iii) it respects the dark baryon number that stabilizes the dark baryons. These assumptions streamline our calculations significantly. However, it is worth considering the possibility of models in which we can relax one or more of these assumptions; we leave this for future work.
The current work is merely the start of a broader program of studying such models in more details. The phenomenology of all such models should be revisited in light of the dramatic change in the relic abundance calculation. Depending on the gauge group 2 This depends on the representation of the quarks under the dark confining gauge group. For instance, if the quarks were in the adjoint representation (similar to the model in Ref. [29]), they could combine with the surrounding gluons, form color-neutral hadrons, and move into the confined phase. Similarly, if in addition to the heavy quarks the spectrum were to include light quarks in the fundamental representation as well, the heavy quarks would not remain trapped within the pocket as they could make a color-neutral bound state with one of many light quarks surrounding them to escape the pocket. 3 See also Refs. [59,60] for more recent studies of the unitarity bound on thermal DM models.
under study, the quarks' representation, and the portal, different models (with vastly different phenomenology) can be constructed. Our study indicates a natural window of DM masses between 1-100 PeV for such a setup. While conventional searches may lose sensitivity for such a high DM mass, the stochastic gravitational wave background due to the first-order phase transition in this scenario can be detected by planned future facilities. (See Ref. [61] for the projected reach of such facilities.) While this signal depends on the UV parameters that control the thermodynamics of the phase transition, it does not depend on the nature of the portal.
The rest of this paper is organized as follows. In Sec. 2 we provide an overview of the cosmology of our dark sector. In Sec. 3 we write down and solve the Boltzmann equations that determine the relic abundance of the dark matter. In Sec. 4 we provide an overview of the possible phenomenological implications of our dark sector before concluding in Sec. 5. We also provide three appendices to supply more details. In App. A we provide more details on the thermodynamics of first-order phase transitions in an expanding universe and detail a simulation we performed to fix the phase transition parameters that enter the relic abundance calculation. In Apps. B and C we review some results in the literature for the cross sections and binding energies of heavy quarks and their bound states.

A Qualitative Overview of the Cosmology
We consider a dark sector with a non-abelian SU (3) gauge group and a single flavor of heavy quarks q in the fundamental representation where G µν is the dark gluon field strength and m q is the dark quark mass with m q Λ (in practice we consider m q 100Λ), where Λ is the dark confinement scale at which a phase transition takes place.
Given that m q Λ we expect that such heavy quarks decouple from the thermal bath well before the phase transition, so that the phase diagram of this model is very close to that of pure Yang-Mills for T m q . Since the heavy quark regime can be well-approximated by the pure-gauge regime, the phase transition behavior is almost independent of the number of heavy quark flavors [62]. The only constraining condition on the number of quark flavors is then asymptotic freedom, or in special cases asymptotic safety [63]. Various lattice gauge theory studies have established that the SU (3) phase transition takes place at a critical temperature very near the confinement scale, T c ≈ Λ, and is first-order [41][42][43][44]. We will therefore assume that this dark sector features a first-order phase transition exactly at T c = Λ.
The effects of this phase transition on the relic abundance of DM have been relatively unexplored in previous studies of such a confining gauge sector, e.g. [27]. We will discuss in this section how this phase transition dramatically changes the relic abundance calculation by causing a second stage of significant DM annihilation.
We remain agnostic about how the dark quark mass is generated as it will not affect our study. We also do not commit to any specific portal between the dark sector and the SM. We merely assume such a portal exists and keeps the dark sector in kinetic equilibrium with the SM. The portal enables the decay of the dark glueballs and mesons to the SM while respecting the dark baryon number symmetry, thus stabilizing the dark baryons. These baryons, which are three-quark bound states, are the DM candidate in this setup. We will also assume a symmetric initial condition, n q = nq.
In this section we provide an overview of the cosmology of such a sector, focusing primarily on the effect that the phase transition has on the DM relic abundance. The goal is to provide the reader with a broad picture of the various moving parts in this study, while leaving some of the more detailed calculations for later sections and App. A.

Pre-confinement epoch
For high enough temperatures, T > T c , the dark sector exists in a deconfined thermal state in which quarks move freely within a gluon bath. Naively this setup seems at odds with confinement, which requires that colored objects not propagate freely over distances greater than the confinement length, Λ −1 . However, qualitatively, these colored quarks can move freely because they are connected to a network of thermal gluons [64]. These gluons screen a quark's color charge so that the quark effectively behaves like a color neutral object, in a process analogous to Debye shielding in plasmas. More quantitatively, lattice simulations have shown that when T ≥ T c , the potential between two heavy quarks flattens when they are separated by a distance of roughly more than Λ −1 [41]. In other words, distant quarks in a gluon bath do not influence one another.
In the deconfined phase, the quark relic abundance calculation proceeds analogously to a standard WIMP relic abundance calculation. For large enough dark quark masses, m q 20Λ, which will be satisfied for all the parameter space we consider in our analysis, the dominant number changing process qq ↔ gg freezes out before confinement. The Boltzmann equation governing this freeze-out is simplẏ n q + 3H(T )n q = − σv n 2 q − n eq q 2 , (2.2) where H(T ) is the Hubble constant at temperature T , σv is the thermally averaged annihilation rate for qq ↔ gg, n q is the quark number density, and n eq q is its value in thermal equilibrium with a thermal bath of temperature T and with zero chemical potential. Since T Λ in this epoch, σv can be calculated perturbatively [65], where α(m q ) is the dark, strong coupling constant evaluated at the dark quark mass scale and ζ is a prefactor encapsulating plasma effects and non-relativistic enhancements (with numerical values presented in Fig. 14). We will find that the exact size of this cross section does not qualitatively change the main results of this paper. For further elaboration about this cross section and others, see App. B. 4 For the running coupling constant we use [66] α(m q ) = 12π where N f (µ) is the number of light flavors contributing to the beta function at mass scale µ. We set N c = 3 and N f (m q ) = 1.
In the left panel of Fig. 1 we show the resulting quark number density evolution for specific choices of quark mass and confinement scale. A generic obstacle in the study of strong sectors is the uncertainty in determining cross sections. To characterize this uncertainty, we vary the cross section within an order of magnitude around the central value in Eq. (2.3), which produces the green bands.
Importantly, we find that heavy quarks are well-separated just before the phase transition begins. To characterize their separation, we define the typical inter-quark spacing in units of the confinement length, where t n is the time at which extensive bubble nucleation starts, i.e. the onset of the phase transition, and n q (t n ) is the number density of the quarks at this time. This quantity measures, in units of Λ −1 , the typical distance between quarks at the onset of the phase transition. When ξ(t) is large, quarks are separated by much more than a confinement length. In the right panel of Fig. 1 we show ξ(t n ) as a function of m q and Λ. Indeed, quarks are generically further from each other than a confinement length just as the phase transition begins. Were quarks not so well separated, the details of the phase transition would have a less dramatic effect on the DM relic abundance calculation and we would be able to use the combinatorial method of Ref. [27].

Bubble dynamics
Once the universe cools down to the critical temperature, T c = Λ, a first-order phase transition begins. Phase conversion cannot occur right at the critical temperature as both phases have the same free energy, so the temperature of the deconfined phase initially cools slightly below T c . As the deconfined phase supercools further into a metastable state, bubbles of the confined phase begin to nucleate and expand at a nonnegligible rate.
As deconfined phase is converted to confined phase, latent heat is released. In contrast to weakly coupled phase transitions, there is no perturbative parameter suppressing the latent heat, meaning that phase conversion will serve as a significant heat source in the temperature evolution of the universe. As a result, the plasma heats back up to a temperature very close to T c quickly after bubble nucleation becomes efficient. Since the nucleation rate is exponentially sensitive to the degree of supercooling, (T − T c )/T c , subsequent nucleation of bubbles is completely suppressed.
For the phase transition to continue, at least some of the bubbles from the brief period of efficient nucleation must continue to grow. To determine the bubble growth rate, we borrow an argument from [57]. As bubbles grow, the local temperature at the bubble walls increases towards T c , diminishing the free energy difference between the two phases that drives the expansion. The expansion rate is then limited by the rate at which the wall can cool. The cooling rate is controlled by the temperature gradient between the wall and surrounding fluid -if there were no temperature difference, heat would not flow. Since the wall temperature cannot exceed T c without reversing direction, we assume that this temperature difference is bounded above by the small degree of supercooling (T c − T )/T c . By modeling the heat dynamics near the wall in App. A, we estimate that the wall velocity is also bounded above by the degree of supercooling, v w ≤ (T c − T )/T c . For simplicity, we assume that v w saturates this bound.
In Fig. 2 we plot the degree of supercooling as a function of time during the bubble expansion stage of the phase transition. This result comes from a simple simulation that we develop to track the nucleation and growth of bubbles during this epoch. Further details about this simulation can be found in App. A. The stages of the phase transition discussed above are visible in this plot. The universe initially supercools through Hubble expansion until bubble nucleation becomes efficient, leading to quick bubble growth and latent heat injection that reheats the universe. The heating and cooling rates then roughly balance one another, leaving the temperature at a value very near T c , which suppresses further bubble nucleation as explained above.
Eventually, half of the universe converts to the confined phase. At around this so-called percolation time, most bubbles are in contact with one another and start coalescing. Soon after we are left with isolated "pockets" of the deconfined phase submerged in a sea of the confined phase. To properly compute the spectrum of shapes and sizes of these pockets would require a full numerical 3D bubble simulation. Instead, to simplify our analysis, we will assume that soon after percolation there is a characteristic size of a typical pocket, that pockets can be approximated as spherical, and that the details of the spectrum of pocket shapes and sizes will give only sub-dominant corrections to our results.
To determine this characteristic pocket size, we first determine the characteristic size of bubbles just before percolation. Using our simulation from App. A, we find that at percolation the spectrum of bubble radii peaks strongly at where M pl = 2.4 × 10 18 GeV is the reduced Planck mass. We now borrow another argument from [57] to argue that for most values of Λ, these bubbles coalesce quickly until they reach a larger characteristic size, denoted as R 1 .
The central idea is that small bubbles coalesce and merge quickly into bigger bubbles, and that the time scale for two bubbles in contact to merge becomes longer as bubble sizes grow. Intuitively, the larger the coalescing bubbles, the more matter has to be moved via the bubbles' surface tension, which takes more time. Thus there is a special bubble size, R 1 , above which bubbles merge slower than the timescale over which the phase transition takes place. We find this critical size to be (2.7) Figure 3 shows that for Λ > ∼ 1 TeV the typical size of bubbles just before percolation (R 0 ) is always smaller than R 1 . Thus, we assume that, for this range of Λ, at percolation all bubbles quickly coalesce until they reach a size of R 1 . For smaller Λs we assume that all bubbles will have radius R 0 instead. We then make the simplifying assumption that the characteristic size of pockets just after percolation is the same as the characteristic size of bubbles just before percolation, i.e.  The typical radius of bubbles just before percolation (blue line) and the characteristic coalescence radius R 1 (orange line). For any Λ with R 0 ≤ R 1 , we assume that bubbles quickly coalesce and grow to radius R 1 at percolation. Right: The asymptotic velocity of the pocket wall during its contraction as a function of the confinement scale when quark pressure is ignored. In the more realistic case where internal quark pressure is allowed to resist the contraction of the pocket, we expect v w to be much smaller, and to not necessarily asymptote to a constant value at small radii.
where R i is the characteristic initial pocket radius after percolation. It is more complicated to determine the wall velocity of the contracting pocket. The main complication is that quarks are trapped within pockets, which we will show in the next section, and will generically slow down the wall. For now we will assume that we can neglect the effect of the enclosed quarks, which will lead to an overestimate of v w . In Sec. 2.4 we will revisit the effect these quarks have on v w .
In App. A we find that at radii much smaller than R i , the pocket contraction rate asymptotes to a constant value, which is shown in the right panel of Fig. 3 and can be fit by (2.9) It will turn out in Sec. 3 that the relic abundance of DM is set while R R i , so we can neglect the initial stages when v w varies and treat it as a constant. The pockets' radii therefore shrink as a function of time according to where t is the time after percolation.
To the best of our knowledge, the problem of characteristic bubble properties, e.g. the wall velocity and charactristic size at percolation, is not completely settled for first-order phase transitions even in weakly interacting theories (see [67] and the references therein for recent discussions on calculating the wall velocity). In our numerical calculations in Sec. 3, we characterize these uncertainties by varying both v w and R i within one order of magnitude of the results shown in Fig. 3.
We note here that since the quark temperature is fixed near T c throughout the phase transition, the typical quark velocity is v q ∼ Λ/m q . (2.11) For the range of parameters that we will be interested in, we find v q v w . This inequality will become important in the next section when we analyze the effects that the walls have on the quarks.
To summarize, the phase transition begins with an initial, complicated stage of bubble nucleation and growth until bubbles come into contact with one another. It then enters an even more complicated bubble coalescence stage. The space between bubbles is made of pockets of the deconfined phase with the same characteristic size, R i . These pockets become isolated and eventually spherical, then contract initially with a velocity that is determined by the local heat diffusion rate. The contraction rate gradually slows down due to the pressure of the enclosed quarks, and the pockets eventually vanish. The phase transition has completed at this point, and the universe can proceed with its standard expansion history.
Further details of this phase transition, as well as an overview of the relevant thermodynamics can be found in App. A. Our study of the phase transition's effect on the DM relic abundance is insensitive to many of the details of the phase transition; we merely need an expression for the characteristic initial radius of pockets and their wall velocity, which are respectively provided in Eqs. (2.8) and (2.9). We emphasize that this latter expression for v w , which neglects the effect of quark pressure, overestimates the wall velocity during the contraction phase .

Heavy quarks during the phase transition
During the entire process of bubble nucleation and expansion described in the previous section, bubble walls run into quarks and anti-quarks. In this section we study these encounters in detail, and argue that the walls are impermeable to quarks, but permeable to color-neutral bound states. While we will focus on the interaction between walls and quarks, our conclusions hold for anti-quarks as well.
When a wall encounters a quark, the quark can push against it and deform it locally. Whereas in electroweak-like phase transitions a particle is able to penetrate through the wall at the cost of only a finite mass difference [47,48], the energy cost for an isolated quark to enter the confined phase is unbounded [41], preventing it from traveling far into the confined phase. Therefore, a quark can enter a bubble only if it either forms a color-neutral bound state before it enters the bubble, or deforms the wall so that it remains immersed in the color-screening gluon bath (see Fig. 4).
There are two ways in which the quark could form a bound state. First,qq pairs could be spontaneously created, binding with the quark as it passes through the wall. We can imagine a scenario as in Fig. 4 in which the quark pushes into the bubble and is connected to a gluon string [64] starting from its initial point of contact with the wall. If the quark were light enough, at some point this stretched string could break into ā qq pair and theq could bind with the incoming quark to form a color-singlet bound state that enters the bubble (see [56] for an example in which this process is efficient). However, for a heavy quark, the string breaking rate is extremely suppressed; this rate can be approximated [68] using the Schwinger mechanism [69], The exponential of the square of the large ratio m q /Λ makes this string breaking timescale much larger than all other timescales, completely shutting off this process. The inefficiency of string breaking and the quark's inability to pass through bubble walls is the main feature distinguishing our model from those that involve light quarks. The second way a quark could form a bound state is by encountering an anti-quark or two quarks somewhere within the deconfined phase, binding, then escaping into the confined phase before the bound state dissociates. These processes are important and will be analyzed in the next section via Boltzmann equations.
If a quark has not managed to bind into a color singlet state by the time it reaches a bubble wall, it deforms the wall to avoid entering the confined phase. As the wall deforms, its surface area increases, which increases the energy of the system. The surface tension therefore creates a force that opposes this deformation. If we estimate this force to be of order Λ 2 on dimensional grounds, then we find that the timescale for the surface tension to restore the shape of the wall and reverse the quark's velocity is This rebound timescale is much shorter than the string breaking timescale. It is also orders of magnitude smaller than the phase transition timescale, which we find in App. A to be t PT ∼ 10 −2 H −1 ∼ 10 −2 M pl /Λ 2 (also, see Fig. 2). Finally, the pocket contraction Figure 4: A depiction of a quark interacting with a phase boundary. As bubbles of confined phase grow (the wall moves to the right in the figure), their walls run in to quarks that move with typical velocity, v q , which is much larger than the wall velocity, v w (left configuration). The quark locally deforms the wall (center configuration), introducing an opposing force via surface tension. One can think of the quark as connected to the deconfined phase through a gluon string [64]. Since string breaking is shut off, i.e. the quark-gluon configuration does not have sufficient energy to pull a heavyqq pair out of the thermal background, eventually the quark comes to a halt and then rebounds back into the deconfined phase with its initial speed (the right configuration).
timescale is of order t contract ∼ R/v w . Since R > Λ −1 and we find that v w 10 −3 , we have t contract > 10 3 /Λ. Since we only consider quark masses that satisfy m q /Λ ≤ 10 2 in this paper, we have t rebound t contract . Since this rebound timescale is the shortest timescale in the problem, quarks rebound off walls very quickly before any other process can take place. Therefore, the bubble walls act like very stiff surfaces that quickly reflect quarks that come into contact with them.
As these bubbles grow, the walls sweep quarks and anti-quarks into the ever-shrinking deconfined regions, increasing the quark density over time. Moreover, since v q v w , quarks that are swept in can quickly travel through the shrinking deconfined region and maintain homogeneity, meaning that n q is independent of position in the pocket throughout the phase transition. Eventually, these particles end up inside the isolated pockets formed toward the final stage of the phase transition (see Fig. 5).
We note that in models with additional light (though not massless) quarks, there would be no first-order phase transition [70], and even if such a transition did exist, quarks would easily pass through bubble walls and would most likely be unaffected by the phase transition.
Within any fixed volume of the universe, including the isolated pockets, the baryon number is a fluctuating random variable. Although the baryon number averaged over all Once the temperature drops slightly below T c = Λ, bubbles of confined phase begin nucleating everywhere. The nucleated bubbles start growing and push quarks (black dots) around. Top Middle: The bubbles grow to a point that a O(1) fraction of the universe has converted into the confined phase. At this point bubbles start coalescing and quickly grow larger. Top Right: As the bubbles keep growing and combining, eventually we are left with isolated pockets of the deconfined phase submerged in a sea of the confined phase. Bottom Left: A single isolated pocket with quarks trapped inside it. Each pocket contracts as the phase transition continues. Bottom Middle: The particles in the pocket are compressed and their interactions recouple. During this phase the particles can either annihilate or bind with one another. Non-singlet states can not enter the confined phase, but once they form color-neutral bound states (orange dots) they can escape into the confined phase. Bottom Right: In the end, the pockets vanish and only the fraction of quarks that ended up inside color-neutral baryons survive. These particles diffuse away from the original pocket's position due to their local overdensity.
pockets must be zero due to our symmetric initial condition, any given pocket is expected to have an overabundance or underabundance of quarks relative to anti-quarks, which we call the pocket asymmetry, η. We find that the initial total number of quarks in a pocket, N initial q , is large, so that by the central limit theorem the standard deviation of fluctuations above and below the mean is N initial q . Therefore, no matter how efficient qq annihilation processes are in these contracting pockets, on average, at least a fraction of the initial quarks (or antiquarks) in a pocket will survive. This observation will have important consequences for our relic abundance calculation in the next section. Once isolated pockets have formed and their asymmetries have been set, they will contract and compress quarks and anti-quarks until formerly frozen-out interactions turn back on. These interactions includeqq annihilation as well as q + q binding via the attractive anti-triplet channel [65]. As these diquarks build up their occupation number, they can eventually bind with quarks to form color-singlet baryons that can quickly fly out of the pocket. 5 These escaping stable baryons constitute the DM candidate of our model, while the rest of the particles eventually dump their energy into the SM sector through an unspecified portal interaction.
We define a survival factor as the fraction of quarks and antiquarks that escape the pocket within baryons and antibaryons, In the next section we write down the Boltzmann equations governing the quark dynamics within contracting pockets and calculate this survival factor. As remarked above, S is bounded below by the asymmetry of a given pocket, and the expectation value of this lower bound is (2.16) After these surviving baryons escape the pockets and after the phase transition eventually completes, these baryons continue to diffuse away until they re-establish homogeneity in the universe. If the asymmetry bound is not saturated, baryons and anti-baryons can continue to annihilate as they diffuse outside of the pocket. A more detailed study of this final annihilation stage requires integrating inhomogeneous Boltzmann equations, which we leave for future work. In summary, the dark matter undergoes a short squeeze, where collapsing bubbles during the phase transition induce a second stage of rapid annihilation that drastically depletes the universe's pre-existing stock of dark matter. This extra annihilation after freeze-out opens up parameter space that had previously been ruled out due to overproduction of DM. We will show that this effect allows for this thermal DM to be heavier than the conventional unitarity bound of ∼ 300TeV [58].

The quark pressure on the wall
Although we have considered the effect that the wall has on the trapped quarks, we have ignored the effect that the trapped quarks have on the wall. In this section, we will argue that the trapped quarks generically slow down the contraction rate of the wall. 6 Much like a piston, a pocket wall can contract only if it does work on the enclosed gas of heavy quarks. Since we have assumed that this gas is thermally coupled to the rest of the SM bath, which has a much larger heat capacity than that of the dark sector, the quarks contract at constant temperature. Using an ideal gas equation of state, we can write the pressure of this quark gas as p q = n q T . The work that the wall does on the gas when it contracts by an amount dR is therefore p q dV = 4πR 2 n q T dR.
The forces that are responsible for this change of pocket radius are the surface tension and net gluonic pressure, the latter of which is directed inward whenever T < T c . We will show in the next section that during the earliest stages of pocket contraction, quark interactions are inefficient and the total number of quarks in the pocket is initially conserved. As a result, when the pocket shrinks the work required to contract the pocket grows like R −1 . At the same time, the work that the surface tension and net gluonic pressure do when contracting the pocket by dR is proportional to the change in area and volume, respectively, and so shrink like R and R 2 . Altogether, as R contracts, the forces pushing out grow while the forces pushing in shrink. We therefore expect that v w decreases with decreasing R.
As this physics involves non-equilibrium, strong dynamics, we cannot reliably compute v w as a function of R. Instead, in the remainder of this section we will argue that the effect of quark pressure is to slow down v w by orders of magnitude relative to the upper bound of Eq. (2.9) we computed when we neglected the quark pressure. For more details relevant to the following discussion, we refer the reader to App. A.4.
When we simulate pocket contraction while keeping track of the quark density within a pocket (see next section), we find that there always comes a point when the quark pressure has grown to such an extent that, were we to suddenly include it, the quark pressure would exactly oppose all inward-pointing forces. This point of mechanical equilibrium is defined by σ dA + ( p) dV = 0, where σ is the surface tension, dA the change in surface area, dV the change in volume, and p the sum of pressures acting on the wall. (Inward-pointing pressures are defined to be positive while outward facing pressures are negative.) If we were to suddenly include the effects of quark pressure, the motion of the wall would suddenly become calculable, since the state of the wall would be determined by equilibrium physics. The pocket would suddenly slow down and proceed to adiabatically shrink while maintaining mechanical equilibrium. Number changing processes would deplete n q , diminishing the quarks' outward-pointing pressure, and the universe would supercool further, increasing the net gluonic inward-pointing pressure. We find that in this scenario, v w suddenly drops by orders of magnitude when mechanical equilibrium is achieved, and v w steadily decreases many more orders of magnitude as the pocket contracts.
The discontinuous drop in v w signals a breakdown of our assumption that quark pressure was negligible before mechanical equilibrium was achieved. This simulation merely demonstrates that it is inconsistent to neglect quark pressure, and that it can potentially slow down the pocket contraction rate by orders of magnitude. We therefore expect that a more realistic simulation that correctly includes the effects of quark pressure from the very beginning will lead to a v w that gradually decreases from our upper bound of Eq. (2.9), which eventually overestimates v w by orders of magnitude.
We will use the results of our pocket evolution simulations to calculate a few parameters that enter the Boltzmann Equations that govern the abundances of various bound states in the pocket. While the expression for the pocket radius, Eq. (2.8), is robust to the uncertainties introduced by quark pressure, we argued that the wall velocity v w is sensitive to this uncertainty. In the next section, we will study the evolution of the bound state abundances in the pocket in two extreme cases: (i) when the effect of quark pressure on v w is completely ignored, or (ii) when its effect dramatically reduces v w .

Boltzmann Equations During Compression
As described above, toward the end of the phase transition, the deconfined regions form isolated pockets that contain all of the dark quarks. In this section, we describe the dynamics of the dark quarks and their bound states as the contracting pockets compress  them. The Boltzmann equations that we solve keep track of the many processes by which quarks either ultimately annihilate into gluons or form baryons that escape the pockets and become dark matter. We will solve the Boltzmann equations for a typical pocket with initial characteristic radius R i and pocket asymmetry set to its root-mean-square value, N initial q . We assume that the S of this typical pocket is approximately equal to S averaged over the full distribution of initial pocket radii and pocket asymmetries. The total number of DM particles that survive until today will then equal the total number of DM particles entering the phase transition times S.

Ingredients of the Boltzmann equations
We begin by listing the degrees of freedom that we will include in our Boltzmann equations, which have been tabulated in Tab. 1. We have neglected a host of exotic hadronic bound states like tetra-and penta-quark states because we assume that they are unstable and promptly decay to the states listed in Tab. 1. We also do not consider excited states of any of the bound states. To simplify the notation, we label states by their quark number throughout the text (for example, a baryon state is a 3 state while an anti-diquark is a −2 state).
We also neglect the mesonsqq in our analysis. This can be justified by comparing their decay rate to the fastest annihilation rate that we will encounter (see App. B) where the last inequality is obtained because we have heavy quarks (m q /Λ 100) and the interquark spacing in units of Λ satisfies ξ(t) ≥ 1. Such a fast meson decay rate ensures that these states are kept in equilibrium so that their number density is negligibly small. We also have verified numerically that including the mesons in our Boltzmann equations below has a negligible effect on our results.
Let us now look into the Boltzmann equation for the particles in Tab. 1 as they are compressed by the contracting pockets. We start with the Liouville operators. For the colored particles, i.e. ±1 (quarks/anti-quarks) and ±2 (diquarks/anti-diquarks), we have where the second term captures the effect of pocket compression. Notice that we have not included the usual factor of +3Hn i for the dilution of space due to Hubble expansion. As argued above, t PT H −1 . Therefore, the Hubble dilution rate is negligible during the phase transition and can be ignored.
For the color-neutral particles, i.e. baryons and anti-baryons, the compression term will be absent. Unlike the colored particles, the baryons are not constrained by confinement to remain in the deconfined pockets. The baryons formed in the pocket can then be thought of as a gas created in a container without walls. The gas of baryons will thus escape with a rate governed by its internal pressure, or equivalently by the thermal velocities of the baryons. 7 Once the baryons escape the pocket they are no longer tracked by the Boltzmann equations. So we must include baryon escape as a sink term in our Boltzmann equations, which we do by modifying the Liouville operator, To derive this escape rate, consider a small time step, dt. In each time step the pocket radius contracts by v w dt. The typical baryon moves a distance of about v q dt, where we ignore the distinction between the baryon and quark velocities. We then overestimate the escape rate by an O(1) factor by assuming that all baryons at the edge of the bubble move radially outward, giving a total number of escaped baryons of Notice that the justification for why baryons in the pocket are homogeneously distributed is different than that of the quarks and diquarks. Gradients in the baryon density naturally arise as the baryons flow from their high density points of creation to the low density exterior of the pockets. However, a homogeneous component of baryons is constantly being produced within a pocket due to the binding of (homogeneously distributed) quarks and diquarks. We find that the rate of production is faster than the escape rate, so the baryon density in the pocket remains homogeneous to a good approximation.
Combining this with the rate of change for pockets volume gives the density loss rate due to baryon escape used in Eq. (3.3). It will be convenient to track the evolution of the total number of particles in a pocket as opposed to number densities. Define the pocket volume, Multiplying the number density of species i by the volume of a pocket then gives the total number of species i in the pocket, It will also be convenient to replace the time coordinate with R using Eq. (2.10). We can then rewrite the Liouville operators as where N ≡ dN/dR and we have usedṘ = −v w . Now that we have dealt with the Liouville operators we write down the collision operators. We will only be concerned with 2-to-2 processes since n-to-2 processes are Boltzmann suppressed while 2-to-n processes are suppressed by extra factors of α(m q ) and phase space factors. We denote each of these terms with the following notation, (a, b) → (α, β) = σv ab→αβ n a n b − n α n β n eq a n eq b n eq α n eq with a, b, α, β = 0, ±1, ±2, ±3, and f ab,αβ ≡ N eq a N eq b N eq α N eq β . For gluons we have n 0 = n (eq) 0 , i.e. the gluons are always in equilibrium.
Once we have identified all the important interactions to be included in our Boltzmann equations, we can write down the complete system of differential equations for N i (R). We supply these equations with the initial conditions, which were derived in Sec. 2. The initial pocket radius is R i while the initial quark number in the pocket, N 1 , is found by multiplying the number density result of the pre-confinement freeze-out calculation in Eq. (2.2) by 4π 3 R 3 i . We find that the initial conditions for N 2 and N 3 are irrelevant, as they quickly approach an equilibrium value independent of whatever values we initially choose (so long as N 2 , N 3 N 1 initially). All that is left is to write down these equations and solve them.

Complete set of Boltzmann equations
The complete set of Boltzmann equations is: (3.10) where (· , ·) → (· , ·) is defined in Eq. (3.9). The right-hand side consists of all interactions that are consistent with quark number conservation. We also make the approximation that While this equality is not strictly satisfied due to the pocket asymmetry, we are able to make it because only one of three scenarios can occur: either (i) the symmetric component is never depleted to the point that the asymmetry is important, (ii) it is completely depleted and the accidental asymmetric abundance is all that survives, or (iii) the symmetric and the asymmetric components are comparable and our answer is off by an O(1) factor. As argued before, this asymmetry introduces a lower bound on S (Eq. (2.16)). Despite Eq. (3.10) having numerous terms, solving these equations numerically is rather straightforward. For convenience we list the important parameters entering into these equations and their reference values in Tab. 2. We remind the reader that Eq. (2.9)  overestimates v w since it neglects the quark pressure's ability to oppose pocket contraction. As we will discuss below, we can bracket the effect that a slower v w would have on the final DM relic abundance quite robustly, see Sec. 3.5 for further details. We also reemphasize that we have used simple approximations for some of other quantitiesparticularly the bubble radius -and a rigorous determination of them is only possible through more extensive numerical calculations. In Fig. 6, we show the solution of the Boltzmann equations in Eq. (3.10) for a specific quark mass and confinement scale when we neglect quark pressure and use Eq. (2.9) for the pocket wall velocity. There are a number of important observations to be made about this figure. First, the fractions of diquarks and baryons are initially very small, justifying why we did not include them in our calculations prior to pocket formation. Next, as the pockets contract, the number of bound states initially grows while the number of free quarks decreases due to binding or annihilation to gluons. As the number of free quarks decreases, the annihilation or escape of bound states become more important than their production, so their occupation numbers reach a maximum and monotonically decrease from there. Finally, we see that each step in the chain of bound state formation (q + q + q → qq + q → qqq) results in a suppression, i.e. the total number of diquarks is suppressed compared to the total number of free quarks, while the total number of baryons is suppressed compared to the diquarks. We anticipate that had we started with a larger SU (N ) gauge group (N 4), the bound states with higher quark numbers would have been further suppressed and the final DM survival factor would be lower. We leave a more detailed analysis of this scenario to future work.
Finally, to calculate the survival factor, we simply integrate Eq. (3.4) to calculate the total number of baryons that escaped during the contraction of the pocket and normalize to the initial quark number in the pocket. Rewriting Eq. (3.4) in terms of N and R, we find where we have used Eqs. (2.10) and (3.6) to change variables, and the subscript in S symm.
is to indicate that this is the survival factor of the symmetric component of the dark quarks. The factor of 3 in the first equality accounts for the fact that three quarks exists �� ���������� ���� ����� �������� Figure 6: Top: Evolution of the fractional number of free quarks (red), diquarks (green), and baryons (black) inside the pocket, normalized to the initial quark number. We neglect the effect of quark pressure on v w in solving the Boltzmann equations for this plot, meaning that we use Eq. (2.9) for v w . As the phase transition proceeds, the pocket radius decreases. Initially, almost all the quarks are unbound. As the pocket contracts, more bound states are formed and fewer quarks are found as free particles. As their numbers increase, the various states' annihilation rates increase as well. At some point their production and annihilation rates are comparable and the number of bound states inside the pocket reaches its maximum. The accumulative surviving fraction assuming zero pocket asymmetry predicted by Eq. (3.12) is denoted by the dotted purple line. The asymptotic value of this line is equal to S symm. . We denote the asymmetry lower bound on S from Eq. (2.16) by the orange dashed line.
Bottom: The DM abundance evolution for this mass and confinement scale. The T > Λ region is similar to Fig. 1; the confinement takes place at T = Λ and gives rise to the abundance suppression predicted by Eq. (3.13).
within every one baryon that escapes. Note that in deriving this result we assumed no asymmetry exists in the pocket. Combining this result with the lower bound on S from the asymmetry component, In Fig. 6 we show S symm. and η rms as well. We find that for the chosen Λ and m q , while neglecting the effect of quark pressure on the wall velocity, S symm. η rms , i.e. the local pocket asymmetry is not saturated during the contraction, but S is within ∼ 1 order of magnitude of this asymmetry bound. In fact, we find this is true for all the points in the parameter space that we study. In the upcoming section we will describe how we can leverage this observation to bracket the range of parameter space that gives rise to the correct DM relic abundance.

Analytic approximation
While the Boltzmann equations in Eq. (3.10) can be solved numerically, the large number of terms involved can muddle one's intuition. In this section, we develop a simple analytic approximation for solving these equations and determining S symm. . From the full set of interactions included in Eq. (3.10), we identify and neglect all but the most relevant processes that provide a closed set of equations with an analytic, asymptotic solution that shows good qualitative agreement with the numerical treatment. The subset of processes that we include are the formation of diquarks, and the subsequent capture of quarks that lead to the formation of baryons. The reduced set of Boltzmann equations is then The analytic solution for this set of equations is obtained relying on several assumptions.
• The initial dark quark abundance N initial q is determined by the pre-confinement freezeout of the elementary constituents.
• As long as the annihilation rate and the baryon escape rate in the contracting pocket is slower than the pocket contraction rate v w /R, the total quark number is conserved. Once those rates are of the same order, the annihilation process "recouples", and the free quark abundance drops toward zero. The condition Γ ann ≈ N initial .
• The initial number of bound states, N X (X = 2, 3), is negligible. As the pocket contracts, bound states start forming. Thus, we can write N X ∼ R −n , with n > 0, which implies N X ∼ N X /R. Inserting this into the Boltzmann equation shows that there is a small parameter controlling the rate of change, which is proportional to Thus expanding in δ the leading order result is obtained by setting N X ≈ 0, which is the equilibrium condition before the recoupling due to pocket contraction.
Given the above assumptions, before the annihilation process recouples, we have the quark number conservation Applying the equilibrium condition before recoupling and neglecting the escape and annihilation terms for the bound states at that point gives and ∆E i denote the heat released during the above processes. The solution to the above algebraic set of equations provide the abundances of quarks and bound states in the contracting pocket before recoupling as a function of the pocket radius R. The total abundance of produced color-singlet baryons is given by the total baryon abundance N 3 evaluated at the recoupling radius R rec . Notice that e −∆E 1,2 /Tc ∼ e −α 2 mq/Tc e −∆E 3 /Tc ∼ e −mq/Tc , where α is evaluated at the bound state's Bohr radius. Thus, we identify a strong hierarchyf 1 ,f 2 f 3 .
As a result, a simple analytic expression for the baryon fraction that survives the phase transition relative to the initial quark abundance N initial q can be found. In the limit of inefficient bound state breaking reactionsf 1,2 V 1 it is Thus, assuming all the cross sections are of the same order of magnitude, we see that the baryon survival factor is of order one, if deeply bound states dominate the system. This is the case if the scale hierarchy m q Λ is taken to be extremely large. In the regime of efficient bound state breakingf 1,2 V 1, we find stronger DM abundance suppression. To simplify things even further, we assume the terms withf dominate and thatf 1 ∼f 2 . 8 With these assumptions, we find that at the recoupling radius we have Now if we assume in Eq. (3.12) the integral is dominated by the contribution around the recoupling point where the bound states total numbers peak, we find . (3.20) We can better understand from this equation the effects that various parameters have on the survival factor. Increasing the quark velocity v q enhances their escape rate (see Eq. (3.4)) thus increasing S symm. . We also see that by increasing σv 1(−1)→00 , the survival factor decreases, which was expected since by increasing this cross section quarks annihilate more against each other instead of binding in bound states. For shallower bound states the binding processes are less favored, thus we expect the survival factor should decrease. This is exactly what Eq. (3.20) suggests: for shallower bound states, the Boltzmann suppression inf 1 becomes less severe andf 1 increases, thus S symm. decreases.
The initial density of quarks in a pocket is determined via a pre-confinement, perturbative freezeout calculation. Yet, N initial q in Eq. (3.20) depends on the initial pocket radius too. Thus, through N initial q we find that S symm. ∼ R −3 i .
Finally, by decreasing v w , according to Eq. (3.15), the recoupling radius increases, which gives less time for the baryon abundance in the pocket to build up before the interactions become efficient again, see Fig. 6. A larger recoupling radius means a smaller peak value for the N 3 abundance, like the one seen at R Λ ∼ 10 5 in Fig. 6, which in turn decreases the survival factor S symm. . This behavior is exactly what we see in Eq. (3.20).

The effect of quark pressure and summary of assumptions
The v w scaling of Eq. (3.20) helps us better understand how our determination of S would change had we included the effect of quark pressure on v w . This equation suggests that by using Eq. (2.9) for v w and ignoring the fact that quark pressure can oppose pocket contraction, we are actually calculating an upper bound on the survival factor, since we are certainly overestimating v w . Also, this v w scaling, combined with the proximity of S symm. to the asymmetry bound η rms across our parameter space when we use Eq. (2.9), motivates us to believe that when the quark pressure is properly taken into account, we should expect that we saturate the the asymmetry bound S = η rms for every point in the parameter space that we study (see App. A.4 for more empirical evidence of this claim). In the upcoming section we use these two limits to bracket the parameter space of the model that reproduces the observed DM relic abundance. We refer to these two limiting scenarios as the zero quark pressure and the asymmetry scenarios.
Before jumping to the result of solving the Boltzmann equations, it is useful to review all the parameters affecting our calculation of S and the final DM abundance. The UV model has a very limited set of parameters: the confinement scale, Λ, and the quark mass, m q . These parameters feed into the calculation of a few secondary quantities that directly affect the calculation of S and are listed in Tab. 2. A precise calculation of these secondary quantities requires various non-perturbative studies. These quantities can be divided into two broad categories: macroscopic and microscopic.
The macroscopic quantities are those concerning the dynamics of the bubbles and pockets, i.e. their initial radius R i and their wall velocity v w . While our expressions for these quantities in Eqs. (2.8) and (2.9) were based on a simplified simulation of the phase transition (see App. A), there is extensive literature concerned with the detailed calculation of these quantities. Unfortunately, this literature has not yet settled on a single, definitive calculation of these quantities, which is why we content ourselves with simple order of magnitude estimates. (See, for example, Refs. [71][72][73]67] and references within for various calculations of the wall velocity.) The microscopic quantities include various cross sections and binding energies. They also determine the dimensionless inter-quark spacing, ξ, which directly affects our final results as well. We use the results from [65] for the cross sections and the binding energies. We summarize the relevant quantities in Apps. B-C. It is also worth reiterating a few important assumptions that significantly streamlined our analysis. Recall that in Sec. 2.2 we argued that the wall velocity is controlled by the amount of supercooling and quark pressure during the phase transition. Following that assumption, we found that the typical velocity of quarks in Eq. (2.11) is much faster than the wall velocity even when the quark pressure effect is neglected in Eq. (2.9). Therefore, any density gradient within a pocket caused by the compression of the walls can be quickly smoothed out by the thermal motions of the quarks. As a result, we assume that the particles within the pockets are homogeneously distributed, which simplifies our analysis significantly.
We also neglect the abundance of the bound states before the phase transition. Furthermore, as suggested in Fig. 1, we assume the quarks are initially well-separated inside the pockets, and that they rebound off the wall surface promptly. As we will argue later, all of our assumptions determine the parts of the parameter space where our analysis is valid.

Results and discussion
We now turn to the central results of this paper. We scan over a range of Λ and m q /Λ values, solving the Boltzmann equations at each point to calculate the survival factor S. As mentioned above, we use Eq. (2.9) for the wall velocity when solving the equations and finding the viable part of the model's parameter space that produces the correct present-day abundance of DM. We argued that this zero quark pressure scenario and the asymmetry scenario, in which we assume S = η rms , are the two limiting cases that bracket the uncertainties in our DM relic abundance calculation. We will find that these two scenarios only give rise to an O(1) difference in the DM mass range that can explain the observed relic abundance.
In Fig. 7 we show contours of constant survival factor for both these scenarios and for different values of Λ and m q /Λ. The asymmetry scenario plot shows the smallest survival factor possible while the zero quark pressure scenario gives an upper bound on the survival factor for every point in the parameter space, see the discussion in Sec. 3.3. In the asymmetry scenario the only sources of uncertainty are those affecting the preconfinement calculation and the initial pocket size, while in the zero quark pressure scenario the uncertainty in determining the wall velocity v w should also be included.
The available parameter space in the asymmetry limit scenario is shown in Fig. 8. Figure 7: Contours of constant survival factor S (green contours) in the two limiting scenarios that we consider: (i) assuming the asymmetry bound on S is saturated (on the left), or (ii) neglecting the quark pressure effect on the pocket wall velocity (on the right). The contours of constant DM mass in TeV are shown as well (black dashed). To obtain the right plot, we solve the full set of Boltzmann equations in Eq. (3.10) using the values for the initial pocket radius R i and the pocket wall velocity v w based on our simulation results discussed in Sec. 2 and App. A. In the light blue region on the right plot we find that the suppression is so severe that only the accidental asymmetric abundance of quarks in each pockets survives after the phase transition. For the left plot we simply assume S = η rms for every point in the parameter space. We observe orders of magnitude suppression in the DM abundance due to the second stage of annihilation during the phase transition in either scenarios. Note that the small difference in the 10 −7 contours in the overlap region is a plotting artifact. Equation (2.7) shows that as Λ increases, R 1 , and so the number of trapped quarks inside the pocket, decreases. Thus, as expected from Eq. (2.16), we find that the larger the initial radius, the smaller the survival factor.
We should keep in mind that many simplifying approximations were made about the dynamics of the phase transition in App. A in order to obtain Eq. (2.8) for the bubble radius. This, inevitably, introduces some uncertainty in our calculation. To characterize this uncertainty, in Fig. 8 we introduce a fudge factor for the bubble radius denoted by f R , to be multiplied against the values from Eq. (2.8). The observed relic abundance line moves within the light purple band as we vary f R between 0.1 and 10. Any point above and to the right of the relic abundance line, including the entire red region, is ruled out.
Since the asymmetry limit scenario was the lowest attainable S in our setup, the relic abundance line in this scenario is an upper bound on the possible masses in our model.
In the other limit, the zero quark pressure scenario provides us with a lower bound on the range of DM masses in this setup that can explain the observed DM abundance. In Fig. 9 we show the available parameter space in this scenario. The calculation can Figure 8: The produced abundance of dark baryons, the DM candidate, in the asymmetry abundance scenario. The black dashed lines are contours of constant DM mass in TeV. The relic abundance line (with the initial radius fixed to its central value, i.e. f R = 1) is plotted (purple line) along with its uncertainty (light purple shades) corresponding to an order of magnitude variation of the initial radius. The shaded red region is excluded, as it produces too much DM, while the unshaded region produces too little DM. The baryons can therefore constitute a sub-component of the DM within the unshaded regions of parameter space. The survival factor is determined by the accidental asymmetry of the pocket, which is independent of the wall velocity as long as the asymmetry bound is saturated. Thus, the major source of uncertainty in the location of the relic abundance line is the initial pocket radius. We also find that the uncertainty from microscopic quantities is sub-dominant to those of initial pocket radius. This figure clearly shows that the baryon masses accounting for the observed DM abundance can be much heavier than the unitarity bound [58]. now be affected by a change in both the initial pocket radius R 1 and its wall velocity v w . To characterize this uncertainty, in Fig. 9 we introduce a fudge factor for both the bubble radius and the wall velocity, denoted by f R and f v , respectively. As expected, for any fixed Λ the observed DM relic abundance in this scenario is obtained by smaller DM masses than that of the asymmetry scenario in Fig. 8.
In this figure we also show the relic abundance line when we use the analytic approximation of Eq. (3.20) to calculate the survival factor. We find reasonable agreement between our analytic approximation (the orange curve) and the full numerical result (the purple curve).
In both these limiting scenarios studied in Figs. 8-9, we find a similar range of DM masses and Λ that can account for the present-day DM abundance. We expect that these two scenarios bracket the true location of the relic abundance line when the effect of the quark pressure on the pocket wall velocity is appropriately included. The figures indicate to characterize the uncertainty in the final relic abundance calculation stemming from these quantities. The relic abundance line using the analytic approximation of Eq. (3.20) for the survival factor is denoted by the orange curve too. Even with the quark pressure neglected, we still find a substantial suppression in the DM abundance during the phase transition. We still find that the baryon masses accounting for the observed DM abundance can be much heavier than the unitarity bound [58]. that, depending on the macroscopic parameters, the region of parameter space that produces the observed DM abundance predicts m DM ∼ O(1) − O(100) PeV, well above the thermal relic unitarity bound of m DM 300 TeV [58]. Even with various sources of uncertainty, our results predict a confinement scale roughly in the O(1) − O(100) TeV range, in contrast to [27], which predicts a much wider range of confinement scales in such models. The parameter space above this range is ruled out, while the remaining parameter space is allowed, producing a sub-component of DM.
The values of the cross sections and the binding energies entering the Boltzmann equations can be found in appendices B and C, respectively. We find that the uncertainty in our results due to these microscopic quantities is sub-dominant to the uncertainty from the macroscopic bubble dynamics parameters discussed above. For further details about these parameters and the uncertainties in determining them see the aforementioned appendices and the references therein.
Determining the exact position of the relic abundance line requires more precise calculations of both macroscopic and microscopic quantities. Nonetheless, such calculations will not change our qualitative conclusions: that the phase transition gives rise to a new stage of annihilation that reduces the relic abundance by orders of magnitude and shifts the DM mass to well above the unitarity bound.
We can also understand the expected results for the parts of parameter space not plotted. For larger Λs than were plotted, Figs. 8-9 suggests that this model always overproduces DM and is ruled out. For smaller Λs than were plotted, our assumption that the pre-confinement abundances of bound states are negligible breaks down. Since so many baryons are produced before the start of the phase transition, the survival factor becomes comparable to 1. For low enough Λ, we should use the combinatoric calculation of the relic abundance described in [27]. Further investigation of this region is left for future works.
As we go to larger values of m q /Λ, our assumption that v q v w breaks down. In this case, local inhomogeneities appear in the distribution of particles in the pockets and the entire homogeneous system of equations in Eq. (3.10) must be modified. Furthermore, we find that, for higher m q /Λ than is shown in Figs. 8-9, the quark separations during the contraction epoch can become as low as ∼ 1/Λ (due to the small cross sections allowing for a greater degree of compression). In this case, the picture of well-separated quarks that rebound off the stiff bubble wall (before they run into other colored particles) must be modified. Non-perturbative effects become more relevant in this case. It is also possible that at such high densities quarks bind into more stable and massive dark nuclear states such as nuggets, see Ref. [74] for a study of dark quark nuggets in the light dark quark limit and Ref. [75] for similar macroscopic objects.
Finally, as we go to lower values of m q /Λ, eventually the first-order phase transition turns into a second order one and then a cross over, see e.g. [44]. In this regime there will be no bubble walls to compress quarks into a second stage of annihilation. Even for lower values of m q /Λ for which there still exists a first-order phase transition we run the risk of breaking our assumption that the string breaking rate is negligible, so that quarks and diquarks can escape from the pocket before significant annihilation takes place.
All in all, outside of the window shown in Figs. 8-9, either the parameter space is already ruled out, or at least one of the simplifying assumptions we made fails and our analysis becomes unreliable.

Extensions of our analysis
So far we have focused on a confining SU (3) gauge group with a single generation of heavy fermions in the fundamental representation. Nonetheless, it is conceptually straightforward to repeat our analysis for slightly different setups. In this section we comment on the differences that we expect would have arisen had we varied the number of colors, N c , or the quark representation.
Had we chosen a gauge group with a larger number of colors, SU (N c 4), we expect that we would have found a smaller S symm. since the stable DM candidate in such a theory (the analogue of the baryon) requires more constituent quarks to bind together in more steps. (Notice that in Fig. 6 as the quark number of a state increases, its abundance decreases within a pocket.) However, if even with N c = 3 we find S symm. η rms , we expect to saturate the asymmetry bound for larger gauge groups as well. The additional S symm. suppression would not change the final survival factor S.
The quark representation under the dark gauge group has a slightly more complicated effect on our results. For any quark representation, one must first identify the list of all possible bound states and then write down the Boltzmann equations with all possible interactions. As explained in Sec. 3.3, the binding energies of these bound states can also have a significant effect on the solutions of the Boltzmann equations. As an example, consider the case in which quarks are in the adjoint representation of the group. These quarks can bind with gluons to form color-neutral gluequarks [29]. Since the gluons can be found abundantly, we expect that quarks can easily pass through pocket walls by binding with a nearby gluon. Thus, pocket walls will not compress the quarks, and there will be no second stage of annihilation due to the phase transition.
Besides changing the model under consideration, our work would also benefit from improving our order of magnitude estimates and simplifying assumptions. Dedicated numerical simulations that more carefully model the bubble dynamics and non-perturbative physics could reduce the uncertainties in both macroscopic and microscopic quantities listed in Tab. 2, narrowing down the uncertainty on the relic abundance line in Figs. 8-9. Finally, the glueball lifetime can have profound effects on our results and should be studied in more details in a concrete model. If they live long enough to dominate the energy budget of the universe, after their decay they inject entropy into the bath, further diluting the dark baryon. This pushes the relic abundance line in Figs. 8-9 to even heavier dark baryon masses. The lifetime of the glueballs depends strongly on the dimension of the portal operator that enables their decay to SM particles [76][77][78][79].
When this operator is dim-8, the glueballs are long-lived and a dilution in the dark baryon abundance is expected. When the operator is dim-6 instead, the glueballs decay fast and do not affect the final baryon abundance as long as m q /Λ 10 4 . Models that contain such dim-6 decay operators can contain a Higgs portal coupling to the dark sector quarks, see Ref. [27] for explicit models. The dark baryon relic abundances in these models should be re-evaluated taking into account the new effects discussed in this work.

Potential Experimental Signals
Our study so far only relied on fairly general properties of a dark sector. We only assumed the dark sector under study is a confining SU (3) gauge theory with a single generation of heavy fermions; we also assumed a portal exists between the sectors that keeps them in kinetic equilibrium and allows the glueballs and mesons to decay to the SM. All the conclusions drawn in the previous sections were independent of further details of the portal and the origin of the heavy dark quark mass.
A detailed study of all the phenomenological signals of such a sector has to be carried out in a model-dependent way with a specified portal. As a result, here we merely list the signals and constraints that should be expected from this broad class of models.
• The main feature of our setup is a first-order phase transition in the early universe.
Such a phase transition can also give rise to a stochastic gravitational wave (GW) background that can be detected in a host of different future experiments, e.g. see Refs. [80][81][82][83] for a recent study of the GW signals of confining dark sectors. The characteristics of the resulting GW, such as the frequency and the strength, depend on a handful of thermodynamical parameters, see [84] for a brief review.
This GW signal is independent of the portal to the SM. A naive estimate 9 shows that different parts of our parameter space could potentially be probed in future experiments like DECIGO [86,87] and BBO [88]. Early universe phase transitions can also give rise to anisotropies in the GW spectrum, which can potentially be detected at future facilities (see e.g. [89]). Given the extremely high mass range of the DM candidates in our model, the GW signals could have the highest discovery potential in such sectors. We leave the further study of GW signals in this class of models for future work.
• The glueballs and the mesons are unstable due to the portal to the SM. Stringent bounds from BBN require that these relics have a short lifetime. See for instance [90][91][92] for recent studies. As a rule of thumb, one can avoid various constraints by assuming all these bound states decay before the BBN, i.e. their lifetime is τ 1s. This bound on the lifetime introduces a lower bound on the strength of the portal. This lower bound can vary substantially depending on the details of the portal. Our requirement that both sectors are in kinetic equilibrium also imposes a lower bound, though we expect the BBN bound to be more stringent.
• The portal to the SM introduces possible direct and indirect detection signals. However, the DM number density in the universe and in our galaxy is very suppressed due to this model's heavy DM mass. A naive estimation suggests that our model's indirect signal from DM annihilation within the Milky Way is severely suppressed and undetectable. The direct detection signal, however, depends on the details of the portal and should be studied model-dependently. We note that in this heavy mass range even very large DM-SM elastic cross sections are allowed, but within the reach of upcoming and ongoing experiments [93].
• A separate indirect detection signal comes from the observation that our composite DM model admits excited states. De-excitations from these excited states might lead to radiation that could be detected. Excitations could be produced in the early universe or via interactions with matter today.
• Yet another indirect signal could come from the capture of DM in celestial bodies see for example Refs. [94,95]. 10 As DM accumulates at the bottom of these potential wells, it can begin to annihilate at a significant rate, possibly affecting the evolution of these celestial bodies in an observable way or enhancing a potential annihilation signal [98,99].
• For the O(PeV) and above DM masses predicted in our model, direct production of DM at collider facilities is not possible in the foreseeable future. Yet, if the portal is substantially lighter, it can be directly observed at collider experiments. While the dark quarks are too heavy to produce at collider facilities, the glueballs of the new dark sector, whose mass is O(10Λ), could potentially be produced at future colliders.
• Various studies suggest an upper bound on the DM self-scattering [100][101][102][103][104][105][106]. As a rough estimate It is straightforward to check that for the high confinement scales we are studying, this upper bound is easily satisfied.
• One can also search for signals coming from the inhomogeneities in the DM density that were produced during the phase transition when DM was compressed by contracting pockets, but this seems unlikely. By performing a Jeans stability analysis we find that the internal baryon pressure easily overcomes the self-gravity of these overdensities. Pockets therefore do not seed self-gravitating DM clumps. One might also look for modifications to the matter power spectrum due to these overdensities, but initial estimates indicate that the matter power spectrum would only be modified at unobservably small mass scales if at all. Specifically, the total DM mass within a horizon radius soon after the phase-transition epoch (after which the comoving abundance is fixed) can be estimated as the DM density multiplied by H −3 , with a DM density crudely approximated (ignoring changes in the number of relativistic degrees of freedom over time) as ∼ (Λ/T CMB,0 ) 3 × the present cosmological density of DM, where T CMB,0 ∼ 2×10 −4 eV is the present-day temperature of the radiation bath. This gives an enclosed mass: Thus for phase transitions at the TeV scale and above, we would expect phasetransition-induced inhomogeneities to affect DM clumps at the kg scale and below. Even if these clumps survived, this mass scale is vastly lower than can be probed by any possible observational constraints on the matter power spectrum, which are currently exploring halo masses of order 10 7−8 M (e.g. [107][108][109]).
Because of its low number density, the dark matter in our setup can have significant interactions with the SM particles and still have escaped detection so far. Creative new search strategies will be needed to explore this possibility. Novel ideas for direct detection of such a scenario have been put forward in Ref. [93], and interesting signals in heavy isotope searches [110] could arise if our dark baryons can bind to SM atoms and nuclei.
In addition to the above signals, which should exist for any specific realization of the DM-SM portal, there may exist additional portal-dependent signatures. We also find, using the results of Ref. [27], that depending on the type of the portal to the SM the glueballs lifetime could be larger than the Hubble time at T = Λ. In such a scenario, the delayed decay of the glueballs can further dilute the DM abundance [29,32] in the parameter space that we have studied, thus pushing the relic abundance line in Figs. [8][9] to even higher DM masses. A proper study of this effect, as well as other signals from any specific portal, is left for future works.

Conclusion
In this work we studied the consequences of a first-order phase transition in a confining dark sector with a single heavy quark in the fundamental representation. We assumed a portal exists between our sector and the dark sector that keeps the two sectors at kinetic equilibrium at the time of the phase transition and respects dark baryon number conservation. The arguments we presented do not depend on further details of the portal.
We argued that the bubbles of the confined phase, after nucleation, expand very slowly. Soon after the bubbles come in contact and coalesce, pockets of the deconfined phase form, and are submerged in a sea of the confined phase. The quarks are trapped inside these isolated and ever-contracting deconfined phase pockets. There is always an accidental asymmetry in the net dark baryon number in a given pocket, due to local stochastic fluctuations in the number of quarks and anti-quarks at the onset of pocket formation. As the pockets contract, the enclosed quarks compress until formerly frozen-out interactions recouple, giving rise to a second stage of annihilation.
We wrote down the complete set of Boltzmann equations with all 2 → 2 interactions between lowest lying bound states of the heavy quark. By solving these equations, we were able to calculate the fraction of dark quarks that survive the second annihilation event. These surviving quarks bind into stable, color-singlet states that comprise the DM abundance we see today. We find that these Boltzmann equations predict a dramatic suppression in the DM relic abundance. This suppression is sensitive to the initial size of the pocket, the density of the quarks trapped within, and the pocket wall velocity. While there is a large uncertainty in determining these parameters, we showed that for virtually any reasonable values of these parameters there is a significant suppression in the DM relic abundance.
We find the effect of quark pressure on the pocket wall velocity difficult to model. However, we do know that this effect will further slow down the pocket wall, which we showed will imply a smaller survival factor. We calculated the relic abundance of DM in this setup in two extreme scenarios: (i) the zero quark pressure scenario, and (ii) when we assume the quark pressure is so severe that the asymmetry bound on the survival factor is saturated. These two limiting scenarios bracket the range over which the relic abundance line can move when the quark pressure effects are properly taken into account. We found that for a fixed dark confinement scale Λ, the DM mass in this setup only changes by O(1) factors between these two scenarios.
After identifying the parts of the m q − Λ parameter space that predict the observed present-day DM abundance, we found that this large suppression opens up parts of the parameter space that were previously thought to be ruled out. In particular, we found a DM mass scale well above the often-quoted unitarity bound. Our calculation also suggests an upper bound on the dark confinement scale Λ ∼ O(1) − O(100) TeV. For any Λ above this bound DM is overproduced, despite the dramatic suppression of its abundance during the phase transition. Depending on the value of Λ, the dark baryon mass that can explain the observed DM abundance varies roughly between 10 3 to 10 5 TeV.
There are many possible signals that our setup can give rise to. With the exception of gravitational waves, all the other potentially detectable signals depend on the specific form of this model's portal to the SM. It will be interesting to investigate the signatures of specific portals and their constraints, which we leave to future work.
There are numerous ways in which our analysis can be improved. To decrease the uncertainties in our results it will be important to perform more detailed numerical simulations of the macroscopic bubble dynamics during the phase transition and the microscopic strong dynamics that determine the particle interactions. The most important quantities to be calculated would be the initial pocket size and its subsequent contraction rate, and the cross sections and binding energies included in our Boltzmann equations. Additionally, it will also be interesting to study the relic abundance calculation for other gauge groups and different quark representations. Furthermore, for any specific portal we should study the potential DM dilution due to a delayed glueball decay after the phase transition.
Confining sectors are natural dark sector candidates. In this paper we focused on such sectors with only a single species of heavy dark quark. We have pointed out the dramatic effect that this model's first-order phase transition has on the relic DM abundance of such a sector. The dynamics lead us to a sharp prediction about the natural mass scale of such DM candidates, 10 3−5 TeV. It is of paramount importance to study variations of these minimal dark sectors in greater detail and their potential signatures in various upcoming experiments.

Acknowledgments
We thank Andrei Alexandru, Tom Cohen, Daniel Hackett, Julian Muñoz, Jessie Shelton, and Xiaojun Yao for useful discussions. We especially thank Benjamin Svetitsky for his collaboration in early stages of this work. We also thank Yonit Hochberg, Graham

A Thermodynamics of a First-order Phase Transition
In this Appendix we collect results that are helpful for understanding the dynamics of a first-order phase transition. We also describe a numerical simulation of the confining phase transition from the main text, which we perform to fix the initial pocket radius and its contraction rate. The code we used to perform this simulation can be found at https://github.com/gridgway/ConfiningPT Bubbles.

A.1 Standard thermodynamics
The defining characteristic of a first-order phase transition is that a first derivative of the free energy, or the free energy density f , is discontinuous at a critical point. In our case, the entropy density, s = − ∂f ∂T , is discontinuous, corresponding to a non-zero latent heat release due to phase conversion. In contrast, f itself is continuous, meaning that the free energy of the two phases are the same at the critical point, f deconf = f conf .
To better understand the latent heat, first notice that in the absence of chemical potentials the critical point is specified by a single parameter, the critical temperature T c , and the free energy density is minus the pressure, f = −p. Since the free energy of both phases are the same at the critical point, their pressures are the same. Using the Euler relation, ρ = T s − p, where ρ is the energy density, one can take the difference in energy densities of the two phases at the critical temperature to find ∆ρ = T c ∆s ≡ l , (A.1) which defines the latent heat density. Whereas the free energy density is continuous, the energy density is not. Another important quantity that characterizes first-order phase transitions is the surface tension, σ. The surface tension is the energy cost per unit area of an interface separating the two phases. The latent heat density and surface tension are calculable via lattice simulations. We assume that at temperatures below the quark mass, the thermodynamics of the system are insensitive to the heavy quark field and we need only consider lattice simulations of pure SU (3) Yang-Mills theory. In particular, the authors of [46] calculate a latent heat density and surface tension in the infinite volume limit (see their tables 7 and 15) These quantities have been computed elsewhere [111][112][113][114], and we find that the uncertainties on l and σ among these different lattice calculations are not large enough to qualitatively change our results. As the universe expands, our thermal system, still in the deconfined phase, supercools. 11 Once the system is supercooled, the free energy of the confined phase is lower than the free energy of the deconfined phase, meaning that it is energetically favorable for a phase transition to take place. For any non-zero amount of supercooling, there now exists a critical radius, R c , at which the energy cost of increasing a spherical bubble's surface area is exactly compensated by the free-energy decrease due to phase conversion, where ∆f is the confined phase free energy density minus that of the deconfined phase, and so is positive.
To relate ∆f to the latent heat density, first recall that at T c the entropy difference is given by ∆s = l/T c . If we then use the thermodynamic relation ∂f ∂T = −s we find ∂∆f ∂T Tc = −l/T c . If we assume small supercooling we can use a Taylor expansion, which at leading order gives We then find The total free energy of a bubble at the critical radius is Thermal fluctuations will randomly convert regions of varying shapes and sizes from one phase to the other. When the system is supercooled, there is a chance that a converted region will be large enough that it expands rather than contracts. The probability per unit time and volume of converting a region with free energy F is determined primarily by the Boltzmann factor e −F/T . Using dimensional analysis we have [117], where A is assumed to be some O(1) number that is roughly constant with respect to temperature and determined by the microscopic theory. Provided A is indeed an O(1) number, we find that the exponential is by far the more important factor for determining the behavior of the phase transition, so we set A to 1 without qualitatively changing our results. For an alternative though similar expression for Γ, See [116] §99. Let us make the simplifying assumption that all bubbles can be approximated as spherical. Then bubbles that nucleate with radii below R c are ephemeral, quickly shrinking due to surface tension, while bubbles with radii well above R c are exponentially less likely to nucleate than critical bubbles by Eq. (A.7). Then to a good approximation, we can assume that only bubbles at the critical radius nucleate. Combining Eqs. (A.6) and (A.7), we find the nucleation rate of these critical bubbles is where in the last line we have used the lattice results from Eq. (A.2). Because the latent heat density is an order one number in units of T c while the surface tension is a small number in units of T c , κ turns out to be a very small number (which Ref. [46] indicates is generically true for SU (N ) theories). As a result, very little supercooling is required before bubble nucleation becomes efficient.

A.2 First half of the phase transition: bubble growth
After bubble nucleation begins, the phase transition proceeds via the nucleation of new bubbles and the expansion of old bubbles. To keep track of the progress of the phase transition, we define the fraction of the universe that is in the confined phase [117], where t c is the time at which the universe first reaches the critical temperature and R(t, t ) is the radius at time t of a bubble nucleated at time t . Applying a time derivative yieldsẋ where the first term corresponds to the nucleation of new bubbles while the second corresponds to the expansion of old bubbles. The temperature evolution of the early universe plasma includes the usual adiabatic cooling term due to Hubble expansion, but now it also includes a new heating term due to the steady release of latent heat as the deconfined phase converts to the confined phase. By Eq. (A.1), converting a fraction dx from the deconfined to the confined phase releases dρ = ldx energy per unit volume. The released energy is absorbed in each phase with a temperature increase determined by each phase's respective specific heat, dρ/dT . For now, we will focus on the large scale average temperature so that we do not have to deal with small scale temperature gradients between points near and far from the sites of latent heat release. If we assume that the portal interaction between the SM bath and dark sector leads to frequent enough interactions between the two sectors per Hubble time, then in each phase the specific heat is dominated by the many degrees of freedom of the high temperature SM bath. For example, the deconfined phase of the dark sector is found to contribute about 5% to the specific heat at the critical temperature [114]. We can then use for both phases, where g * (T ) is the effective number of relativistic degrees of freedom and is g * (T ) ≈ 106.75 for all temperatures of interest in our analysis. We therefore approximate both phases as having the same specific heats and find The total temperature evolution of the universe during the phase transition is then given byṪ = −HT + 10 −2 T cẋ , (A.14) where the first term comes from the adiabatic cooling of relativistic species due to Hubble expansion. Had we considered a model in which few interactions take place between the standard model and dark sector baths per Hubble time, then the two sectors would be thermally decoupled. T would refer to the dark sector's temperature, which would heat relative to the SM temperature, and we would have divided by the dark sector's specific heat, eliminating the factor of 10 −2 in Eq. (A.14).
There is an important distinction between the temperature evolution of a weakly first-order phase transition and a strongly first-order phase transition. In a weakly first-order phase transition the latent heat in units of T c is typically small (see [117]) so that the heating term in Eq. (A.14) is negligible and the amount of supercooling does not change much due to the added latent heat. In a strongly first-order phase transition the latent heat in units of T c can be order one or larger (see Eq. (A.2)) so that T can be driven back up to T c before the phase transition completes. Furthermore, since κ in Eq. (A.9) scales inversely with l 2 , the large value of the latent heat decreases the amount of supercooling needed to achieve efficient bubble nucleation as is seen in Eq. (A.8), making it easier for the universe to reheat to a point where nucleation is negligible. We will show in a simulation below that this reheating scenario is achieved in our phase transition.
Notice also that T cannot ever reheat all the way up to T c . If it did, the critical radius would diverge and all bubbles would shrink. The second term in Eq. (A.14) would change sign and become a cooling term since latent heat would be absorbed, and hence T would be driven back below T c . Instead, |T − T c | decreases from its maximum to an equilibrium value very close to zero at which the heating and cooling terms in Eq. (A.14) nearly balance one another, as we will show below. This equilibrium phase coexistence is exactly the regime described by the Maxwell construction for first-order phase transitions (see [115], §84).
As explained in the previous section, we assume that bubbles nucleated at time t 0 are of size R c (T (t 0 )), giving the initial condition R(t = t 0 , t 0 ) = R c (T (t 0 )). To determine the radius at a future time, we need an expression for the bubble wall velocity, R(t, t ). An accurate treatment of the bubble wall velocity, requires full 3+1 dimensional numerical simulations of bubble dynamics during the phase transition. However, even in simplified settings, various numerical simulations have not converged on a single, definitive answer [71][72][73]67]. Instead, we will use a convenient, basic model ofṘ(t, t ). We require that critical bubbles not change their radius, and that larger bubbles expand while smaller bubbles contract. To capture this behavior, we use the simple functional formṘ where we define sign(0) = 0.
To determine v w (t) would require a better understanding of the underlying strong dynamics. Instead, we will estimate an upper bound on v w (t) based on thermodynamic arguments. As argued in the main text, the larger v w is, the less the DM relic abundance will be suppressed. So we will always set v w to its upper bound, assuming it to be a conservative choice.
As a bubble expands, it releases latent heat near its wall, locally heating the plasma at the wall to a temperature T wall > T and reducing the free energy difference at the interface. Since the free energy density is minus the pressure, this local heating reduces the net pressure acting on the wall (See Eq. (A.4)). Since the degree of supercooling is so small, the wall could potentially heat up to a temperature at which the net pressure balances against the surface tension. If the wall reached this temperature, bubble growth would no longer be thermodynamically favorable, and so the wall motion would come to a halt. By an argument completely analogous to the one that lead to Eq. (A.3), except we evaluate Eq. (A.4) at the wall temperature, we find this equilibrium wall temperature to be 16) We assume that as the wall temperature approaches T eq wall its growth slows down gradually. Before the wall temperature reaches T eq wall it will have slowed down to a steady state at which the rates of wall heating and cooling cancel one another. By estimating the rates of wall heating and cooling then setting them equal, we will determine an approximate expression for v w (t).
We start with the cooling rate. We will assume that T wall is very close to T eq wall , which is in turn very close to T c since the second term in Eq. (A.16) is very small compared to 1 for bubbles larger than Λ −1 . This assumption will lead to a faster v w . The fractional temperature difference between points near and far from the wall is then (T c − T )/T c . We assume that the heat loss rate is given by a diffusion equation,Ṫ cool ∼ −K∇ 2 T , and that the transport coefficient at T c is of order K ∼ Λ −1 . If we further assume that the length scale of the density gradient is Λ −1 , then we findṪ cool ∼ −Λ 2 (T c − T )/T c . Now we move on to the heating rate. We start with the energy injected per unit wall area and time, l v w . If we assume that this energy is injected within a typical length Λ −1 of the wall, then the energy injected per unit volume is lv w Λ. As before, dividing by the specific heat, dρ/dT , converts the energy increase into a temperature increase. Rather than assume that the specific heat is dominated by the SM degrees of freedom as we did when deriving Eq. (A.13), we assume that for this process the SM degrees of freedom are irrelevant and the specific heat is dominated by the dark sector degrees of freedom, though this choice will not affect our final result. We do so because we anticipate that for most models the portal interaction between the SM and dark sector will generically take place on timescales much slower than Λ −1 , and so will be inefficient compared to the interactions within the dark sector. Hence, we assume that the dark degrees of freedom disperse latent heat on a fast timescale of order Λ −1 , and only on much longer timescales does this heat find its way into the SM degrees of freedom. 12 Then we can use dρ DS /dT ∼ T 3 [114], which gives a heating rate ofṪ heat ∼ Λ 2 v w . 13 The wall velocity at which both rates are in balance is then As argued in Ref. [57], we see that v w is suppressed by the small degree of supercooling during the phase transition.
We apply a similar argument to the wall velocity of contracting bubbles, which are smaller than the critical radius. As bubbles contract, latent heat is absorbed near the wall, decreasing the wall temperature, increasing the net pressure on the wall, and thus opposing further contraction. The contraction rate is therefore limited by the rate at which the cold wall can heat due to heat flow from the hotter surrounding plasma. Balancing the heating and cooling rates as before leads to a wall velocity given by Eq. (A.17) up to a relative sign, justifying the symmetric functional form of Eq. (A.15).
Again, the above expression for v w is in reality an upper bound. An expanding bubble wall cannot move any faster than Eq. (A.17), at least for an extended period, because then it would locally overheat the wall. There is a similar consideration that the global rate of increase in temperature from bubble expansion should not significantly outpace the Hubble cooling, lest the whole universe be heated above T = T c . Consequently, we will see that the value of v w derived in Eq. (A.17) evolves so that it always lies below or close to this "global threshold".
The wall velocity (Eq. (A.17)) and global threshold are plotted in Fig. 10. At the start of the phase transition, when bubbles are rare and their radii are small, the global heating from bubble expansion is very small even for v w → c, and so the global threshold 12 We assume that the timescale over which interactions between the SM and dark sector baths exchange energy is much faster than the Hubble rate, justifying why the SM degrees of freedom are included in the heat capacity in Eq. (A.13). 13 If the portal interaction is actually efficient enough to keep the two sectors in equilibrium on this short timescale, then the heating rate would be suppressed by a factor of g . However, the cooling rate would also be suppressed by a factor of g since K is inversely proportional to the number density of interacting degrees of freedom, [118] §3.9. These two factors would then cancel in the expression for v w . We therefore expect Eq. (A.17) to be independent of the specific portal interaction. The wall velocity (blue) and "global threshold" (black dashed) above which phase conversion is so fast that the universe experiences net heating. The right panel displays the same data as the left panel, but zooms in to the very end of the phase transition when pockets have contracted significantly. The discontinuity in the middle is an artifact of our modeling. It occurs when the spectrum of bubbles, which peaks at R 0 , discontinuously jumps to a delta function spectrum of pockets centered at R 1 .
velocity goes to infinity. The same is true at the end of the phase transition, when there is very little volume available for phase conversion. However, during the phase transition when the rate of phase conversion is rapid, this threshold becomes relevant.
Let us consider what happens when bubble coalescence and growth causes v w = (T c − T )/T c to exceed this global threshold (which depends on the bubble density and typical bubble radius) for the first time. If v w overshoots the threshold, then the universe will begin to heat up on average, thus reducing v w . The net effect is for the degree of supercooling to evolve such that v w tracks the global threshold velocity. In practice, we observe that there is at first an abrupt drop in both v w and the global threshold velocity, associated with a sharp increase in the bubble number density due to nucleation; in this epoch, v w slightly exceeds the global threshold velocity and this supports a fairly rapid increase in T . Once nucleation becomes inefficient due to the rising temperature, the global threshold velocity evolves more slowly, driven by the expansion of the largest existing bubbles. In this epoch v w tracks the global threshold closely andṪ ≈ 0, with a slow adiabatic increase in temperature toward T c driven by the slow decrease in the global threshold velocity (which requires a corresponding decrease in v w and hence in the degree of supercooling). We can even derive the scaling of the global threshold (and v w ) with t during this period. By Eq. (A.14),Ṫ ≈ 0 impliesẋ ≈ 100H. Using the fact that the spectrum of bubbles is strongly peaked at a single radius R, and the nucleation of bubbles is so suppressed that the number density of bubbles is constant, we have x = 4πR 3 n bub /3. Combining the two equations gives dR/dt ∝ R −2 so that R = 3A(t − t 0 ) 1/3 for some constants A and t 0 . Therefore, v w = dR/dt = A(t − t 0 ) −2/3 . Indeed, we find that the slope of the line in Fig. 10 between times .001/H and .005/H is precisely −2/3.
The Hubble cooling can be relevant here, even though the phase transition takes place on timescales much smaller than a Hubble time, because the bubble expansion is so sensitive to the degree of supercooling; in contrast, we can freely drop e.g. density dilution terms corresponding to the Hubble expansion, as there is no comparably small density difference relevant to our calculation (see e.g. Eq. (3.2)). After coalescence, when the heating comes from pockets of shrinking radius rather than bubbles of expanding radius, the reverse process occurs, with a slow adiabatic decrease in the equilibrium temperature due to the decreasing size of the pockets. The temperature evolution eventually switches over to the standard Hubble cooling once v w = (T c − T )/T c can no longer reach the global velocity threshold (and henceṪ ≈ 0 cannot be maintained).
With Eqs. (A.11) and (A.14) through (A.17) in hand we are able to simulate the first half of the phase transition. This system of equations models the initial bubble nucleation and accompanying latent heat release, and the proceeding equilibrium regime up until percolation when bubbles begin to overlap and new dynamics must be included. Our initial conditions are that the universe has supercooled a little bit, T −Tc Tc = −10 −4 , the universe is fully in the deconfined phase, x = 0, and that no bubbles have nucleated yet. We evolve forward in small time steps, ∆t = 10 −6 /H(T c ). In each time step we nucleate Γ(T (t))(1 − x(t))∆t bubbles per unit volume at radius R c (T (t)) and add them to a list. We allow all other bubbles from previous steps to expand or contract by an amountṘ(t, t )∆t, which only depends on the bubble size and temperature within that time step. This procedure produces bubbles with radii less than or equal to 0, so we set such bubbles' radii to zero and remove them from our list. Additionally, many time steps result in an additional number density of bubbles that is so exponentially small that the computer sets the number density to zero. We remove these bubbles from our list, too. Each bubble nucleation and all bubble expansions increase x byẋ∆t and the temperature byṪ ∆t according to Eqs. (A.11) and (A.14). We finish our evolution once x = 1 2 . Outputs of our simulation are shown in Fig. 11. The left plot shows the degree of supercooling before percolation and the middle plot shows the fraction of phase converted, x. These plots make clear that the first half of the phase transition can be divided into three distinct stages. In the first stage, the degree of supercooling is so small that the bubble nucleation rate is too suppressed to have a significant effect on the simulation. During this period x = 0 and the universe cools through Hubble expansion. In the second stage, the supercooling reaches a point at which nucleation becomes efficient. These nucleated bubbles quickly grow and inject heat, corresponding to the sudden jump in the temperature and x a little before t = 0.001/H. The temperature reaches a point very close to T c at which nucleation of new bubbles becomes inefficient again, leading to the third stage. This stage is exactly the equilibrium phase coexistence regime described above. There are a fixed number density of large bubbles that grow and inject latent heat at such a rate as to cancel Hubble cooling. The net effect is thatṪ ≈ 0. By Eq. (A.14), we haveẋ Since the temperature is constant during this stage, the Hubble rate is as well, meaning that x grows linearly in time, which can be seen in our middle plot. This equation explains why the phase transition occurs over a time scale of 10 −2 /H. Had we assumed instead that the portal interaction between the SM bath and dark sector was very weak and led to few scatters per Hubble time, then the dark sector would not have access to the SM heat capacity and the above factor of 10 2 would be replaced by an O(1) factor instead. The rightmost plot in Fig. 11 shows the spectrum of bubble radii at the time of percolation, x = 1/2. The shape of this spectrum is a product of the preceding stages. The earliest bubbles were nucleated during the first stage. They had the longest time to grow, but were nucleated at a time of relatively small supercooling, meaning that their number density was exponentially suppressed. So as R increases away to the right of the peak of the spectrum, dn bub /dR decreases. Just as the second stage begins the universe is maximally supercooled. The bubbles produced at this point are the most numerous and constitute the peak of the distribution. In the rest of the second stage, the supercooling quickly diminishes, producing an exponentially suppressed population of bubbles that are smaller than the peak radius since they have less time to grow, leading to the sharp decrease of the spectrum to the left of the peak.
We define the peak in dn bub /dR to occur at R 0 and empirically find that it is very well fit by the function R 0 , however, is not the only relevant length scale for the bubbles. An additional length scale, R 1 , emerges from the dynamics of bubble coalescence. At percolation, bubbles frequently come into contact with one another and begin to coalesce. To model the coalescence dynamics, we borrow another argument from [57]. When two spherical bubbles of radius R coalesce into one larger-radius bubble, they decrease their surface area and therefore reach a more energetically favorable configuration due to surface tension. The energy difference between the two configurations is ∆E ∼ 4πR 2 (2 − 2 c . This change in energy is achieved by applying a force to the mass M in the bubbles over a characteristic distance R. So F ∼ ∆E/R ∼ M a, where a is the acceleration of the material in the bubbles. This acceleration can also be estimated as a ∼ R/t 2 coalesce , where t coalesce is the coalescence timescale. If we then use that the total mass in the two bubbles is M = 2 × 4 3 πR 3 ρ deconf ≈ 8 3 πR 3 T 4 c , where we used that ρ deconf (T c ) ≈ T 4 c [114], then we find The above equation shows that small bubbles coalesce quicker than large ones. Therefore, at percolation small bubbles will quickly coalesce until they reach a size R 1 past which t coalesce takes longer than the timescale of percolation. From Fig. 11 we can estimate the timescale of percolation as the time it takes x to change by a couple percent at around x = 1 2 . We find t perc ∼ 10 GeV is the reduced Planck mass and G is Newton's constant. Setting the percolation timescale equal to t coalesce then yields the critical bubble size of, in units of T c = Λ. In Fig. 3 we plot R 0 and R 1 as a function of the confinement scale. We find that for Λ 1 TeV , R 0 is less than or equal to R 1 . Therefore, for this range of Λ, once our simulation finishes at x = 1 2 we assume that all bubbles begin coalescing and quickly grow to radius R 1 . Since x does not change during this process, we assume that T remains fixed, too. For Λ < 1 TeV, we assume that coalescence is inefficient and bubbles remain at radius R 0 .

A.3 Second half of the phase transition: pocket contraction
After percolation, most bubbles are in contact with one another. The confined regions form a web and the deconfined regions form isolated pockets. We assume these pockets quickly attain spherical symmetry due to surface tension, and also that the typical size of a bubble before percolation is equal to the typical size of a pocket after percolation, R 1 .
Since all pockets are at the same initial radius, we can solve for the initial density of pockets. The number density of pockets that are all of radius R satisfies Since x = 1 2 at percolation, we have n pocket = 3 8πR 3 1 . (A.24) We will find that the degree of supercooling continues to be so small during pocket contraction that the nucleation of more bubbles is completely suppressed. Therefore n pocket remains constant and we find As before the contraction rate of the pockets is limited due to the latent heat release near the pocket wall. However, the wall velocity estimate for pockets is slightly different than it was for bubbles. Whereas supercooling results in a net pressure outward for bubbles, supercooling results in a net pressure inward for pockets, since the two phases are on opposite sides of the wall in either case. With this caveat in mind we can repeat our argument leading to Eq. (A.17).
As before, we argue that as it expands, the wall quickly heats up, approaching a threshold temperature at which pressure and surface tension are in equilibrium. This threshold now corresponds to slight superheating at the temperature T c 1 + 2σ l R . Before T wall reaches this threshold it achieves a steady state at which heating from latent heat injection cancels against cooling from heat diffusion from the wall. Using the same arguments as before, we find that this steady state corresponds to a velocity of Again, we used that 2σ lR = .03 R Λ 1, so we approximate the temperature difference between the wall and its surroundings appearing in the above equation as T c − T , which leads to an overestimate of the wall velocity.
There is another much more important new effect modifying the contraction rate of a pocket: quark pressure. The density of quarks trapped within pockets increases as they are forced within ever-shrinking volumes. This pressure opposes contraction, slowing down the wall velocity relative to Eq. (A. 26). For now, to build intuition, we will ignore the effect of quark pressure, and we will consider it in the next section.
With Eqs. (A.14), (A.25), and (A.26) we have all we need to simulate the second half of the phase transition. We use the same method and parameters as we did for the first half of the phase transition. In Fig. 12 we show the results of our simulation. The first instant of the simulation features a discontinuity in the temperature (and thus v w ). This discontinuity results from our discontinuous change of the spectrum of bubbles from the smooth form of dn/dR in the right panel of Fig. 11 to a delta function at R 1 . As the pockets get smaller, their heating rate diminishes. Hubble cooling becomes relatively more important over time leading to the steady decrease in T over time. By the very end of the phase transition, the temperature evolution is purely determined by Hubble cooling.
The right panel shows the pocket wall velocity as a function of the pocket radius. As a function of R, v w asymptotes to a well-defined value near 2 × 10 −4 for the choice of Λ shown in the figure. Naively, this plot seems to be at odds with the left panel, since the two plots are equivalent up to a minus sign. Whereas the velocity seems to asymptote to a constant value at late times, the degree of supercooling seems to vary greatly at late times. The apparent discrepancy is a result of the different x-axis scales. Whereas the x-axis in the left panel is linearly scaled in t, the x-axis in the right panel is log scaled in R. The majority of the simulation takes place when the pocket radii are very large. The pocket radii are much smaller than their initial value only for a very short time at the end of the simulation, at time scales much shorter than the 1/H. In fact, the small timestep 10 −6 /H was chosen to resolve this smaller timescale. At such small timescales and pocket radii, very little Hubble cooling or latent heat injection takes place, leading to the plateau in v w . Since most of the quark interactions recouple at pocket radii much smaller than R 1 (see e.g. Fig. 6), we are justified in treating v w as constant within the Boltzmann equations of Sec. 3. 14 The asymptotic v w value as a function of Λ is plotted in Fig. 3

A.4 The effect of quark pressure
Up to this point we have neglected the effect of quark pressure on the phase transition, For the first half of the phase transition, this approximate treatment was justified. At the start of the phase transition, p q is suppressed compared to the gluon pressure. During the phase transition n q grows, and so the quark pressure could potentially become large enough to affect the bubble dynamics. During the first half of the phase transition, however, the quark density, and hence p q , grows by only a factor of two. Including this factor of two enhancement, we find that the quark pressure is sub-dominant compared to the net gluon pressure during the bubble growth stage of the phase transition, and hence can be neglected. On the other hand, during the second half of the phase transition, the quarks are compressed much more. We find that for every point in the DM parameter space we consider, the quark pressure eventually becomes comparable to the other forces acting on the wall. Most likely, the increased quark pressure would oppose further contraction and slow v w down. Unfortunately, the process by which p q gradually grows and, in response, v w gradually shrinks, is a non-equilibrium, strong physics problem for which we have no solution. Nevertheless, we can still use thermodynamic arguments as before to understand the possible limiting behavior of v w .
Consider the scenario in which the enhanced quark pressure forces the pocket into a state of mechanical equilibrium. Mechanical equilibrium is achieved when the four forces on the wall -gluon pressure inside the pocket, glueball pressure outside the pocket, surface tension, and quark pressure -are in balance. At this equilibrium point, we must have where ∆f still refers to the confined phase minus the deconfined phase gluonic pressures. The temperatures are all evaluated at the local wall temperature and we have used N q = 4π 3 R 3 n q . If the system ever reached this equilibrium point, it would be a stable equilibrium: if R were to contract the increased quark density would oppose it; if R were to expand the surface tension would increase, the quark pressure would decrease, and the wall would absorb latent heat and increase the net gluon pressure difference, all of which oppose further expansion. After achieving mechanical equilibrium, the wall would proceed to contract adiabatically as quark annihilation decreases the quark pressure and wall cooling increases the net gluon pressure compressing the wall. By differentiating Eq. (A.29) with respect to time and defining v w = −Ṙ, we find Since we do not keep track of T wall , nor can we calculate v w as a function of T wall when the wall is out of mechanical equilibrium, we have no way of knowing when or if mechanical equilibrium is achieved.
However, to develop some intuition for the possible effects of quark pressure, we can perform a crude approximate calculation. For this calculation, we perform our pocket contraction simulation while simultaneously keeping track of the quark abundance within  Figure 13: The pocket wall velocity (left) and particle abundances (right) within the pocket using a model that crudely incorporates the effects of quark pressure on v w . We choose a confinement scale of 1 TeV and dark matter mass of 10 3 TeV. We begin the simulation neglecting the quark pressure, allowing us to use Eq. (A.26) to determine the contraction rate (red-dotted line). Eventually the quark pressure is so large that it is able to come into mechanical equilibrium with the other forces acting on the wall, which happens at the discontinuity near R Λ = 10 5 . At this point, we switch over to a contraction rate given by Eq. (A.30), leading to the discontinuous drop in v w in the left panel and the sudden depletion of all particle abundances in the right panel.

Evolution of Pocket Abundances
each step using the Boltzmann equations of Sec. 3. Since the quark pressure is subdominant initially, we start with v w given by Eq. (A.26). Eventually there comes a radius when the quark pressure is so large that it is able to oppose the combined forces of surface tension and the net gluonic pressure, even when the latter pressure is at its maximum (which is attained when the wall temperature is at its minimum, T wall = T ). At this radius, we say that the wall has suddenly attained mechanical equilibrium. We then switch over to a wall velocity of Eq. (A.30) for the rest of the simulation. We assume that the system maintains mechanical equilibrium and T wall = T until the end of the simulation. In Fig. 13 we plot the velocity and particle abundances within the pocket as functions of R for this new simulation. One can see that the velocity discontinuously drops at a radius of R ∼ 10 5 Λ −1 when the pocket abruptly achieves mechanical equilibrium. When this happens, further contraction is allowed only by subsequent quark annihilations; the quark depletion processes immediately recouple and dominate the density evolution due to the sharp drop in the contraction rate v w /R. At this same radius, the baryon abundance is at a maximum and we find that the DM relic abundance is set. Afterwards, the pocket slowly contracts and all particle numbers deplete until the pocket vanishes. Importantly, we find that the pocket asymmetry is saturated in this simulation (see Eq. (2.16)), and is also saturated for every other point in the DM parameter space that we consider. Of course, this is a crude approximation -a realistic scenario could have a smoother evolution of v w , or if the wall velocity falls sharply, this could induce plasma shock waves which modify the pocket evolution. However, this explicit example supports our intuition that turning on quark pressure will drive the system rapidly into the regime where the asymmetry is saturated, and once in this regime the details of the evolution do not affect the final relic abundance.

B Cross Sections
For the computation of the baryon survival factor, multiple processes are relevant. We distinguish between three classes of interactions: • Annihilation process, i.e. a direct annihilation of free quarks into dark gluons: 1 + (−1) → 0 + 0.
Similar to Tab. 1, in writing these equations we use each relic's quark number.
For the values of the annihilation and capture cross sections explicit calculations taking into account group theory factors have been performed in [65,27,28]. 15 The cross sections scale as σ ann./cap. v = ζ πα 2 m 2 q ≡ ζ σ 0 , (B.1) where ζ is a numerical factor that depends on the number of colors and flavors in the theory, and the coupling strength α is evaluated at the scale of the momentum transferred in the annihilation process, which is m q . In addition at low interaction energies the bound state formation and the annihilation process experience enhancement due to nonperturbative Sommerfeld corrections. The cumulative effect of those non-perturbative effects at finite temperature can be taken into account by the effective cross section σ eff v rel defined in [65,27,28]. In Fig. 14 the factor ζ is shown for different annihilation and capture processes in our full set of Boltzmann equations. Thermal masses of the dark gluons prevent bound state formation at large temperatures, an effect that has been confirmed by additional investigations of the process at hand in a non-equilibrium field theory treatment [65,121,122].
The rearrangement process is more complex and requires taking into account nonperturbative effects. Here simulations and comparisons to hydrogen−anti-hydrogen annihilation have been performed in [27]. The resulting cross section scales with the area set by the Bohr radius of the colliding bound states and contains a suppression factor which becomes effective once the kinetic energy exceeds the available binding energy: This results in a constant σv cross section, which is expected for an exothermic reaction with C N being the quadratic Casimir of the dark quark representation (C N = 4/3 for quarks in the fundamental representation of SU (N )), that controls the interquark attraction in a non-abelian theory. The overall scaling is thus σ RA v ∼ σ ann./cap. v/α 3 , in agreement with the numerical results of [123]. This non-perturbative enhancement results from taking into account the finite size of the colliding bound states.  Table 3: The processes and cross section classes involved in the annihilation and baryon formation. Bound states are donated by their baryon number. Direct annihilation takes place if multiple gluons need to be emitted in the final state. If one gluon is radiated, then perturbative capture cross section calculations apply. In the case that the final states have no gluons that would be needed for momentum conservation, we use the geometric rearrangement cross sections discussed below. The 0 in the rearrangement processes denotes a pionqq that promptly decays to gluons. Generally, all processes above have equivalent reactions, where all particles are replaced by anti-particles. For those we assume the same cross sections, i.e. assuming CP is conserved.
We summarize all the cross sections entering Eq. (3.10) in Tab. 3. Notice that in this table some processes involving gluons (denoted by 0) are listed as rearrangement. These processes, in fact, have an intermediate step in which the quarks rearrange and make pions (qq), which can promptly decay into gluons, see Eq. (3.1).

C Binding Energies
The binding energies of several types of dark states enter the full set of Boltzmann equations. For the two quark states exact results are available. Since we work in the limit m q Λ DC the Coulomb potential approximation is valid. For the baryon binding energy variational methods are needed, and numerical evaluations were performed in [27]. We focus on the case that N = 3. The resulting binding energies are: • Binding energy of a singlet diquark, or mesonqq: Eq q B = 1 4 α 2 C 2 N m q .
• Binding energy of a bound non-singlet diquark state qq in a binding configuration is: E qq B = 1 4 Eq q B .
• Binding energy of a baron qqq singlet sate is: E qqq B = 0.26α 2 C 2 N m q .
Here α is the gauge coupling of the confining group given by Eq. (2.4). The relevant scale at which the coupling should be evaluated is the Bohr momentum αm q , which can be determined iteratively, starting from the value of α(m q ). It has been shown that the corrections due to the linear (confining) part of the potential between quarks has negligible effect on these binding energies (see Eq. (4) of [28]).