String fragmentation in supercooled confinement and implications for dark matter

A strongly-coupled sector can feature a supercooled confinement transition in the early universe. We point out that, when fundamental quanta of the strong sector are swept into expanding bubbles of the confined phase, the distance between them is large compared to the confinement scale. We suggest a modelling of the subsequent dynamics and find that the flux linking the fundamental quanta deforms and stretches towards the wall, producing an enhanced number of composite states upon string fragmentation. The composite states are highly boosted in the plasma frame, which leads to additional particle production through the subsequent deep inelastic scattering. We study the consequences for the abundance and energetics of particles in the universe and for bubble-wall Lorentz factors. This opens several new avenues of investigation, which we begin to explore here, showing that the composite dark matter relic density is affected by many orders of magnitude.


Introduction
The possible existence of new confining sectors is motivated by most major failures of our understanding of Nature at a fundamental level. First, the stability of particle Dark Matter can be elegantly achieved as an accident if it is a composite state of a new strongly-coupled sector, similarly to proton stability in QCD, see e.g. [1]. The hierarchy problem of the Fermi scale is solved via dimensional transmutation by new confining gauge theories, whose currently most appealing incarnation is that of composite Higgs models [2,3]. Analogous composite pictures can UV-complete [4][5][6] twin-Higgs scenarios [7], and so ameliorate also the little hierarchy problem. A rationale to understand the SM hierarchies of masses and CKM mixing angles is provided by partial compositeness of the SM fermions [8]. Finally, new confining sectors play crucial roles in addressing the strong CP problem [9,10], the baryon asymmetry [11,12], etc.
Given their ubiquity, it makes sense to look for predictions of confining sectors that do not depend on the specific way they address a given SM issue. Cosmology naturally offers such a playground, in association with the confinement phase transition (PT) in the early universe. The low-density QCD phase transition would for example be strongly first-order if the strange or more quarks had smaller masses [13], with associated signals in gravitational waves [14,15]. New confining sectors could also well feature a similar PT. In addition, the confinement transition could be supercooled, a property that for example arises naturally in 5-dimensional (5D) duals of 4D confining theories [16][17][18].
Generically, supercooling denotes a PT in which bubble percolation occurs significantly below the critical temperature. Here we are interested in the case where a cosmological PT becomes sufficiently delayed so that the radiation energy density becomes subdominant to the vacuum energy. The universe then experiences a stage of inflation until the PT completes [19]. This implies a dilution of any pre-existing relic, such as dark matter (DM), the baryon or other asymmetries, topological defects, and gravitational waves, see e.g. [20][21][22].
In this paper we point out an effect that, to our knowledge, had been so far missed: when the fundamental quanta of the strong sector enter the expanding bubbles of the confined phase, their relevant distance can be much larger than the inverse of the confinement scale, thus realising a situation whose closest known analogues are perhaps QCD jets in particle colliders or cosmic ray showers. We anticipate that our attempt to model this phenomenon implies an additional production mechanism of any composite resonancestring fragmentation followed by deep inelastic scattering -which introduces a mismatch between the dilution of composite and other relics. This opens new model building and

JHEP04(2021)278
phenomenological avenues, which we begin exploring here in a model independent manner for the case of composite DM. The application of our findings to a specific model, namely composite dark matter with dilaton mediated interactions, will appear elsewhere [23].

Synopsis
Due to the numerous effects which will be discussed in the following sections, it is perhaps useful for the reader that we summarise the overall picture in a few paragraphs. We begin in the deconfined phase in which the techniquanta TC of the new strong sector (which we will call quarks and gluons) are in thermal equilibrium. Their number density normalised to entropy takes a familiar form where g TC (g s ) are the degrees of freedom of the quarks and gluons (entropic bath) respectively. Next a period of supercooling occurs, in which the universe finds itself in a late period of thermal inflation, which is terminated by bubble nucleation. As is known from previous studies, such a phase will dilute the number density of primordial particles. The dilution factor is given by where T nuc is the nucleation temperature, T start ∝ f is the temperature at which the thermal inflation started, T RH is the temperature after reheating, and f is the energy scale of confinement. We assume reheating to occur within one Hubble time, so that T RH ∝ f . The supercooled number density of quarks and gluons then becomes For completeness, the details entering eq. (2.3) will be rederived in section 3. When the fundamental techniquanta are swept into the expanding bubbles, they experience a confining force. Because f T nuc in the supercooled transition, the distance between them is large compared to the size of the composite states ψ (which we will equivalently call 'hadrons'). The field lines attached to a quark or gluon then find it energetically more convenient to form a flux tube oriented towards the bubble wall, rather than directly to the closest neighbouring techniquantum, which is in general much further than the wall (see figure 2). The string or flux tube connecting the quark or the gluon and the wall then fragments, producing a number of hadrons inside the wall. Additionally, because of charge conservation, techniquanta must be ejected outside the wall to compensate (see figure 3). The process is conceptually analogous to the production of a pair of QCD partons at colliders, and we model it as such. The details are explained in section 4. The result is an increase of the yield of composite particles, compared to the naive estimate following directly from eq. (2.3), by a string fragmentation factor K string ,

JHEP04(2021)278
where γ wp > f /T nuc 1 is the Lorentz factor of the bubble wall at the time the quarks enter.
The Lorentz factor is estimated in section 5. In section 6 we show that our picture can be relevant already for T nuc /T start 1. The quarks ejected from the bubbles are treated in detail in section 7. We find they enter neighbouring bubbles and confine there into hadrons. Acting as a cosmological catapult, string fragmentation at the wall boundary gives a large boost factor to the newly formed hadrons, such that their momenta in the plasma frame can be f . The composite states and their decay products can next undergo scatterings with other particles they encounter, e.g with particles of the preheated 'soup' after the bubbles collide. Since the associated center-of-mass energy can be much larger than f , the resulting deep inelastic scatterings (DIS) increase the number of hadrons. We explore this in detail in section 8. The resulting effect on the yield can be encapsulated in a factor K DIS , and reads where M Pl is the Planck mass and m * = g * f is the mass scale of hadrons. The last proportionality holds in the regime of runaway bubble walls, relevant for composite DM.
Finally the late-time abundance of the long-lived and stable hadrons, if any, evolves depending on their inelastic cross section in the thermal bath σv rel , and on Y SC+string+DIS ψ as an initial condition at T RH . We compute it in section 9 by solving the associated Boltzmann equations.
By combining all the above effects we arrive at an estimate of the final relic abundance of the composite states. Our findings impact their abundance by several orders of magnitude, as can be seen in figure 9 for the concrete case where the relic is identified with DM. The formalism leading to this estimate can readily be adapted for other purposes. For example, if ψ instead decays out-of-equilibrium, it could source the baryon asymmetry. The estimate of Y SC+string+DIS ψ would then act as the first necessary step for the determination of the baryonic yield.

Strongly coupled CFT
Although striving to remain as model independent as possible in our discussion, we shall be making a minimal assumption that the confined phase of the strongly coupled theory can be described as an EFT with a light scalar χ, e.g. a dilaton. The scalar VEV, χ , then parametrizes the local value of the strong scale. It can be thought of as a scalar condensate of the strong sector, such as a glueball-or pion-like state. The scalar VEV at the minimum of its zero-temperature potential is identified with χ = f , where f is the confinement energy scale, while χ = 0 at large enough temperatures. In order to have strong supercooling, we require the approximate (e.g. conformal) symmetry to be close to unbroken, thus justifying the lightness of the associated pseudo-Nambu-Goldstone boson (e.g. the dilaton [24]). That supercooling occurs with a light dilaton is known from a

Thermal history
The vacuum energy before the phase transition is given by with some model dependent c vac ∼ O(0.01) constant. The radiation density is given by where g R counts the effective degrees of freedom of the radiation bath. We define g R ≡ g Ri (g Rf ) in the deconfined (confined) phase. Now consider the case of strong supercooling. The universe will enter a vacuum-dominated phase at a temperature T start = 30 c vac g Ri π 2 1/4 f, (3.3) provided the phase transition has not yet taken place beforehand. The vacuum domination signals a period of late-time inflation. The phase transition takes place at the nucleation temperature, T nuc , when the bubble nucleation rate becomes comparable to the Hubble factor. Following the phase transition, the dilaton undergoes oscillations and decay, reheating the universe to a temperature

4)
At this point the universe is again radiation dominated. We have assumed the decay to occur much faster than the expansion rate of the universe such that we can neglect a matter-dominated phase [20].

Dilution of the degrees of freedom
Now consider some fundamental techniquanta of the strong sector, e.g. techniquarks or technigluons (for simplicity we always refer to them as quarks and gluons). Prior to the phase transition the number density of techniquanta follows a thermal distribution for massless particles where g TC denotes the degrees of freedom of the quanta under consideration. The entropy density is given by

JHEP04(2021)278
where g s are the total entropic degrees of freedom. 1 The number density normalized to entropy before the phase transition, remains constant up to the point when the phase transition takes place. The entropy density then increases during reheating giving when we find ourselves back in the radiation-dominated phase. The dilution factor from the additional expansion during the vacuum-dominated phase can be derived by finding the increase in entropy between T nuc and T RH . It reads If the quarks and gluons were non-interacting following the phase transition, the yield today would be given by the above formula. (In the presence of interactions the above would be taken as an initial condition at T RH for the Boltzmann equations describing the effects of number changing interactions between reheating and today.) The picture would then be analogous to that studied, in a theory without confinement, in [20]. The picture is completely changed, however, for supercooled confining phase transitions, which we elucidate next.

Where does confinement happen?
Bubble wall profile. The expanding bubble is approximately described by the Klein-Gordon equation [44] where s 2 = t 2 − r 2 is the light-cone coordinate and V is the scalar potential. A sketch of a typical bubble profile for close-to-conformal potentials is shown in figure 1. The key point here is that the wall thickness is as shown by numerical computations and analytical estimates, see appendix A for a calculation in an explicit example.

JHEP04(2021)278
χ f s χ * ∼ 1/T n Figure 1. A typical wall profile found in close-to-conformal potentials. After nucleating by tunneling to the exit point, χ * f , the field rolls down and undergoes damped oscillations around the minimum of its potential. The typical wall thickness is L w ∼ 1/T nuc .

Confinement time scale.
The techniquanta (quarks and gluons) constitute a plasma with temperature of order T nuc before entering the bubble. Once they enter the bubbles, they could in principle either confine in a region close to the bubble wall where χ f , or approach as free particles the region where χ has reached its zero-temperature expectation value χ = f . To determine this, let us define a 'confinement rate' and a 'confinement length' as where n TC and v TC are, respectively, the number density and the relative Møller velocity of the techniquanta v TC ≡ |v 1 − v 2 | 2 − |v 1 × v 2 | 2 1/2 [45], and σ conf is a 'confining cross section'. We want to compare L conf with the length of the bubble wall, defined as the distance over which χ varies from its value at the exit point, χ = χ * f , to χ = f . Of course we need to perform this comparison in the same Lorentz frame, so we emphasise our definition of L w as the bubble-wall length in the bubble-wall frame, and L p = L w /γ wp as the bubble-wall length in the frame of the center of the bubble, which coincides with the plasma frame, and where γ wp is the boost factor between the two frames. Let us now move to the confinement timescale of eq. (4.3). Since we expect confinement to happen 'as soon as possible', we assume the related cross section to be close to the unitarity limit [46], Since n 2 TC v TC is Lorentz invariant [45,47], one then has that n TC v TC transforms under boosts as n −1 TC . The boost to apply in this case is γ wp , because by definition the string forms after confinement, so we can treat the plasma frame as the center-of-mass frame of the techniquanta. Combining this with the Lorentz invariance of the cross section, we obtain where in the last equality we have used that the average relative speed and density of the techniquanta in the plasma frame satisfy, respectively, v TC, p 1 and n TC, p ∼ T 3 nuc , because JHEP04(2021)278 they are relativistic. This in turn implies Confinement takes place deep inside the bubble. For the regimes of supercooling we are interested in, the phase transition is of detonation type and the Lorentz factor γ wp is orders of magnitude larger than unity. Therefore, L conf, w L w such that confinement does not happen in the outermost bubble region where χ f . This conclusion is solid in the sense that it would be strengthened by using a confinement cross section smaller than what assumed in eq. (4.4), which is at the upper end of what is allowed by unitarity. The end effect of the above discussion, is that for practical purposes, we can consider the wall profile to be a step-like function between the deconfined phase, χ = 0, and confined phase, χ = f . Furthermore, as we shall discuss below, the quarks will not confine directly in pairs but rather form fluxtubes pointing toward the bubble wall as they penetrate the χ = f region of the bubble.
The ballistic approximation is valid. Equation (4.6), together with the large wall-Lorentz-factors encountered in this study, implies that we can safely neglect the interactions between neighbouring techniquanta during the time when they cross the bubble wall. This is the so-called ballistic regime, see e.g. [48], which will be useful for deriving the friction pressure in section 5.

Fluxtubes attach to the wall following supercooling
A hierarchy of scale. Upon entering the region χ = f of expanding bubbles, the techniquanta experience a confinement potential much stronger than in the region close to the wall. This can be easily understood by taking the long-distance potential of the Cornell form [49][50][51][52][53][54][55][56][57][58] where d c is the techniquanta seperation in their 'center of interaction frame' (or equivalently 'string center of mass frame'), 2 and c TC is an adimensional constant, 3 c qq 10 in QCD [58]. A crucial point regarding the string energy in this context, besides the fact it grows proportionally to χ 2 , is that the inter-quanta distance is large compared to the natural confinement scale, i.e. d c f −1 , due to the supercooling. Indeed the distance between quanta outside the wall, in the plasma and wall frames respectively, scales as outside the wall, and because the quarks and gluons cannot be accelerated upon entering so d w is Lorentz contracted with respect to d p ), one ends up with d c f −1 . What happens then to the techniquanta and to the fields connecting them? 2 Lattice simulations find that the QCD potential at dc fm saturates to a constant, a behavior which is interpreted in terms of pair creation of quarks from the vacuum, see e.g. the recent [59]. Therefore this realises an outcome that, for our purposes, coincides with having ETC ∝ dc to larger distances. Lattice simulations with quarks only as external sources [60], so without sea quarks ('quenched'), find that the linear regime of the QCD Cornell potential extends up to the maximal distances probed, namely dc 3 fm in the results reported in [60]. 3    Flux tubes minimize their energy. In a picture without hierarchy of scales, the fields would compress in fluxtubes connecting different charges, 'isolated' in pairs or groups to form color-singlets. Here, we argue that the fluxtubes have another option, which is energetically preferable: that of orienting themselves towards the direction of minimal energy, i.e. as perpendicular as possible to the bubble-wall, 4 and to keep a 'looser' connection in the outer region where χ f . Indeed, a straight-line connection between techniquanta would result in a much longer portion of fluxtubes in the region χ = f , with respect to our picture of fluxtubes perpendicular to the wall. Via eq. (4.7), this would in turn imply a much higher cost in energy, disfavoring that option. We stress that, in our picture, the fluxtubes are still connecting techniquanta in such a way to form an overall color singlet, just these fluxtubes minimise their length in the region χ = f , and partly live in a region χ χ * f . This picture is visualised in figure 2. Note the nearest neighbour quark from the plasma may also be located outside the bubble.

Condensed matter analogy.
An interesting analogue to the picture above is the vortex string of magnetic flux in the Landau-Ginzburg model of superconductivity. To match onto confinement dynamics a dual superconductor is pictured, in which the external colourelectric field -rather than the magnetic field -is expelled by the Meissner effect [61]. Here the bubble of confining phase corresponds to the superconductor from which the colourelectric field is expelled. Quarks entering the bubble then map onto magnetic monopoles being fired into a regular superconductor.

String energy and boost factors
To possibly be quantitative on the implications of the picture we just outlined, we first need to determine the string energy and the Lorentz boosts among the frames of the plasma, wall and center-of-mass of the string.
String end-points. Let us define as TC i the quark or gluon that constitutes an endpoint, inside the bubble, of a fluxtube pointing towards the wall, and the end-point of the fluxtube on the wall. The energy of the incoming techniquantum in the wall frame is E i,w = 3γ wp T nuc , where for simplicity we have averaged over their angle with respect to the wall. We assume to be at rest or almost, and to carry some O(1) fraction of the inertia of the string. Hence the respective four-momenta are String center-of-mass. Then we define the center-of-mass of the string as the one of TC i and , and find where the second expression is valid up to relative orders (γ wp f /T nuc ) −1 1. By employing a Lorentz boost between the wall and center-of-mass frames, and imposing p i,c = − p ,c , we find On the right-hand side of the equations above we have omitted a factor of 2(q − ), in (4.9), and of 1/ 2(q − ), in (4.10), because for simplicity we take these to be ≈ 1 from now on (as per the benchmark q = 1/2, = 0). Finally we determine the boost between the center-of-mass frame of the string and the plasma frame as which is valid up to a relative order (γ wp f /T nuc ) −1 1.

Hadrons from string fragmentation: multiplicity and energy
The fluxtubes connecting a quark or gluon to the wall will fragment and form hadrons, singlet under the new confining gauge group. We would now like to determine: • The number of hadrons formed per fluxtube.
• The momenta of said hadrons.  Figure 4. Left: average hadron multiplicity per single QCD scattering e − e + → qq. Right: root square mean of the hadron energy per single QCD scattering e − e + → qq, whereĒ hadr = E CM / N hadr is the average hadron energy per scattering. Dots in both plots are extracted via MadAnalysis v1.8.34 [63] by simulations with MadGraph v2.7.0 [64] plus Pythia v8.2 [65], the line in the left-hand plot displays eq. (4.14). Results are expressed as a function of the center-of-mass energy of the scattering in GeV, to export them to our cosmological picture we simply substitute GeV 4πf π by m * = g * f , and use E CM = 3 γ wp T nuc f .

Collider analogy.
We start by noticing that the process of formation of a fluxtube, in our picture, is analogous to two color charges in an overall-singlet state, TC i and , moving apart with a certain energy E CM , where E CM = 3 γ wp T nuc f in the modelling of section 4.3. This physical process appears entirely analogous to what would happen in a collider that produces a pair of techniquanta of the new confining force, starting from an initial singlet state. In light of this observation, we then decide to model the process by analogy with a very well-studied process observed in Nature, that of QCD-quark pair production at electron-positron colliders, where the analogy lies also in the fact that the initial state electron-positron pairs is in a color singlet state. Needless to say, a BSM confining sector needs not behave as QCD in terms of number and momenta of hadrons produced per scattering, see e.g. [62]. However, QCD constitutes a well studied and tested theory, so that we find it reasonable to use it as our benchmark. Moreover, we anticipate from section 8 that our final result for the cosmological abundance of hadrons, in the assumption of efficient-enough interactions between them and the SM, will only depend on the initial available energy E CM . This suggests that, within that assumption, our final findings hold for confining sectors that distribute this energy over a number of hadrons different from QCD.
Numerical simulations. We use Pythia v8.2 [65] interfaced to MadGraph v2.7.0 [64] to simulate the process e − e + → qq for different center-of-mass energies, and MadAnalysis v1.8.34 [63] to extract from these simulations both the total number of hadrons produced per scattering and their energy distribution. We thus recover known QCD results and display them in figure 4. We translate them to our picture by replacing the units of a GeV 4πf π used by Pythia, with the generic mass of a composite state m * = g * f , where 1 ≤ g * ≤ 4π is some strong effective coupling. These results can be summarised as follows: • The number of hadrons produced per fluxtube grows logarithmically in E CM .

JHEP04(2021)278
• The distribution of hadron energies is such that its root square mean coincides, to a percent level accuracy, with the average energy per hadron This will support, in section 8.3, our simplifying assumption that all hadrons produced by the string fragmentation carry an energy of orderĒ hadr .
Results from the literature. The multiplicity of QCD hadrons from various scattering processes has been the object of experimental and theoretical investigation, since the late 1960s [66]. We now leverage such studies both to check the results of our simulation and to obtain analytical control over them. Collider studies have typically focused on the multiplicity of charged QCD resonances per scattering, n ch . In particular, works such as [67,68] have carried out the exercise of collecting the most significant measurements of n ch and 'filling' the missing phase space -not covered by detectors -with the output of MC programs, thus obtaining a full-phase-space quantity. We take as our starting point the result provided in [68] from pp collisions, which reads with (a, b, c, d) = (0.95, 0.37, 0.43, 0.04). Here, as already explained, we substituted the normalisation of a GeV with m * = g * f .
Our modelling. To obtain the total number of hadrons from e + e − collisions we proceed as follows. First, most hadrons coming out from hard scatterings consist in the lightest ones, i.e. the pions. Second, the total number of pions produced is very well approximated by 3 n ch /2, because of isospin conservation. By the first argument, this coincides with very good approximation to the total number of hadrons produced. Third, the multiplicity of composite states from e + e − collisions has been found to roughly match the one from pp collisions, upon increasing the e + e − energy by a factor of 2, see e.g. section 2.2 in [69]. 5 We then model the total number of composite states produced, per string fragmentation, as where we have multiplied by an exponential and added one to smoothen N string ψ (E CM ) to 1 as E CM → m * , because this physical regime was not taken into account in [68]. In the left-hand panel of figure 4 one sees that eq. (4.14) reproduces the results of our Pythia simulation for E CM smaller than a few TeV rather well. This was to be expected since eq. (4.13) was determined in [68] from fits to data up to that energy. It is not the purpose of this paper to improve on this fit, as stated above, we simply use the above results as a check of our Pythia simulation.

Enhancement of number density from string fragmentation
Production of composite states. Prior to (p)reheating, we then have a yield of composite states given by the yield of strings, which can be estimated from eq. (3.8), multiplied by the number of composite states per string light composite state , (4.15) where N string ψ (E CM ) is given by eq. (4.14) and E CM = 3 γ wp T nuc f in eq. (4.9). We have distinguished the cases where the composite state of interest is heavier or lighter than the glueballs (e.g. the analogous of a proton or a pion in QCD). In the former case, the −1 we added to the factor multiplying g g accounts for the fact that, if the final composite states produced by string fragmentation do not undergo other additional interactions, then glueballs decay to the light composite states and do not contribute to the final yield of any heavy composite state of quarks. The yield of composite states ψ then reads The appearance of Y eq TC in eq. (4.16) accounts for string formation from both quarks and gluons. Hence, not only is the number of ψ's enhanced by the string fragmentation, relative to the case with no confinement, but also by the possibility of gluons to form strings. K string and Y SC+string ψ are plotted in figure 8.
Hadrons are highly boosted in the plasma frame. The hadrons formed after string fragmentation schematically consist of two equally abundant groups. Hadrons in the first group, which for later convenience we call 'Population A', move towards the bubble wall with an average energy where we have boosted the energy per hadron of eq. (4.12) to the plasma frame with the γ cp of eq. (4.11), and also used eqs. (4.9) and (4.15). We conclude from eq. (4.17) that the newly formed hadrons have large momenta in the plasma frame. The formation of a gluon string between the incoming techniquanta and the wall acts as a cosmological catapult which propels the string fragments in the direction the wall is moving. Hadrons in the second group move, in the wall frame, towards the bubble wall center, and their energy in the plasma frame is negligible compared to (4.17). Note that if only one hadron is produced on average per every string, then it would roughly be at rest in the center-of-mass frame of the string, with an energy (mass) of order E CM . In the plasma frame, its energy would then read E p γ cp E CM γ wp f /2. As we will see in section 8, the impact of this hadron on the final yield would then be captured by our expressions. Following this first stage of string fragmentation, the composite states, and/or their decay products, can undergo further interactions with remnant particles of the bath, preheated or reheated plasma, and among themselves. Such interactions may change the ultimate yield of the relic composite states. Before taking these additional effects into account JHEP04(2021)278 in section 8, in the next sections we complete the modelling we proposed above, by describing the behaviour of the ejected quarks and deriving the Lorentz factor of the wall, γ wp .

Ejected quarks and gluons and their energy budget
So far we dealt with what happens inside the bubble wall. The process we described apparently does not conserve color charge: we started with a physical quark or gluon with a net color charge entering the bubble, and we ended up with a system of hadrons which is color neutral. Where has the color charge gone?
The necessity of ejecting a quark or gluon. To understand this, it is convenient to recall the physical modelling behind the process of string fragmentation that converts the initial fluxtube into hadrons, see e.g. the original Lund paper [70]. When the fluxtube length, in its center-of-mass frame, becomes of order f −1 , the string breaks at several points via the nucleation of quark-antiquark pairs from the vacuum. Now consider, in our cosmological picture, the quark-antiquark pair nucleated closest to the bubble wall. One of the two -say the antiquark -forms a hadron inside the wall. The only thing that can happen to the quark is for it to be ejected from the wall, because of the lack of charge partners inside the wall. This process, somehow reminiscent of black hole evaporation, thus allows for charge to be conserved. The momentum of the ejected quark, in the wall frame, has to be some order-one fraction of the confinement scale f , because that is the only energy scale in the process. For definiteness, in the following we will take this fraction to be a half. This picture is visualized in figure 3, and it is analogous if TC i is a gluon instead of a quark.
Energy of the ejected quark or gluon. One then has one ejected quark (at least) or gluon per fluxtube, thus per quark or gluon that initially entered. Therefore, the number of techniquanta outside the bubble wall does not diminish upon expansion of the bubble. This population of ejected techniquanta is energetically as important as that of hadrons inside the bubble. Indeed the energy of an ejected quark or gluon (or quark pair), in the plasma frame, reads E ej,p γ wp f. (4.18) This is of the same order as the total energy in the hadrons from the fragmentation of a single string, obtained by multiplying E A,p of eq. (4.17) times half of the total number of hadrons produced per string (i.e. we included only the energetic ones). The population of ejected techniquanta cannot therefore be neglected in the description of the following evolution of this cosmological system.

Bubble wall velocities
The wall boost in the plasma frame, γ wp , affects many key properties of our scenario, from the ejection of techniquanta to the number and energy of the hadrons produced by string fragmentation. It is the purpose of this section to study the possible values it can take over the PT.

JHEP04(2021)278
Final results. As bubbles are nucleated and start to expand, γ wp starts growing as well.
If nothing slows down the bubble-wall acceleration, then γ wp keeps growing until its value at the time of bubble-wall collision, γ runaway wp . Sources of friction that could prevent this runaway regime are given by the equivalent, in this scenario, of the so-called leading order (LO) and next-to-leading order (NLO) contributions of [71] and [72] respectively. We find it convenient to report right away our final result for the maximal possible value of γ wp , where the first entry is associated to γ runaway wp , and the second to the boost as limited by the LO pressure, γ LO wp . γ LO wp is always smaller than γ NLO wp in the parameter space of our interest, so that γ NLO wp does not enter eq. (5.1). We learn that in the regime of very strong supercooling and/or of very large confinement scale f , which will be the most relevant one for the DM abundance, bubble walls run away. The behaviour of γ wp is illustrated in figure 5.
The impact on GW. The behaviour of γ wp also has important consequences for the gravitational wave signal from the phase transition [73,74]. If γ max wp = γ runaway wp then the vacuum energy is converted into kinetic energy of the bubble walls [75]. The gravitational wave (GW) spectrum sourced by scalar field gradient is traditionally computed in the envelope approximation [76][77][78]. However, the latest lattice results [79,80] suggest an enhancement of the GW spectrum at low frequency due to the free propagation of remnants of bubble walls after the collision, the IR slope ∝ k 3 becoming close to ∝ k 1 . This confirms the predictions from the analytical bulk flow model [81,82]. Note that the IR-enhancement is stronger for thick-walled bubbles [79], which is the case relevant for nearly-conformal potential leading to strong supercooling, and thus for the PT considered here. (Instead, for thin-walled bubbles, after collision the scalar field can be trapped back in the false vacuum [11,44,83]. Instead of propagating freely, the shells of energy-momentum tensor remain close to the collision point and dissipate via multiple bounces of the walls.) Irrespectively of whether the IR slope at f β is ∝ k 3 or ∝ k 1 , at much lower frequency, f H, the slope must converge to k 3 due to causality [84][85][86]. Oscillations of the condensate following the PT can provide an additional source of GW [87]. However, instead of β −1 the time scale is set by the inverse scalar mass ∼ f −1 and the signal is Planck-suppressed ∝ β/f [88].
, the vacuum energy is converted into thermal and kinetic energy of the particles in the plasma already prior to the bubble wall collision. The contribution from sound waves or turbulence [73,74], however, in supercooled transitions is not yet clearly understood. Indeed, current hydrodynamical simulations, which aim to capture the contribution of the bulk motion of the plasma to the gravitational wave signal, do not yet extend into the regime in which the energy density in radiation is subdominant to the vacuum [89]. And analytical studies of shock-waves in the relativistic limit have just started [90]. In any case, we expect supercooled transitions to provide promising avenue for detection in future GW observatories.
We now proceed to a detailed derivation of eq. (5.1). In this regime, larger f or smaller T nuc leads to a smaller distance over which the bubble can accelerate. The former because of the smaller Hubble horizon and the latter due to the larger bubble size at nucleation. Therefore γ wp decreases for more supercooling.

Linear growth. The energy gained upon formation of a bubble of radius
where ∆V vac is the difference between the vacuum energy density outside and inside the bubble. The energy lost upon formation of a bubble of radius R is E wall 4πR 2 γ wp σ w , where σ w is the surface energy density of the wall (surface tension) in the wall frame. If a bubble nucleates and expands, its energy E bubble is transferred to the wall energy E wall . As soon as a nucleated bubble contains the region χ f , neither ∆V vac nor σ w change upon bubble expansion. Indeed both are a function of the bubble wall profile, which does not change in that regime (also see figure 1). We thus recover the well-known property that γ wp grows linearly in R, where R 0 is a normalisation of the order of the minimal radius needed for a bubble to nucleate, and where in the second relation we have used R 0 L w ∼ T −1 nuc because we assumed the nucleated bubble to contain the region χ f . A more precise treatment can be found, e.g. in the recent [75], which confirms the parametric dependence of eq. (5.2).

At collision time.
In a runaway regime, i.e. for small enough retarding pressure on the bubble walls, γ wp at collisions then reads where β −1 is the average radius of bubbles at collision, H Λ 2 vac /( √ 3M Pl ), and the value β/H 10 is a benchmark typical of supercooled phase transitions [17,23,27,35,36,41,91], which we employ from now on.

JHEP04(2021)278
The bubbles swallow most of the volume of the universe, and thus most techniquanta, when their radius is of the order of their average radius at collision β −1 . Therefore, in the regime of runaway bubble walls, the relevant γ wp for all the physical processes of our interest (hadron formation from string fragmentation, quark ejection, etc.) will be some order one fraction of γ runaway wp . For simplicity, in the runaway regime we will then employ the simplifying relation γ wp = γ runaway wp . This will not only be a good-enough approximation for our purposes, but it will also allow to clearly grasp the parametric dependence of our novel findings. Moreover, a more precise treatment, to be consistent, would need to be accompanied by a more precise solution for γ wp than that of eq. (5.2), i.e. we would need to specify the potential driving the supercooled PT and solve for γ wp . As the purpose of this paper is to point out effects which are independent of details of the specific potential, we leave a more precise treatment to future work.

LO pressure
Origin. By LO pressure we mean the pressure from the partial conversion -of the quark's momenta before entering the bubbles -into hadron masses [71], plus that from the ejection of quarks. We use the subscript LO in reference to [71,72], because this pressure is of the form P LO ∼ ∆m 2 T 2 , where ∆m is the rest energy of the flux tube between the incoming techni-quanta and the wall. However, in contrast to [71,72], here the pressure arises from non-perturbative effects.
Momentum transfer. The momentum exchanged with the wall, upon hadronization of a single entering quark plus the associated quark ejection, reads in the wall frame where E in 3 γ wp T nuc is the energy of the incoming quark, ∆m 2 in is the fraction of that energy that is converted into 'inertia' of the string, and E ej f /2 is the energy of the ejected quark or gluon. In the second equality, we have used ∆m 2 in E 2 CM 3 γ wp T nuc f from eq. (4.9) and γ wp f /T nuc . Note that ∆p LO is independent of p in .
Pressure. In light of section 4.1, we can safely consider a collision-less approach and neglect the interactions between neighboring quarks. The associated pressure is given by where g a is the number of internal degrees of freedom of a given species a of the techniquanta. Upon using eq. (5.4), we get where we remind that g TC = g g + . This result can be understood intuitively from P LO ∼ n TC,w ∆p LO , where γ wp enters through n TC,w [71]. Note that, in the absence of ejected particles, the pressure would have been a half of our result in eq. (5.6).

JHEP04(2021)278
Terminal velocity. The resulting upper limit on γ wp is obtained by imposing that the LO pressure equals that of the internal pressure from the difference in vacuum energies, and reads (5.8) We finally remark that P LO grows linearly in γ wp , unlike in 'standard' PTs where it is independent of the boost. The reason lies in the fact that the effective mass ∆m 2 in grows with γ wp , whereas in 'standard' PTs it is constant in γ wp . Our results then imply that, in confining phase transitions, the LO pressure is in principle enough to ensure the bubble walls do not runaway asymptotically. This is to be contrasted with non-confining PTs, where the asymptotic runaway is only prevented by the NLO pressure. 6

NLO pressure
Origin. The NLO pressure comes from the techniquanta radiating a soft gluon [72] which itself forms a string attached to the wall in the broken phase.
Result. We derive it in detail in appendix B. We find, cf. eq. (B.20) where C 2 [g, q] are the second Casimirs of the representations of gluons and quarks under the confining group (if SU(N ), , g conf is the gauge coupling of the confining group, ps ≤ 1 encodes the suppression from phase-space saturation of the emitted soft quanta g, important for large coupling g conf , m g is an effective mass of the soft radiated gluons responsible for this pressure, and k * the IR cut-off on the momentum radiated in the direction parallel to the wall.
Vector boson mass. As we model the masses of our techniquanta as the inertia that their fluxtube would gain inside the bubble, these masses increase with increasing momentum of the techniquanta, in the wall frame. The NLO pressure is caused by emission of gluons 'soft' with respect to the incoming quanta. Their would-be mass m g upon entering the wall cannot, therefore, be as large as that of the incoming quanta that emit them, ∆m in 3 γ wp T nuc f . At the same time, the effective gluon mass should at least allow for the formation of one hadron inside the wall, therefore we assume it to be of the order of the confinement scale, m g ∼ f . The fact that m g does not grow with γ wp while ∆m in does, is the reason why unlike in non-confining phase transitions, we find here that P NLO and P LO have the same scaling in γ wp and in the amount of supercooling.

JHEP04(2021)278
NLO pressure is sub-leading. By making the standard [72] choice k * m g , and assuming ps g 2 conf < 1, we then find that P NLO P LO in the entire parameter space of our interest. Thus, for simplicity, we do not report the NLO limit on γ max wp in eq. (5.1). Recently, ref. [92] performed a resummation of the log-enhanced radiation that leads to the scaling P NLO ∝ g 2 conf γ 2 wp T 4 nuc . By using the analogue of that result for confining theories, we find that P NLO dominates over P LO in some region of parameter space, and therefore that the values of the parameters for which bubble walls run away slightly change. Still, even by using that resummed result, we find that the region relevant for DM phenomenology corresponds to the region where bubble walls run away, so that the difference between the results of [72] and [92] does not impact the DM abundance. As observed in [93], the pressure as determined in [92] does not tend to zero when the order parameter of the transition goes to zero, casting a shadow on that result. Therefore, both for this issue as well as for the limited impact on the DM abundance that we will discuss later, we content ourselves with a treatment analogous to [72] in our paper.
Summary and runaway condition. At small supercooling (i.e. not too small T nuc /f ) the bubble wall velocity reaches an equilibrium value set by the LO pressure. At larger supercooling bubble walls collide before reaching their terminal LO velocity, and γ wp is set by the runaway value eq. (5.3). By comparing eq. (5.8) with eq. (5.3), we find that bubble walls run away for The bubble wall Lorentz factor is plotted in figure 5 against the amount of supercooling.

Ping-pong regime
Condition to enter. For even a single hadron to form inside the bubble, one needs E CM ≥ m π , where π is the lightest hadron of the new confining sector (e.g. a pseudogoldstone boson). Via eq. (4.9), this implies Contribution to the pressure. For γ wp γ enter wp , which holds at least in the initial stages of the bubble expansion, the quarks and gluons are reflected and induce a pressure This is to be compared with eq. (5.7), P expand = c vac f 4 , which implies the bubble wall could in principle be limited by this pressure to , this pressure ceases to exist at an earlier stage of the expansion, namely once γ wp = γ enter wp . Hence the maximum Lorentz factor remains encapsulated by eq. (5.1).
Ping-pong regime. In some extreme regions of parameter space, however, one could have γ max wp < γ enter wp , so that all techniquanta in the plasma are reflected at least once before entering a bubble. We leave a treatment of this 'ping-pong' regime to future work.

JHEP04(2021)278
6 Amount of supercooling needed for our picture to be relevant Intuition about the limit of no supercooling. In the limit of no supercooling, one does not expect the fluxtubes to attach to the bubble wall, but rather to connect the closest charges that form a singlet and induce their confinement. In other words, in the limit of no supercooling one expects the picture of confinement to be the one of 'standard phase transitions'. By continuity, there should exist a value of T nuc , smaller than f , such that the our picture ceases to be valid, and one instead recovers the more familiar confinement among closest color charges. We now wish to determine it. In order to do so, we note that the absence of ejected techniquanta is a necessary condition for the above to hold, therefore we now phrase the problem in terms of absence of ejected techniquanta.
Rate of detachment of . We propose and analyse some effects that could lead to fluxtubes detaching from the bubble walls without ejecting particles. To take place, these effects need to happen before the end-point of the fluxtube on the wall, , ceases to exist, i.e. when the string breaking inside the bubble has already taken place and a quark is ejected. So we start by computing the rate Γ det of detachment of , the point where the fluxtubes is attached to the wall, from the wall itself. To estimate it, we again borrow the modelling of the classic paper on string fragmentation [70].
The distances between the several points of breaking of a given string (that connects in our case TC i and ) are space-like. In the frame of each point of breaking, that breaking is itself the first to happen, a time of order N/f after the string formation (we adopt the scaling for strong sector gauge groups SU(N ) [94,95]). This time therefore also applies to the outermost breaking point in our picture, i.e. that closest to the wall, whose frame approximately coincides with the wall frame. We remind the reader that the outermost breaking is the one that nucleates the quark or gluon that is eventually ejected. The rate we need can therefore be estimated as the inverse of the nucleation time of the outermost pair, We now enumerate and model effects that could lead to fluxtubes detaching from the bubble walls without ejecting techniquanta, and compare their time scales with eq. (6.1).

Flux lines overlap.
The faster a bubble-wall, the denser and thus the closer together in the wall frame are the quarks and gluons entering it. Eventually, they could get closer than the typical transverse size of a fluxtube d tr f −1 [96]. When that happens, the fluxtubes between different color charges have a non-negligible overlap. We expect that in this situation it will not be clearly preferable energetically for these strings to attach directly to the wall. Thus there would be no ejected techniquanta. This situation is of course realised also in the case of small supercooling f /T nuc , in addition to and independently of the case of fast bubble-walls.
We then obtain a rate of 'string breaking by fluxtube overlap', Γ overlap , as follows. We define an effective associated cross section as the area of a circle on the wall, centered on any and with radius d tr ,

JHEP04(2021)278
The associated rate then reads where n TC,w = γ wp n TC,p is the density of techniquanta in the wall frame, g TC = g g + 3g q /4, and we have used that they are relativistic v = 1. The condition of no ejected techniquanta then reads The entire fluxtube connecting real color charges, so including its portion in the region χ χ * f (see figure 2), could enter the region χ = f before its portions in the region χ = f break and form hadrons, and eject particles. We see two ways this could happen.

Attractive interaction between neighboring flux lines. The points
are not static, because they move by the force exerted by the part of the string which is outside the wall, in the layer where χ χ * . Defining y as the transverse distance, on the wall, between two points connected by a fluxtube, one has where, consistently with our previous treatments, we have assigned to an inertia m ∼ f . If y goes to zero in a time shorter than the breaking time τ det ,w ∼ N f −1 , then the two fluxtubes connect and become fully contained in the region χ = f before they break and form hadrons, and thus there are no ejected techniquanta. To determine this condition, we assume initially static points , and thus we only need the initial distance between them y (t = 0) (γ wp n TC,p ) −1/3 . We then obtain The resulting condition for no ejected quarks reads (6.7) 2.2 Limit of no distortion of the flux lines. When the string portion in the region χ χ * f has a small enough length d , the possibility that it is pulled inside the region χ = f could be energetically more convenient than the one of our picture, where it stays outside and instead energy goes in increasing the length of the strings that are perpendicular to the wall. The energy price, for the string portion in the region χ χ * f to enter the region χ = f , reads in the wall frame

JHEP04(2021)278
where we stress that the length of the string portion d is transverse to the bubble-wall velocity and therefore is not Lorentz contracted in the process of being pulled into the bubble. In the wall frame, it reads d (γ wp n TC,p ) −1/3 The transition between χ χ * f and χ = f is exponentially fast in the proper coordinate s (see appendix A), and happens over an interval (a distance, in the wall frame) L f ∼ f −1 . The energy price of eq. (6.8) should therefore be compared with the one to stretch two strings, inside the wall, by an amount L f : where we have used that the string length in the expression for E qq , eq. (4.7), has to be evaluated in the string center-of-mass frame, and that γ wc 3 γ wp T nuc /f from eq. (4.10). Therefore, it is energetically more convenient to pull the fluxtube inside the region χ f , and so to have no ejected quarks, if Contrary to the previous two possibilities to have no ejected quarks, eqs. (6.4) and (6.7), the possibility in eq. (6.10) imposes an upper limit on γ wp . We anticipate that, in the regimes of supercooling interesting for our work T nuc /f 1, eq. (6.10) cannot be satisfied consistently with γ wp > 1, so that it is not relevant for our work.
Summary of required supercooling. In the regime where T nuc f , we expect that neither ejection of techniquanta nor string fragmentation should take place, and that the standard picture of quarks and gluons confining with their neighbors should be recovered (which we dub the 'standard phase transition'). More precisely, if any of eqs. (6.4), (6.7) and (6.10) hold, we depart from our picture in at least one regard. By demanding none of these inequalities hold, we expect the new effects of our study, namely flux line attached to the wall, string fragmentation, quark ejection and deep inelastic scattering, to take place. In the non-runaway regime, we require In light of this figure, we conclude that some new effects pointed out in our study are also relevant in confining phase transitions where T nuc ∼ T start ∼ T c (see appendix A for the definition of the critical temperature T c ), e.g. [97][98][99][100][101][102][103], provided c vac is small enough. A possible impact on the QCD phase transition, e.g. [104][105][106][107][108][109][110][111][112][113], remains to be investigated.
Averaged quantities only. We conclude this section by also stressing that all the conditions above refer to averaged quantities, and therefore do not take into account the leaks from tails of distributions. These leaks could for example imply that there are a few strings JHEP04(2021)278 that hadronise without ejecting particles, even if all conditions eqs. (6.4), (6.7) and (6.10) are violated. As these strings constitute a small minority of the total ones, these effects have a negligible impact on the phenomenology we discuss. They could however be important in studying other situations of supercooled confinement. Though certainly interesting, the exploration of these effects goes beyond the scope of this paper.

Density of ejected techniquanta
In the wall frame, since we have one ejected quark or gluon per each incoming one, we find n ej,w = n TC,w (r ej ) = γ wp (r ej )n TC,p , where n TC,p = g TC ζ(3)T 3 nuc /π 2 is the density of the diluted bath in the plasma frame. The density of ejected techniquanta then depends on the time passed since bubble wall nucleation, or equivalently on the bubble radius at the time of ejection r ej , via γ wp (r) (see section 5). In the plasma frame, and at a given distance D from the center of the bubble, we then have 7 where we have included the surface dilution from the expansion between the radius at which a given quark has been ejected, r ej , and the radius D where we are evaluating n ej,p .
Radial dependence. It is convenient to express n ej,p as a function of the radial distance x from the bubble wall in the plasma frame, where for definiteness x = 0 denotes the position of the wall and x = L ej,p the position of the techniquanta ejected first (which constitute the outermost layer). In order to do so, we determine the relation between the position x of a quark and the radius r ej (x) when it has been ejected. We assume that the bare mass of the quarks is small enough such that they move at the speed of light, like the gluons. The wall at x = 0, instead, moves at a speed v wall 1 − 1/(2γ 2 wp ) (we have used the relativistic limit γ wp 1), dependent on its radius. The coordinate x of a given layer of ejected particles can then be found by integrating the difference between the world line of an ejected particle and that of the wall, where we defined t D and t ej as the times when the bubble radius is respectively D and r ej , and we used γ wp (t) T nuc t, cf. eq. (5.2), valid up to relative orders 1/γ 2 wp 1. It is convenient to rewrite eq. (7.3) as

JHEP04(2021)278
We finally obtain nuc D x 4 n TC,p , (7.5) where the last equality is valid as long as the bubbles run away, i.e. as long as eq. (5.2) γ wp T nuc r holds.
Thickness of the layer of ejected techniquanta. Our result eq. (7.5) implies that the highest density, of ejected techniquanta, is located in the shell within a distance of the bubble wall L eff ej,p The density of ejected quarks n ej,p (x) extends to x = L ej,p , i.e. to the outermost ejected layer, that we now show to be much larger than L eff ej,p . Indeed, L ej,p can be related to the time t first of ejection of the first techniquanta (corresponding to γ wp m π /T nuc , eq. (5.11)). Using t i t first and t first ∼ m π /T 2 nuc , we find where for simplicity we have assumed m π f as in QCD. As long as L ej,p L eff ej,p , as it holds for our estimate eq. (7.7), the value of L ej,p does not affect any of the results of this paper. 8 The density profile of eq. (7.5) is shown in figure 7.
Sanity check. As a check of our result eq. (7.5), we verify that one has one ejected quark or gluon per each one that entered the bubble. Indeed, we compute where we have assumed D f /T 2 nuc , i.e. we have placed ourselves deep in the regime where hadrons can form inside bubbles (see eq. (5.11)). Equation (7.8) guarantees that the number of ejected techniquanta in the layer of thickness L ej,p is equal to the total number of techniquanta that entered the bubble up to radius D.
Interactions between ejected quarks. Let us finally comment why, we think, interactions among the ejected techniquanta cannot much alter their density. The density of the particles in the incoming bath does not change out of their own interactions. In the wall frame, both the density and the relative momentum of the ejected techniquanta are of the same order of those of the particles in the incoming bath. Therefore, we analogously expect that the density of the ejected techniquanta would also not change after ejection. Since what will matter for the following treatment is the energy in the ejected techniquanta, rather than how this energy is spread among the various degrees of freedom, we content ourselves with this qualitative understanding and leave a more precise treatment to future work.

Scatterings of ejected quarks and gluons before reaching other bubbles
Before possibly reaching other expanding bubble-walls and their ejected techniquanta, ejected quarks and gluons could undergo scatterings with particles from the supercooled bath at temperature T nuc , and with techniquanta ejected from other bubbles. In this section we study the effects of these scatterings.
Ejected techniquanta are energetic. As soon as a bubble occupies an order one fraction of its volume at collision, the total energy in ejected particles is much larger than that in the supercooled bath outside the bubble. Indeed, we have seen that for each quark or gluon in the supercooled bath that enters a bubble, there is at least an ejected one, and that the energy ejected per each incoming particle is much larger than the energy per each particle in the bath, E ej,p γ wp f T nuc , eq. (4.18). Assuming the degrees of freedom in quarks and gluons are not an extremely small fraction of those in the diluted medium, then the diluted medium outside the bubbles does not have enough energy to act as a bath for the ejected particles. This implies that most ejected particles keep most of their energy upon passing through the supercooled bath.
Energy transfer between ejected techniquanta and diluted bath. By reversing the logic above, the ejected particles can deposit in the supercooled bath an energy much larger than its initial one. Pushing this to the extreme, the ejected techniquanta could make the bath move away from the bubble wall, thus making our treatment so far valid only in the first stages of bubble expansion. In order to assess this, we estimate the rate of transferred energy between ejected techniquanta and particles from the bath outside the bubbles,

JHEP04(2021)278
where n ej is the density of ejected techniquanta, ∆E is the energy transferred per single scattering, and where in the second equality we have taken the limit of relativistic particles v 1 and small energy transfer per single scattering ∆E, so that the Mandaelstam variable t can be expressed as t −∆E 2 . The quantity dσ/dt depends on the specific model under consideration, in particular it depends both on whether the ejected particle is a quark or a gluon, and on the identity of the scatterer in the bath outside the bubbles. For definiteness, we model it as the cross section for fermion-fermion scattering mediated by a light vector with some effective coupling √ 4πα eff , We then obtain Γ ej-bath is of course not Lorentz invariant, it depends on the frame via the density of ejected techniquanta n ej determined in section 7.1.
Impact on diluted bath. The average energy transferred to a particle in the diluted bath at position D, when this particle goes across the layer of ejected techniquanta (so before it reaches the wall and initiates the processes described in section 4), then reads where we remind that the spatial coordinate x is the distance between a given layer of ejected techniquanta and the wall at x = 0. Upon use of eqs. (7.11) and (7.5), we can then evaluate the average energy transferred to an incoming particle from the diluted bath, eq. (7.12), as Note that the product D n TC,p is Lorentz-invariant, so that Q ej-bath is indeed a Lorentzinvariant quantity. To learn whether particles from the diluted bath are prevented from entering the wall, because of the interaction with the ejected techniquanta, we compare the energy they exchange with them upon passing their layer with their initial energy in the wall frame, 9 E i,w 3 γ wp T nuc , 14) The novel physical picture we described in sections 4 and 7 is valid as long as Q ej-bath /E i,w 1. As seen in section 5, γ wp initially grows linearly with the bubble radius, γ wp T nuc D, until the retarding pressure possibly becomes effective. It will turn out in section 9 that the runaway regime of linear growth is the one relevant for the phenomenology we will discuss. In that regime, the condition Q ej-bath /E i,w 1 translates into T nuc / √ −t IR 1.

JHEP04(2021)278
IR cut-off. The quantity −t IR is the IR cutoff of the scattering, −t IR ≡ m 2 V , with m V some effective mass of the mediator responsible for the interactions that exchange momentum. In the absence of mass scales, which is the case for example for the SM photon and for the gluons, the effective mass m V is equal to the plasma mass of these particles in the thermal bath. If the only bath was the diluted one, one would have m 2 V,therm ∼ α eff n TC,p / E TC,p ∼ T 2 nuc (see e.g. [114]). However, the process of our interest here happens in the much denser bath of ejected techniquanta, n ej,p n TC,p , so that we indeed expect m 2 V,therm T 2 nuc , so that Q ej-bath /E i,w 1 and our picture so far is valid. More precisely, the screening mass for non-equilibrium systems scales as [115] (f (p) is the non-equilibrium phase space distribution of the particles in the system) where we have used E ej,p ∼ γ wp f and n ej,p ∼ (D/L ej,p )n TC,p ∼ D 2 T 5 nuc ∼ γ wp (D) 2 T 3 nuc . Equations (7.14) and (7.15) teach us that, in the regions of parameter space where γ wp f /T nuc , the energy received by each particle in the diluted bath, from scatterings with the ejected techniquanta, is much smaller than their energy in the wall frame E i,w 3 γ wp T nuc . 10 Since E i,w was the crucial input quantity for our treatment in section 4, the picture that emerged there is not affected by these scatterings.

Energy transferred to techniquanta ejected from other bubbles.
Finally, before ejected techniquanta can possibly enter another expanding bubble, they also have to pass through the layer of the techniquanta ejected from that other bubble. To investigate this, one can use the result derived above, eq. (7.13), with the specification that now D is the maximal radius reached on average by expanding bubbles, because the shells of ejected quarks and gluons meet just before the bubble walls do. We then find that the average energy transferred is much smaller than the energy of an ejected techniquanta in the plasma frame γ wp f , One could be worried that in the outer shell of size Lej,p ∼ 1/f , the thermal mass m V,therm is much smaller than its value in the densest region in x 0, explicited in eq. (7.15), such that Q ej-bath E i,w becomes larger than 1. We can check that it is not the case by including the x-dependence of m V,therm , eq. (7.15), in the integral in eq. (7.12) where nej,p is defined in eq. (7.5), Lej,p ∼ 1/f in eq. (7.7), D ∼ γwp/Tnuc in eq. (5.2), and Eej,p ∼ Eej,p ∼ γwpf . We compute the integral in eq. (7.16) and obtain This confirms that the energy of incoming particles, Ei,w 3γwpTnuc, is not affected by the shell of ejected quarks.

JHEP04(2021)278
Hence, for the purpose of determining the average energy of ejected quarks when they enter another bubble, one can safely ignore the interactions between the two shells.

Ejected techniquanta enter other bubbles (and their pressure on them)
Ejected techniquanta are squeezed. In the plasma frame, all ejected techniquarks are contained within a shell of length given by eq. (7.7) L ej,p ∼ 1/f , and most of them lie within a length given by eq. (7.6) L eff ej,p ∼ 1/(T 2 nuc D) 1/f . In the frame of the wall of the bubble they are about to enter, these lengths are further shrunk, so that ejected techniquarks are closer to each other than 1/f by several orders of magnitude. Therefore we expect no phenomenon of string fragmentation when they enter other bubbles. So each ejected particle, upon entering another bubble, forms a hadron with one or more of its neighbours. This also implies there is no further ejection of other techniquanta. Each of these hadrons carries an energy equal to that of the techniquanta that formed it, of order γ wp f in the plasma frame.
Contribution to the retarding pressure. This conversion of ejected techniquanta into hadrons results in another source of pressure on the bubble walls, that acts for the relatively short time during which the bubble wall swallows the layer of ejected techniquanta. In the frame of the bubble wall that they are entering, the energy of each ejected quark or gluon reads E ej,w2 2γ 2 wp f . We then proceed analogously to what done in section 5.1, and compute (7.20) where we have used g TC = g g + 3g q /4, ∆m 2 in f 2 , n ej,w2 2γ wp n ej,p and, for simplicity, the peak value n ej,p 2 γ 2 wp n TC,p of eq. (7.5). The population of techniquanta ejected from other bubbles thus exert, on a given bubble wall, a pressure comparable to that exerted by the techniquanta incoming from the bath at LO, cf. eq. (5.6). Therefore, the pressure from ejected techniquanta does not alter the picture described so far -a fortiori -because it is exerted only just before bubble walls collide and not throughout their entire expansion.

Ejected techniquanta heat the diluted SM bath
In section 7.2 we found that the scatterings between ejected techniquanta and the diluted bath do not quantitatively change the picture of string fragmentation described in section 4. These scatterings may however affect the properties of the particles, in the diluted bath, that do not confine. These particles include all the SM ones that are not charged under the new confining group, so that for simplicity we denote them as 'SM'. By a derivation analogous to the one that lead us to eq. (7.13), we find that the average energy they exchange with the ejected quarks reads

JHEP04(2021)278
where we have used T nuc D γ wp and −t IR ∼ T 3 nuc γ wp /f , cf. eq. (7.15). We have denoted by α SM an effective coupling between SM particles and the techniquanta, which is modeldependent. Now assume the techniquarks carry SM charges, e.g. as expected in composite Higgs models. Then, in the wall frame, the fractional change of energy is of course similar to that derived in eq. (7.14) for the incoming techniquanta. However the incoming techniquanta next undergo string fragmentation, and eq. (7.14) does not affect that energy balance for γ wp f /T nuc . In other words, string fragmentation renders this energy transfer irrelevant for the techniquanta, while the SM particles neutral under the confining group just proceed undisturbed so they keep track of it. In particular, Q ej-SM , is much larger than the latter energy in the plasma frame ∼ T nuc , and may even be slightly larger than the confinement scale f . 11 This need not be the case, however, as the new techniquanta may be very weakly interacting with the SM. As they cannot interact too weakly, otherwise our assumption of instantaneous reheating would not hold, for simplicity we ignore this case in what follows and we assume that some techniquarks carry SM charges.

Deep inelastic scattering in the early universe
The physical picture described so far results in a universe that, before (p)reheating from bubble wall collisions, contains three populations of particles.
• Population A. Arises from hadronisation following string fragmentation. It consists of N string ψ /2 hadrons per quark or gluon in the initial bath, each on average with energy in the plasma frame, and of roughly the same number of hadrons with much smaller energy. (The latter can be thought as coming from the half of the string closer to the center of the bubble wall.) The physics resulting in this population is described in section 4, see eq. (4.17) for E A and eq. (4.14) for N string ψ (E CM = 3 γ wp T nuc f ).
• Population B. Comes from the hadronisation of the ejected techniquanta. This population consists in ∼ one hadron per quark or gluon in the initial bath, each with energy So this population carries an energy of the same order of that of population A. Its physics is described in section 7, the energy E B is that of the initial quark or gluon, eq. (4.18). 11 As already anticipated, in the regime of interest for DM phenomenology we will find that bubble walls run away, so that γ max wp is (much) smaller than ∼ 10 −3 (f /Tnuc) 3 , see eq. (5.1). Note also that, for eq. (7.21) only, gTC = 3gq/4, i.e. the gluon contribution to heating the SM is negligible because they cannot carry SM charge.

JHEP04(2021)278
• Population C. Consists of the particles that do not feel the confinement force, that we denote 'SM' for simplicity, each with a model-dependent energy given by eq. (7.21), and whose total energy is much smaller than that in populations A and B.
The direction of motion of all these populations points, on average, out of the centers of bubble nucleation.
Hadrons from both populations A and B have large enough energies, in the plasma frame, that showers of the new confining sector are induced when they (or their decay products) scatter with the other particles in the universe and/or among themselves. These deep inelastic scatterings (DIS): • Increase the number density of composite states.
• Decrease the momentum of each of these states with respect to the initial one | p ψ |.
Hence, such effects need to be taken into account to find the yield of any long-lived hadron.
The evolution of our physical system would require solving Boltzmann equations for the creation and dynamics of populations A, B and C in a universe in which preheating is occurring, and of the interactions of populations A, B and C among themselves and with the preheated particles produced from bubble wall collisions. While certainly interesting, such a refined treatment goes beyond the purpose of this paper. In this section, we aim rather at a simplified yet physical treatment, in order to obtain an order-of-magnitude prediction for the yield of long-lived hadrons.

Scatterings before (p)reheating
We begin by considering the interactions among populations A, B and C.
Number densities of scatterers. Let us define L X , with X = A, B, C, the effective thickness of the shells containing populations A, B, and C respectively. For example, L B,p = L eff ej,p of eq. (7.6). We know that population A(B) consists on average of N string ψ /2 hadrons (one hadron) per each quark or gluon in the initial diluted bath, and that population C is the initial diluted SM population. By conservation of the number of particles, we then obtain the number densities where D is the average radius of a bubble at collision and we have used L X D.
Energy transferred between scatterers. We now determine the average momentum, transferred to a particle from population X, upon going across a shell of population Y . In order to do so, we use our result eq. (7.11) for the rate of transferred energy and compute where α X-Y is the effective interaction strength of the scatterings of interest, g Y the number of degrees of freedom in density of population Y (where we include a factor of N string ψ /2

JHEP04(2021)278
for Y = A), and we have used the relation T nuc D = γ wp valid in the runaway regime. We conclude that: • Populations A and B. The energies of the hadrons of population A and B in the plasma frame, respectively γ wp f /N string ψ and γ wp f , are both much larger than the energy they can exchange with any of the other baths among A,B,C, by a factor that scales parametrically as f /T nuc or larger (because for all populations we have −t IR ∼ n/ E > T 2 nuc , see the discussion in section 7.2). Therefore these elastic scatterings are not effective in reducing the energy of the hadrons of either population A or population B.
• Population C. On the contrary, Q A,B→C can be of the same order of the energy of each particle in population C, eq. (7.21), which therefore are significantly slowed down by these interactions. Importantly for our treatment, this does not alter the fact that population C was energetically subdominant with respect to populations A and B.

No significant DIS between populations A, B and C.
Finally, we determine whether any of the scatterings among particles in populations A,B,C could result in significant hadron production, via deep inelastic scattering. A single scattering event potentially results in a shower of the new confining sector if the exchanged momentum is larger than the confinement scale, t 2 > f 2 . This condition is allowed by kinematics, because the centerof-mass energy of the scatterings between any of the populations above is much larger than f . A significant amount of DIS happens if the DIS scattering rate Γ DIS Y→X of a particle from population X, upon going across a shell of population Y , is much larger than the inverse of the length of the shell Y . We then compute where again we have used the runaway relation T nuc D = γ wp and, for definiteness, we have assumed the scattering cross section has the form of eq. (7.10). Therefore, no significant DIS happens in the regions where γ wp (f /T nuc ) 2 . This condition will turn out to be always satisfied in the parameter space of our interest, so we can ignore the DIS among populations A, B and C in what follows.

Scatterings with the (p)reheated bath
By preheating, we intend the stage between the time when bubble walls collide and start to produce particles (e.g. from the resulting profile of the condensate), and the reheating time when these particles have thermalised into a bath. We now discuss the scatterings of populations A and B with the particles produced at preheating, that we have assumed to be efficient. The contribution of population C to the final yield of hadrons is subdominant with respect to the one of populations A and B because, as seen in sections 7.4 and 8.1, the total energy in population C is much smaller than that in populations A and B.

JHEP04(2021)278
Energy of the (p)reheated bath. The preheated particles are produced with energies, in the plasma frame, of the order of the mass of the scalar condensate, 12 Their total energy scales as with V the volume of a large enough region of the universe. For comparison, the total energy in populations A and B scales as which is much smaller than f 4 V because γ wp (f /T nuc ) 3 , eq. (5.1). So the preheated particles can act as a thermal bath for all the other populations A, B and C, because the energy of A, B, and C is subdominant in the energy budget of the universe.
Inelastic versus elastic scattering. Scatterings of hadrons (or their decay products) with the preheated bath will, therefore, eventually slow down and thermalise populations A and B. However, these scatterings can also exchange energies much larger than f , thus inducing deep inelastic scatterings. Indeed their center-of-mass energy squared reads where ) is the result of our simplifying assumption to neglect masses and to average to zero scattering angles with particles in a bath: define p E = E(1,Ê), p prh = m χ (1,m), then s = (p E + p prh ) 2 2E m χ (1 −Ê ·m) 2E m χ . We now determine if those center of mass energies are entirely available for particle production via DIS, or if instead they are reduced by several low-momentum-exchange interactions. In order to do so, we evaluate the rate of energy loss of a particle from population A or B, Γ loss A,B , as the ratio between the rate of energy it exchanges with the preheated bath, that we evaluate analogously to eq. (7.11), and its initial energy E A,B . We then compare this quantity with the rate for a deep inelastic scattering to happen with the full energy available s In the last equality, we have again used the screening mass for non-equilibrium systems [115] − t IR ∼ n prh E prh ∼ c vac f 4 m 2 χ , (8.11) 12 In the picture we have in mind, non-perturbative effects such as Bose enhancement or parametric resonance (see e.g. [116]) are not relevant: the first because the SM particles are interacting, thus they exchange momentum and do not occupy the same phase space cells; the second because the variation of their masses from the dilaton's oscillations is smaller than their mass at the minimum. Note that, unlike what occurs in many inflationary scenarios, we expect only a small hierarchy T RH E prh .

JHEP04(2021)278
where we have used that by conservation of energy n prh ∼ ρ RH / E prh , and where we have expressed the energy density of the reheated bath ρ RH using the results of section 3.2. We conclude that, if m χ f 2 c 1/2 vac , (8.12) the full center-of-mass energies s A,B are available for deep inelastic scattering, i.e. populations A and B do not lose a significant amount of their energy via interactions with the preheated bath. For simplicity, in what follows we assume this model-dependent property to hold.

Enhancement of hadron abundance via DIS
The picture: a cascade of DIS. The number of composite states arising from a hard scattering depends on how the strings fragment, so on the same physics that set the abundance of the composite states when the techniquanta cross the bubble walls, discussed in section 4.4. Each scattering, depending on its center-of-mass energy, produces a number N string ψ of hadrons ψ, that we model in the same was as in eq. (4.14). Given the large initial energies s A,B , the daughter hadrons typically still have enough energy to themselves induce further deep inelastic scatterings with the particles in the preheated bath, and hence additional hadron production. Analogously, SM particles produced in such DIS typically have large enough energies to also initiate showers of the new confining force with their subsequent scatterings. This process iterates until the average energy of scatterings drops below the confinement scale.
Number of hadrons produced per scattering. For reasons given in section 4.4, together with simplicity, we assume that the available energy √ s at each scattering splits equally among all the outcoming particles. We then write the average of this number as where the factor of 2 in the argument of N string ψ arises because eq. (4.14), which defines N string ψ , assumes that √ s is the center of mass energy of the scattering of two particles neutral under the new confining force. If a hadron is included among the two scatterers, then QCD studies find that the final number of hadrons can be obtained by just halving the energy in the center of mass frame [69], also see footnote 5. 13 Energies of produced hadrons. Explicitly, we assume E com = √ s/N DIS , where E com is the energy of any outgoing particle (SM and/or composite) in the center-of-mass frame of the scattering. To iterate to many scatterings, we write E com in the plasma frame as 13 Note that if a hadron instead decays to two SM particles before it scatters, which is model-dependent, then √ s/2 is again the good argument for the function N string ψ , because then one has two particles each with half the initial energy, but both neutral under the new confining force. In this case, however, eq. (8.13) becomes N DIS (s) = 2N string ψ ( √ s/2). When iterating the treatment to many scatterings, we find that this extra factor of 2 does not impact the final abundance of hadrons, which can be understood by thinking that the same initial energy is spread faster to zero.

JHEP04(2021)278
where γ andv are the associated Lorentz boost and its direction, andv is the direction of motion of the outgoing particle in the center-of-mass frame of the scattering. By averagingv ·v to zero for simplicity, we obtain (8.14) We then determine γ by observing that the energy of each particle, in the center-of-mass frame of the scattering, is both E com = √ s/2 and E com = γ E prh (1 +v ·Ê com ), where E prh is the energy in the plasma frame of the particles in the preheated bath. By averaginĝ v ·Ê com to zero for simplicity, we obtain the Lorentz boost Using eq. (8.9) for s we finally obtain (If we did not average over angles, we would have obtained . So, after a hard scattering the energy of each outgoing particle in the plasma frame is roughly the initial energy divided by a factor N DIS . The subsequent s is then reduced by the same factor, ensuring a convergence of N DIS (s) to unity, via eq. (4.14), after only a few iterations. This also teaches us that the average energy of the particles, produced this way, quickly decreases to values lower than about m * .

Number of hadrons produced by a chain of DIS.
Let us now estimate the yield of final hadrons by following the above arguments. Assuming interactions are fast enough, also those following the first one happen with preheated particles of the same average energy E prh . Now define the number of states (both composite and not) N k produced at the k th interaction. This can be expressed as where we remind the reader that the function N DIS is obtained from eqs. (4.13) and (8.13).
Starting from a single resonance produced from the fragmentation of strings between quanta inside the bubble, after this chain of scattering processes one obtains a total number of resonances given by the product k N k (s). We find numerically that this product can be expressed as In other words, the iterative process we described converts the initial available energy into the rest mass of hadrons m * . Since our aim here is not to achieve a more precise treatment, we refrain from refining the assumption that the momenta are distributed evenly among the particles coming out of a scattering process. In the same spirit of building a physically-clear picture without drowning in model-dependent details, we do not cover here the possibility

JHEP04(2021)278
that every scattering produces, in addition to the composite states, a comparable or larger amount of SM particles. (That would result in N DIS > N string ψ and in a faster degrowth of the available scattering energy to m * at each step.) In addition to simplicity, this can be justified by observing that, in the limit of large number of degrees of freedom in the dark sector, our assumption that they carry SM charges will make their production dominant with respect to the one of SM particles.
Additional comments. We conclude our derivations with two comments concerning its validity. • We have ignored the production of heavy particles from the collisions of bubble walls [36,[117][118][119]. This is justified as it has been shown that it only occurs when the minima of the potential are nearly degenerate and seperated by a sizable barrier [120,121], which is not the case for the close-to-conformal potentials we have in mind.
Hence we expect only particles lighter than the scalar condensate to be produced during reheating following the wall collision.

DIS summary
The yield of hadrons, resulting from the processes of deep inelastic scattering described above, receives contributions from: • Population A. That is, the hadrons produced from string fragmentation as described in section 4. Their contribution reads where we have used K DIS A = s A /m 2 * , cf. eq. (8.18), with s A 2m χ E A 2m χ γ wp f /N string ψ (E CM ), cf. eqs. (8.9) and (8.1). Note that the above expression captures also the regime where each string fragmentation produces on average one hadron, because the energy of that single hadron is roughly γ wp f /2, see the related discussion in section 4.5.
• Population B. That is, the hadrons produced out of the techniquanta ejected from the bubbles, described in section 7. Their contribution reads

JHEP04(2021)278
Thus, the combined contribution to the total hadron yield is given by where we have defined Note finally that, in the regime of runaway bubble-walls, one obtains the parametric scaling Y SC+string+DIS ∝ (T nuc /f ) 3 γ wp . Which is much larger than the simple supercooling dilution, ∼ (T nuc /f ) 3 , in the regions of parameter space where our analysis holds, namely for γ wp > f /T nuc .
9 Supercooled composite dark matter

Initial condition for thermal evolution
Finally, all unstable resonances decay either to SM or to the long-lived or stable hadrons, which we take to form DM. To obtain the yield of any such hadron i at the onset of reheating, one should use the expression where K DIS , D SC and Y eq TC are defined, respectively, in eqs. (8.22), (3.9) and (3.7). BR i is a pseudo-branching ratio, of the energy available to the confining techniquanta, into ψ i particles. Estimates of BR i for the cases where ψ i is a meson and a baryon are given in appendix C, which show a broad range of underlying-model dependent values are possible, albeit with a large uncertainty. For example, in a QCD-like theory where ψ i is a baryon with mass ∼ 4πf and the pions have mass ∼ f , one obtains values BR i ∼ 10 −6 , while larger values BR i are obtained for baryon-pion mass rations closer to one, or if ψ i is a meson. Hence, we will take BR i to be a free parameter.
For completeness, the supercooling plus string and supercooling yields read where K string is defined in eq. (4.15). We have included a factor 3 4 g q /g TC in Y SC i to account for the fact that, in the case of no string fragmentation nor DIS, gluons do not contribute to the final abundance of heavy composite states of quarks. It would be absent if one was interested in light composite states of quarks. The yield of the various contributions is shown in figure 8.
It will turn out that the measured DM abundance is achieved in the regime of runaway bubble walls. In that regime, the resulting expression for the DM yield has a simple parametric form that eventually results in the DM abundance being independent of the DM mass, if it is to match onto observation Y DM 0.43 eV/m DM [122], which we find convenient to report here. By using eqs.

Thermal contribution
To complete our discussion, we must still determine the effects on the yield of any DM interactions with the thermal bath after supercooling, DIS, and reheating. The importance of thermal effects following reheating was already pointed out in [20] (therein dubbed the subthermal contribution). Following the phase transition and particle production through DIS, the SM bath and the DM have returned to kinetic equilibrium. The scattering energy is now insufficient to break the resonances, but these may still annihilate into SM particles or be produced in the inverse process. Thus, just after the reheating, the DM abundance evolves according to the well known Boltzmann equation [123] dY where we use x ≡ m DM /T as the time variable, and M pl is the reduced Planck mass. For simplicity we only consider velocity independent cross sections here. As an intitial condition we take the relic abundance at the reheat temperature, Y DM (T RH ) = Y SC+string+DIS DM , estimated following string fragmentation and DIS enhancement in eq. (8.21). For our plots we solve the Boltzmann equation numerically. If the cross section and reheating temperatures are sufficiently large the system will be driven back into equilibrium. The relic density is then largely set by freezeout dymanics, albeit with somewhat different initial conditions. On the other hand, if the cross section and reheat temperatures are small enough, the relic density is set by dilution, string fragmentation and DIS, with only negligible thermal corrections following reheating. Using the dilution mechanism of the PT, of course, we JHEP04(2021)278 Below this mass, the relic density is necessarily suppressed compared to the observed DM density, due to efficient DM annihilation after reheating. In the purple region 3 γ wp T nuc f the quarks are reflected by the first wall they encounter, but may enter the bubbles in following stages of their evolution, and the DM abundance lines ignore possible modifications arising from this 'ping-pong' effect. They also ignore that, for values of γ wp only slightly larger than f /T nuc and depending on other model-dependent parameters, the energetics of our treatment may be more complicated, see eqs. (7.14) and (7.15). The dashed gray line delimits the area T nuc < O (100) MeV where the supercooled phase transition could happen because of QCD dynamics. The dashed light blue line indicates the regimes where bubble walls run away, cf. eq. (5.10). The dashed purple line indicates the regime where γ wp < (f /T nuc ) 2 , and the fact it lies above the horizontal part of the DIS line confirms that our treatment has been consistent when ignoring the DIS of eq. (8.5).
can avoid the usual unitarity constraint on the maximum thermal relic DM mass [46] (see e.g. [124,125] for recent appraisals).

Dark matter relic abundance
We now combine all our results together and determine the amount of supercooling required to match the observed relic abundance Y DM 0.43 eV/m DM . Examples are shown in figure 9 for some representative choices of the parameters. From these figures we can draw a number of conclusions.

JHEP04(2021)278
i) If we assume σv rel ∝ 1/m 2 DM , thermal effects will necessarily dominate if the DM is light enough. This occurs because T RH cannot realistically be arbitrarily suppressed below f , for sensible choices of g Ri and g Rf . This regime corresponds to the point in which the contours turn vertical in figure 9. At which value of m DM this occurs depends on the precise choice for σv rel . For definiteness, in figure 9 we choose σv rel = 4π/m 2 DM as typical of baryon scatterings in a strongly coupled sector. Thermal effects can of course be further suppressed if we depart from the efficient reheating assumption made here [20].
ii) String fragmentation and DIS lead to large corrections to the composite DM relic density, compared to the naive supercooling dilution. This implies a mismatch between the relic abundances of primordial elementary and composite relics alluded to before. Whether the composite or elementary relic would have the greater abundance depends on the details of confinement (for elementary relics BR i = 1). If the composite relic is say, a light meson which is produced abundantly, the multiplicative DIS process can be highly efficient in populating these states following the PT. This implies we require much more supercooling to match onto the observed DM relic abundance. On the other hand, if the composite relic is some heavy state, perhaps a baryon, it could be produced in a highly suppressed rate both in string fragmentation and DIS. In this latter case, the required amount of supercooling to match onto the DM relic density is also reduced. The two cases are illustrated with two different assumptions for the branching ratios in figure 9. 14 iii) In some cases, we find T nuc 100 MeV, as delineated in figure 9. Thus QCD effects could assist in completing the PT [20,37,127,128]. On the other hand, if QCD effects help the transition to occur, they can also suppress the eventual gravitational wave signature [91] (simply because the QCD effects increase the tunneling probability and thus will act to shorten the timescale of the PT). The details will depend on the physics entering the effective potential of the scalar χ and need to be studied in a model dependent way.
Together with the gravitational wave signal from the PT, there may also be model dependent collider, direct, and indirect detection signatures associated with the DM from the strongly coupled sector. We will investigate these further, together with their interplay with the novel string fragmentation and DIS effect, in a concrete realisation of such a confining sector in a companion paper [23]. 14 We checked that the DIS line is unaffected if we use the treatment of the NLO pressure of [92], instead of the one of [72] that we have employed in this paper, cf. section 5.2 and appendix B. On the other hand, the run-away wall dashed line, and hence the string fragmentation line, could be affected by this choice. For simplicity as well as in light of the criticism of [92] appeared in [93,126], we employed the results of [72] in this paper.

Discussion and outlook
The possible existence of a new confining sector of Nature is motivated by several independent problems of the Standard Model of particle physics and by cosmology. This encourages the identification of predictions of confining sectors, that are independent of the specific problem they solve. One such prediction is the possibility that the finite temperature phase transition in the early universe, between the deconfined and confined phase, is supercooled. This possibility has received a lot of attention in recent years, see e.g. [20, 27-31, 37, 75, 91, 128].
In this paper, we have pointed out and modelled a novel dynamical picture taking place in every supercooled confining phase transition, that (to our knowledge) had been missed in the literature. This novel picture stems from the observation that, when fundamental techniquanta of the confining sector are swept into expanding bubbles of the new confining phase, the distance between them is large with respect to the confinement scale. Therefore the energy of the fluxtubes connecting techniquanta is so large that string breaking produces many hadrons per fluxtube, with large momenta in the plasma (CMB) frame, in a sense analogously to QCD hadrons produced in electron-positron collisions at colliders. These hadrons and their decay products subsequently undergo scatterings with other particles in the universe, with center-of-mass energies much larger than both the confinement scale and the temperature that the universe reaches after reheating. The dynamics just described is partly pictured in figures 2 and 3.
The processes of string fragmentation and 'deep inelastic scatterings in the sky', synthetised above, have a plethora of implications. A key quantity to study them is the pressure on the bubble walls induced by this novel dynamics, which we have determined in section 5, see eq. (5.1) and figure 5 for the resulting bubble-wall velocities. An interesting aspect of our findings is that the so-called 'leading-order' pressure is proportional to the boost factor of the bubble wall, unlike in the case of non-confining supercooled PTs [71,72].
We then quantified the values of supercooling below which one recovers the 'standard phase transition', where confinement happens between nearest charges. By relying on the modelling we proposed in section 6 we found, interestingly, that the PT does not proceed in the 'standard' way already for minor supercooling, i.e. if bubbles are nucleated and expand just after vacuum energy starts to dominate. Our proposed dynamics should not only be employed in the large supercooling region, but also in the minor supercooling one depending on the value of another model-dependent parameter, see figure 6. The regimes in between these regions (one being the 'ping-pong' regime of section 5.3) will be studied in future work, to not charge this paper with too much content.
Next, we have focussed on the implications of our dynamical picture for the abundance of long-lived or stable particles that are composite states of the new confining sector. They are summarised in the Synopsis, section 2, and a quantitatively accurate expression of the final yield of a given composite particle is given in eq. (9.4), for concreteness in the regime where bubble walls run away. Compared to the simple dilution of relics induced by supercooling of non-confinement transitions, these processes enhance their abundance by parametrically large factors. Therefore they have to be taken into account whenever a JHEP04(2021)278 property of the universe, e.g. the DM and/or the baryon abundance, depends on the final yield of hadrons. As an example, their dramatic impact on the abundance of supercooled composite DM can be seen in figure 9.
Concerning DM in particular, this study constitutes a novel production mechanism of DM with mass beyond the unitarity bound [46]. It would be interesting and timely to study its experimental signals, given the new wave of telescopes that is starting to take data of high-energy neutrinos and gamma rays (e.g. KM3NeT, LHAASO, CTA) and given their potential in testing heavy DM, e.g. see [129]. One such study will appear in a forthcoming publication [23].
During the course of carrying out this study we have made a number of simplifications, for the purpose of obtaining a general and clear enough picture of the physics involved. For example, the various populations of particles created by this novel dynamics, such as the ejected techniquanta and the hadrons that follow the bubble walls, could be better described by Boltzmann equations, by the use of simulations etc., rather than with our simple treatment that focused on their average properties.
Finally, this dynamics opens broader and exciting avenues of investigation, that we think deserve exploration. For example, it would be interesting to study its interplay with recent interesting ideas regarding phase transitions [11,16,17,30,32,36,120,[130][131][132][133][134], or its impact on the production of gravitational waves in supercooled confining phase transitions. As for the latter, our study of the bubble wall Lorentz factor in section 5 constitutes a necessary first step. spontaneous breaking of scale invariance generates a pseudo Nambu-Goldstone boson, the dilaton which we parameterize as [135] where f is the confining scale and where σ(x) transforms non-linearly σ(x) → σ(λx) + log λ under the scale transformation x → λx. Its potential is given by [27] V T=0 where m χ is the dilaton mass, and c χ is a constant of order 1, which we fix to c χ = 1. The dilaton coupling constant g χ is chosen to reproduce the glueball normalization with N being the rank of the confining gauge group. The validity of the EFT relies on the smallness of the parameter |γ | 1 (here taken negative) which controls the size of the explicit breaking of scale invariance, and thus of the dilaton mass.
Note that in the limit where |γ | 1, the dilaton potential at zero-temperature reduces to the Coleman-Weinberg potential [136], i.e.
Thermal corrections. To model thermal effects, we follow [17,27], and consider the finite-temperature corrections generated by the particles charged under the confining force (the CFT bosons) The total number of CFT bosons n is fixed to 15 in order to recover the free energy of N = 4 SU(N ) large N super-YM dual to an AdS-Schwarzschild space-time [16] By writing eq. (A.8), we have neglected the contribution from the fermions present in the plasma. For simplicity, we suppose that the dilaton degree of freedom χ still exists in the deconfined phase, such that the total potential for the dilaton is where V χ (χ) and V T (χ, T ) are given by eq. (A.2) and eq. (A.6). We plot the potential in figure 10. The supercooling stage starts when the energy density becomes vacuum- Here the damping is purely geometrical, reminiscent of the SO(3, 1) symmetry and we do not consider the damping due to the dilaton decay or due to the interaction with the plasma (see e.g. [139]). In right panel of figure 11, we display the scalar field JHEP04(2021)278 profile obtained after integration of the time-like equation in eq. (A.24), using the initial condition χ(s = 0) = χ * given by the bounce solution in eq. (A.15).
The full bubble wall profile. The full bubble wall profile is obtained after matching the profile in the space-like region, left panel of figure 11, with the profile in the time-like region, right panel of figure 11. One can see that the first confining scale, that the incoming techniquanta are subject to upon entering the wall, is the exit scale Our explicit computation also shows that the length of the section of the bubble wall where χ = χ * , in the wall frame, satisfies as we assumed in eq. (4.2) in the main text. Then, χ transits to its zero-temperature value f over a length, in the wall frame, of order f −1 .

B NLO pressure on the bubble walls
Transition splitting. In section 5.1, we have presented the retarding pressure due to the change in inertia of the system incoming-quark + gluon-flux-attached-to-the-wall when entering inside the confined phase, as well as the retarding pressure due to the ejected quark. In this section we introduce a possible correction which arises in presence of a finite gauge coupling constant. The correction term, which is called NLO pressure, arises from the possibility for the incoming particle to radiate a soft boson which gets a mass in the broken phase [72] We have introduced the variable x ≡ E c /E a and assumed x 1 in the last equality. We can now split the integral over z across the wall in eq. (B.5) into a contribution from the broken phase and a contribution from the symmetric phase. We assume that the vertices V and WKB phases A on each side of the wall are z-independent and we denote them by (V h , A h ) and (V s , A s ), such that we obtain Radiation of a soft transverse boson. It can be shown [72] that the most important process contributing to the pressure at large E a is likely to be X(p a ) → V T (p c ) X(p b ) where V T is a transverse vector boson. The corresponding vertex function is phase-independent, V h = V s , and equal to where g is the gauge coupling constant and C 2 [R] is the second casimir of the representation R of the incoming particle with respect to the gauge group. Therefore, eq. (B.12) becomes (B.14) where we have replaced m V ≡ m c . The k ⊥ integral becomes where k * is the IR cut-off on k ⊥ . It is expected to be of order of the vector mass k * ∼ m V .
Final NLO pressure. Finally, injecting the last two equations into eq. (B.8) yields . (B.16) The Pauli blocking or Bose enhancing factor 1 ± f b is of order 1, while 1 ± f c sums to 1 when considering both absorption and emission processes. Hence, the result simplifies to where b a = 1 (3/4) for bosons (fermions) and α ≡ g 2 /4π. The Lorentz factor γ wp between the wall and the plasma comes from d 3 p a . We have introduced ps ≤ 1 to encode the suppression from phase-space saturation of the emitted soft techni-gluon,which is important for large coupling g, and which we justify in the next paragraph.

JHEP04(2021)278
Phase-space saturation. At order g 4 , the emitted gauge boson can interact among each other. These processes are weighted by g 2 f (k) with respect to NLO case studied in the last section. The occupancy function f (k) can be estimated to be of order where ∆k ∼ m 3 V is the available phase space and n ∼ P 1→2 /(p z a − p z b − p z c ) with (p z a − p z b − p z c ) ∼ m V . Hence, we can not consider the individual transition splitting processes as independent from each other as soon as At such large γ wp , we expect the NLO pressure to change behavior. See [72] for more details and particularly about some hints of P NLO going from ∝ γ wp to γ 4/7 wp . For simplicity, we just encode this effect into the coefficient ps ≤ 1 in eq. (B.17).
Case of a SU(N ) confining sector. In the scenario we are interested, the deconfined phase contains g q techni-quark and g g techni-gluons, and the NLO pressure would be induced by the possibility for these techni-quanta, to radiate a soft techni-gluon acquiring a mass m V = m g in the confined phase. Hence, eq. (B.17) becomes where g conf is the gauge coupling of the confining group, and where C 2 [g] = N , C 2 [q] = (N 2 −1)/2N if the confining gauge group is SU(N ). Note that in the parameter space which we consider (c vac = 0.01, g TC = 78) the LO pressure in eq. (5.8) prevents the condition in eq. (B.19) to be satisfied such that we expect ps to be close to unity.

C Example estimates of the string to DM branching ratio
In section 4.2, we have discussed that, after supercooling, the quarks enter inside the confined phase, with a typical seperation ∼ T −1 nuc , much larger than the confining scale f , such that a highly energetic fluxtube forms. We have shown that this string, which is unstable under quark-anti-quark pair nucleation, breaks into K string pieces. The dynamics of strings is then also relevant in the processes of deep inelastic scatterings of section 8. In this section, we estimate the branching ratio of a string to a given hadron i, introduced in eq. (9.1), in two different cases. First, when i is a light meson, in which case we expect the yield of i to be independent of its mass and given by a combinatoric factor implying the number of flavors. Second, when i is a heavy baryon in which case one expects the yield to be Boltzmann suppressed.

JHEP04(2021)278
Light meson -combinatorics. In the limit of large string energy, E CM f , one expects the fragmentation of the string to be democratic with respect to the different bound-states if they are light enough. In that case, the string-to-i branching ratio is given by a combinatoric factor depending on the number of flavors N f and the number of quark constituents (either 2 for meson and N TC for baryons). In the particular case of a light meson q 1q2 , one obtains Heavy baryon -Boltzmann suppression. For this example a useful model for us will be the thermal one [140][141][142][143], which was able to fit LEP data of particle yields up to a 10% error [144], even with an initial state far from thermal equilibrium. In this model, the yield of heavy mesonic or baryonic resonances is suppressed by a Boltzmann factor [140][141][142][143], in which the strong scale plays the usual role of temperature. The yield of heavy resonances can be modelled by where M i and J i are the mass and spin of the state i respectively. Here A i is an overall normalisation, which will depend on whether the particle is a pseudoscalar meson, vector meson, or baryon etc. In QCD it was found to differ by 10 between vector mesons, tensor mesons, and baryons [143]. For these particles B i was found to be a common factor between the groups, B i ≡ B ∼ 150 MeV [143]. Note the pseudoscalar mesons in QCD, however, which are lighter, follow a softer spectrum.
Following the above discussion, we shall construct a toy model for the baryonic particle yield from our string fragmentation. In order to retain some simplicity in our model we will consider all particles to share a common B i = m * = g * f . In our toy model we consider SU(N c ) theories, with techniquarks in the fundamental representation, in which baryons will contain N c quarks. Mesons on the other hand will contain a quark-antiquark pair independent of N c . In order to take into account the reduced probability of creating a baryon as opposed to a meson it is therefore suitable to include an additional suppression in the prefactor A i for baryons [145] Other than this we take a common A i = p Bi A. Applying energy conservation, we thus find the average number of the composite state i produced per string breaking to be where the sum runs over all the states in the spectrum, and we remind that π denotes the lightest composite state(s). In this case it is clearly possible to have a highly suppressed BR i , e.g. BR i 10 −6 for m i = m * 4πf , m π f , N c = 10, N π = 3.

JHEP04(2021)278
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.